The impact of molecular dynamics on drug design: applications for the characterization of ligand–macromolecule complexes

Highlights

• Which molecular dynamics (MD) techniques are available for drug design?

• How is MD used to investigate ligand–macromolecule complexes?

• How are MD studies applied to human and non-human therapeutic targets?

• What are the latest advances in the field of MD?

Among all tools available to design new drugs, molecular dynamics (MD) simulations have become an essential technique. Initially developed to investigate molecular models with a limited number of atoms, computers now enable investigations of large macromolecular systems with a simulation time reaching the microsecond range. The reviewed articles cover four years of research to give an overview on the actual impact of MD on the current medicinal chemistry landscape with a particular emphasis on studies of ligand–protein interactions. With a special focus on studies combining computational approaches with data gained from other techniques, this review shows how deeply embedded MD simulations are in drug design strategies and articulates what the future of this technique could be.

Introduction

Molecular dynamics (MD) has become a new major technique in the arsenal of tools developed to design novel bioactive molecules and investigate their mode of action. The main experimental technique to elucidate the molecular structures of macromolecules is the X-ray diffraction of crystallized protein. This technique has led to great achievements in the field of structural biology; but, in the best cases, X-ray crystallography can only provide a static snapshot of a fully functional state. In the past decade, NMR has become an increasingly important technique for protein structural investigations, giving access to the flexibility of a system by revealing an ensemble of conformations. But, despite giant steps achieved lately in the field, NMR spectroscopy remains challenging and time consuming for large protein complexes. Consequently, complementary tools have always been required to allow a dynamic insight into biological targets and ligand binding. MD is an in silico method that can utilize structural data gained experimentally to extrapolate the possible conformations of molecular systems and the different paths between each of them.

Computationally driven simulations of constrained molecular systems were initially used to investigate models with a limited number of atoms. But an impressively fast technological development made possible an astonishing increase of computational speed and data storage volume over the past decade. Today, investigations of larger systems, such as receptors fully incorporated within their biological environment, are routinely conducted. The simulation time being constantly extended, current studies report simulations in the microsecond range, enabling more relevant events to be observed in silico for the rationalization of experimental data.

The scientific literature published in the past four years led to a surprisingly high amount of studies reporting MD simulations (>750 articles – using the very restrictive keywords ‘drug’ and ‘molecular dynamics’). If this analysis of the current situation acknowledges that MD has become a major technique, it remains unclear how MD simulations are used and what unique information they give access to. To answer these questions, we have reviewed and discussed peer-reviewed articles published since 2011 in various fields of medicinal chemistry. In the next section, a short overview on the tools and methods available in the field of MD is given. Then, studies involving MD at a key step of the investigation of therapeutically relevant targets with a ligand are reported. The reviewed strategies are discussed and compared for each main class of enzyme structure investigated in the literature, including separated sections for nucleic acids and targets from pathogenic microorganisms. We finally close this review with a discussion on the weaknesses of the technique and exceptional studies focusing on the improvement of these specific points. This review aims at giving an overall picture of the actual impact of MD on the current landscape in medicinal chemistry and to anticipate what the future of this technique could be.

The most established computational tools available for classical MD simulations are the software suites GROMACS [1], NAMD [2], CHARMM [3] and AMBER [4] (Table 1). MD simulations are a numerical solution of Newton's equations of motion (Eq. (1)):

$${\text{f}}_i=m_i{\text{a}}_i=-\frac{\partial V({\text{r}}^N)}{\partial {\text{r}}_i}$$

where m_i is the mass of atom i, a_i the acceleration of atom i and f_i is the force acting on atom i given by the partial, spatial derivative of the potential energy function V which is dependent on the positions r^N = (r₁, r₂, …, r_N) of all N particles in the system. MD simulations generate a time series of configurations of a system, which allows studying dynamic processes. The time series can also be viewed as a statistical-mechanical ensemble of configurations. By averaging over this ensemble, thermodynamic properties of a system can be accessed.

Four basic elements define the model that is simulated by MD: (i) resolution of the model; (ii) description of the interactions between atoms; (iii) generation of configurations; and (iv) definition of boundary conditions. The resolution defines the ‘elementary particles’ of the model [5]. These can be nuclei and electrons in quantum-mechanical (QM) simulations, atoms in classical MD simulations or coarse-grained (CG) beads where multiple atoms are described by a single particle. The computational demand decreases and thus the accessible time and length scale increases in this order. All three resolution levels and mixtures of them (i.e. QM/MM or hybrid atomistic/CG simulations) are relevant for computer-aided drug design. QM/MM simulations allow the study of enzymatic reactions, reaction pathways, energy barriers and proton and electron reshuffling in larger systems, where the active site is treated quantum-mechanically while the other parts of the system can be simulated at the atomic level 6, 7, 8, 9, 10. Coarse-graining of atomic models is a popular technique to access longer time scales in MD simulations. By embedding multiple atoms in a single CG bead, the number of particle–particle interactions to be calculated can be reduced considerably. However, several issues arise with this technique. First of all, coarse-graining involves a loss of information, therefore only degrees of freedom that are less important for the question of interest should be coarse-grained. Second, no general procedure to select the type and number of atoms per CG bead exists, which leads to a large number of CG models being proposed [11]. For drug design, where a high resolution of the ligand and active site is likely to be important, atomistic/CG hybrid or multiscale models are particularly interesting [12].

In classical MD simulations, the potential energy function commonly denoted as the ‘force-field’ describes the interactions between the atoms (or generally particles). It consists of so-called bonded terms between covalently bound atoms (i.e. bonds, angles, torsions) and nonbonded terms (i.e. van der Waals interactions and electrostatic interactions). The parameters of the force-field are typically fitted to reproduce data from higher-level calculations and/or experiments. For the majority of biomolecules such as proteins, DNA, lipids and sugars, a relatively small number of building blocks can be used, which have to be parameterized only once. Traditional biomolecular force-fields such as AMBER [13], OPLS [14], CHARMM [15] and GROMOS [16] have been continuously improved and tested over the years. The diversity of small organic molecules by contrast prevents the use of building blocks, and ligands have to be parameterized individually. Parameters for the bonded terms and the van der Waals interactions are typically taken from a generalized force-field, whereas the partial atomic charges are derived from a QM calculation. The force-fields most commonly used for this purpose are the general AMBER force-field (GAFF) [17], the OPLS all-atom force-field (OPLS-AA) [18], the CHARMM general force-field (CGenFF) 19, 20 and the GROMOS automatic topology builder [21].

The resulting MD trajectory of a ligand–macromolecule complex is analyzed to extract information such as distances and interactions between atoms or residues of interest. Additional values are measured to control the overall stability of the complex, such as the root-mean-square deviations (RMSD) from a reference configuration, the atom-positional root-mean-square fluctuations (RMSF) or the torsion angle distributions of the protein backbone. Additionally, an important application of MD simulations – especially in the context of drug design – is the estimation of free energy differences, ΔG. The free energy of binding, ΔG^bind, is the difference between the free energy of the ligand free in solution and bound to the protein, and is directly correlated to the binding constant, K_i (Eq. (2)):

$$\Delta G^{bind}=G^{complex}-G^{free}=RT\ln (K_i)$$

Because the direct estimation of ΔG^bind from unbiased MD simulations is difficult, the difference in ΔG^bind between two ligands A and B is often calculated. This quantity, ΔΔG_BA^bind, is accessible by MD simulations through so-called ‘alchemical’ transformations between ligand A and B. The rigorous methods for free energy calculation using MD such as thermodynamic integration (TI) [22], free energy perturbation (FEP) [23] and Bennett acceptance ratio (BAR) [24] belong to the most accurate, but also most time-consuming, techniques. Among these, TI is probably the most popular and robust method. In TI, the system is perturbed in small steps along an artificial ‘reaction coordinate’ λ from state (ligand) A to B, and the resulting curve is integrated to yield ΔG_BA (Eq. (3)):

$$\Delta G_{BA}=G({\lambda }_B)-G({\lambda }_A)={\int }_{{\lambda }_A}^{{\lambda }_B}{〈\frac{\partial H(\lambda )}{\partial \lambda }〉}_{{\lambda }^{\prime }}d{\lambda }^{\prime }$$

where H(λ) is the λ-dependent Hamiltonian of the system.

Alternative approaches that are more efficient but more approximate are the linear interaction energy (LIE) [25] and the combination of molecular mechanics (MM) calculations and continuum solvation models such as Poisson–Boltzmann surface area (PBSA) and generalized Born surface area (GBSA) 26, 27. In addition to their lower computational demands, MM/PBSA and MM/GBSA allow the direct estimation of ΔG^bind. The free energy of the ligand, the protein and the protein–ligand complex are each estimated as a sum of four terms (Eq. (4)):

$$G=E_{MM}+G_{solv}+G_{non-polar}-T{S}_{MM}$$

where G_solv is the polar solvation free energy (estimated by the PB or GB equation), G_non-polar the nonpolar solvation free energy (estimated by the solvent accessible surface area), E_MM the enthalpy and S_MM the entropy of the system (estimated by molecular mechanics). In this equation, the entropy term, S_MM, is the major limitation of this approach because it is approximated by normal mode analysis of the harmonic frequencies obtained from the MM calculation [28].

Because most MD methods assume a system in equilibrium, insufficient equilibration before the production simulation can have severe effects on the accuracy of the obtained results. This is especially important in the context of drug–receptor complexes, which tend to be large and can involve long-time motions. Similarly, sufficient sampling during the production simulation is crucial for high-quality results, leading to the development of a large variety of sampling enhancement techniques to overcome barriers and increase the sampling efficiency [29].

The different approaches discussed above can be used in computer-aided medicinal chemistry to understand drug–receptor interactions better and support the design of more-active or more-specific ligands. To gain an overall perspective on the subject, literature in the field of medicinal chemistry was surveyed for the past four years and articles using MD were collected. This analysis enables understanding of the variety of fields in which MD is used and highlights the key information this technique gives access to. In the following section, representative examples of MD are reported and discussed for studies of ligand–target complexes.