Perspectives on High-Throughput Ligand/Protein Docking With Martini MD Simulations
- 1Groningen Biomolecular Sciences and Biotechnology Institute and Zernike Institute for Advanced Materials, University of Groningen, Groningen, Netherlands
- 2PharmCADD, Busan, South Korea
- 3Molecular Microbiology and Structural Biochemistry (MMSB, UMR 5086), CNRS, University of Lyon, Lyon, France
- 4Faculty of Biomedical Sciences, Institute of Computational Science, Università della Svizzera Italiana (USI), Lugano, Switzerland
- 5Department of Pharmacy, University of Naples “Federico II”, Naples, Italy
- 6Department of Physics, Pukyong National University, Busan, South Korea
Molecular docking is central to rational drug design. Current docking techniques suffer, however, from limitations in protein flexibility and solvation models and by the use of simplified scoring functions. All-atom molecular dynamics simulations, on the other hand, feature a realistic representation of protein flexibility and solvent, but require knowledge of the binding site. Recently we showed that coarse-grained molecular dynamics simulations, based on the most recent version of the Martini force field, can be used to predict protein/ligand binding sites and pathways, without requiring any a priori information, and offer a level of accuracy approaching all-atom simulations. Given the excellent computational efficiency of Martini, this opens the way to high-throughput drug screening based on dynamic docking pipelines. In this opinion article, we sketch the roadmap to achieve this goal.
Structure-based drug design has been extensively used by pharmaceutical companies and academic research groups to reduce the cost and time necessary for the discovery of new drugs. The approach relies on the knowledge of the atomistic structure of the biological target, obtained by experiments (e.g., X-ray crystallography, NMR spectroscopy, cryo-electron microscopy) or modeling (e.g., based on homology). Standard pipelines often start with in silico docking experiments, used for virtual screening of thousands of compounds or molecular fragments (Sliwoski et al., 2014; Leelananda and Lindert, 2016; Duarte et al., 2019). After a significant reduction of the chemical space, a selected group of molecules can be optimized by all-atom (AA) molecular dynamics (MD) based simulations (Jorgensen and Thomas, 2008; Jorgensen, 2009; Vivo et al., 2016; Limongelli, 2020). AA MD simulations can be used not only to improve the prediction of the binding pose and affinity, but also to get insight into (un)binding rates and pathways (Dror et al., 2011; Shan et al., 2011; Limongelli et al., 2013; Tiwary et al., 2015; Copeland, 2016; Bruce et al., 2018). As a third step, further selection can be performed considering predictions of absorption, distribution, metabolism, excretion and toxicity (ADMET) (Van De Waterbeemd and Gifford, 2003; Cheng et al., 2013). The obtained lead compounds need to be validated by in vitro assays and structurally improved – lead optimization – to achieve drug candidates, which are tested in animal models and eventually enter clinical trials before final approval. Despite the rapid advances in computer-aided drug discovery methods, the limitations of such approaches are still major. Docking assays remain to date the first option in the drug discovery pipeline thanks to their capability of “virtually” testing thousands of molecules in a short time. However, docking accuracy is poor due to limitations in simplified energy (“scoring”) functions, sampling ligand and protein flexibility (Grinter and Zou, 2014), and representation of the environment – crucial in hydrated binding pockets and in transmembrane proteins, that represent a large fraction of the pharmaceutically relevant protein targets. AA MD simulations can tackle these limitations, but they are still computationally prohibitively expensive (Durrant and McCammon, 2011; Liu et al., 2018; Miller et al., 2020), due to relatively long time scales of conformational dynamics in proteins. Moreover, predictions of dissociation pathways and rates are extremely challenging, and require high performance computing and enhanced sampling techniques (Limongelli et al., 2013; Casasnovas et al., 2017; Brotzakis et al., 2019; Schuetz et al., 2019). Peptide and protein design for biopharmaceutical applications have similar pitfalls, with the current approaches reasonably successful in predicting protein structures (Hutson, 2019; Callaway, 2020) and rigid-body protein-protein interactions (Siebenmorgen and Zacharias, 2020) but with limitations in the design of conformational changes (Feldmeier and Höcker, 2013; Yang and Lai, 2017; Perkel, 2019; D’Annessa et al., 2020).
Coarse-grained (CG) modeling is a computationally cheaper alternative to high-resolution atomistic approaches (Ingólfsson et al., 2014; Kmiecik et al., 2016), as it reduces the computational cost by grouping atoms into effective interaction sites. Numerous CG models have been developed during the past two decades, with different levels of coarsening and different mathematical representations. CG models have been successfully applied to study a large range of processes in biology (Yen et al., 2018; Bruininks et al., 2020; Lucendo et al., 2020) and materials science (Casalini et al., 2019; Alessandri et al., 2020; Li et al., 2020; Vazquez-Salazar et al., 2020). Applications such as structure-based drug design are particularly challenging for CG modeling because of the severe requirements: (1) high chemical specificity (i.e., allowing to distinguish most chemical groups); (2) capability to represent all possible components of the system (proteins, cofactors, nucleic acids, drug candidates, waters, lipids, etc.) in a coherent way; (3) realistic representation of conformational flexibility of each molecule in the system; and (4) accurate thermodynamics and kinetics of binding. Currently, none of the CG force fields available fulfills all the requirements above, but the Martini CG force field fulfills at least some (Marrink et al., 2007; Marrink and Tieleman, 2013), as it allows modeling all main biomolecules (Monticelli et al., 2008; López et al., 2009; de Jong et al., 2013, 2015; Uusitalo et al., 2015, 2017; Wassenaar et al., 2015) with relatively high chemical specificity, and proteins may still retain reasonable conformational flexibility (Periole et al., 2009; Melo et al., 2017; Poma et al., 2017). As in AA MD simulations, most of the details of the environment can be included in Martini CG simulations, for instance an explicit solvent model or a complex bilayer composition (Ingolfsson et al., 2015; Marrink et al., 2019).
Although Martini-based CG MD simulations have been used to study a wide range of biomolecular processes, examples of protein–ligand binding are still scarce (Negami et al., 2014, 2020; Delort et al., 2017; Ferré et al., 2019; Jiang and Zhang, 2019; Dandekar and Mondal, 2020). Studies of protein-protein interactions are more common, although usually restricted to membrane environments (Baaden and Marrink, 2013; Castillo et al., 2013; Lelimousin et al., 2016; Sun et al., 2020). In some cases, binding of lipids to sites deeply buried inside the protein can be obtained by brute force Martini MD (Arnarez et al., 2013; Van Eerden et al., 2017; Corradi et al., 2019). Overall, some limiting factors hampered the use of Martini in small-molecule and protein design: (1) chemical specificity to reproduce the broad chemical space of drugs; (2) the thermodynamics of ligand-protein and protein-protein interactions are generally overestimated (Stark et al., 2013; Javanainen et al., 2017; Alessandri et al., 2019); and (3) introduction of conformational flexibility in proteins requires case-by-case optimization (Negami et al., 2020; Ahalawat and Mondal, 2021). A new version of the Martini force field, named Martini 3 (Souza et al., 2021), partly solves these issues: it can represent a broader variety of chemical compounds, and it features improved molecular packing and optimized molecular interactions (along with specific interactions mimicking H-bonding and electronic polarizability). Recently, Martini 3 was successfully applied to a range of protein-ligand system examples, from the well-characterized T4 lysozyme to members of the GPCR family and nuclear receptors to a variety of enzymes (Souza et al., 2020). In addition, combination of Martini 3 and Gō-like potentials can substantially improve the modeling of protein flexibility (Poma et al., 2017; Souza et al., 2019). Combined, these new features open the possibility of computer-aided drug design based on CG models.
In this perspective, we sketch a possible roadmap for a drug design pipeline using Martini, where no a priori information about the target pocket is necessary. Competition between ligands for different pockets and environments can be included in the screening. Protein flexibility can be incorporated to a certain degree, allowing the possible discovery of cryptic (hidden) pockets (Kuzmanic et al., 2020). Ligand (un)binding pathways are accessible via enhanced sampling techniques (Raniolo and Limongelli, 2020), and enable for the first time the possibility of a “dynamic” drug screening based not only on ligand binding modes, but also on kinetically relevant states – that is considering binding affinity and dissociation rates (i.e., drug residence time). The next sections detail the key steps of this pipeline.
Ligand Databases: Coarse-Graining the Ligands
The very first step to develop a Martini drug design pipeline is to create curated and validated databases containing hundreds to thousands of small-molecule models. This CG database needs to include molecular moieties usually found in drugs, such as halogens, heterocycles, and sulfamides. Alternatively, the databases of low-molecular-weight molecules (∼150 Da) can also be created for fragment-based drug discovery campaigns (Rognan, 2012). Parameters for molecules/fragments of pharmaceutical interest need to be validated by comparison between CG, AA and, if available, experimental data for a subset of relevant target systems. Once validated, all the models will be made available via the open-access Martini Database (MAD) web server1. The initial CG databases are also the foundation to develop and calibrate automatic tools to generate parameters for new CG models. Such automatic tools should perform AA to CG mapping [as performed by auto-martini (Bereau and Kremer, 2015)], bead assignment (i.e., the choice of the CG interaction parameters), and determination of the bonded parameters [as PyCGTOOL (Graham et al., 2017) or Swarm-CG (Empereur-mot et al., 2020)], allowing further coverage of chemical space. The creation of accurate databases and integration of automatic tools is currently one of the main bottlenecks hampering high-throughput screening with Martini.
Virtual Screening: Martini Dynamic Docking
Virtual screening is the core of the drug design pipeline, and usually relies on docking algorithms. The use of Martini CG models will enable a new approach: dynamic docking with no a priori knowledge of the binding pocket in the target structure. The concept here is to sample protein–ligand interactions with CG MD simulations, which is around 300 to 1,000 times faster than atomistic MD (Souza et al., 2020). A practical example of such speed up can be given for propranolol binding to β2 adrenergic receptor, which has been simulated in atomistic (Dror et al., 2011) and coarse-grained (Souza et al., 2020) resolution. Atomistic simulations showed one binding event every 11.9 μs, which for a single simulation would take 84 days of computing time (using the 4 CPUs and 1 GPU in a computer/conditions described in the performance tests of Souza et al., 2020). The same system in CG simulations showed roughly the same number of binding events per μs (considering a normalization based in the different concentration of ligands), and would take 2 to 7 h of computing time, on the same hardware. We remark that a fair comparison between coarse-grained and atomistic simulation time is not trivial, since this should consider the different simulation conditions (e.g., ligand concentration) and parameters (Souza et al., 2020).
Multiple strategies are possible to accelerate sampling even more, with different computational costs and different levels of sophistication. Unbiased MD simulations could be applied in certain cases, to obtain not only binding poses but also estimates of binding affinities, as recently demonstrated for T4 lysozyme (Souza et al., 2020). However, for a general approach to virtual screening, faster methods are necessary. One possibility is to combine CG models with enhanced sampling techniques that do not depend on prior knowledge of the binding pathways. Examples are Gaussian accelerated molecular dynamics (GaMD) (Miao et al., 2015; Pang et al., 2017), and Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) (Wang et al., 2013; Luitz and Zacharias, 2014). Computational performance can be straightforwardly increased by optimizing ligand concentration, to increase the probability of binding. The approach was already tested with atomistic simulations in a variety of systems (Dror et al., 2011; Shan et al., 2011; Decherchi et al., 2015; Schneider et al., 2016; Mondal et al., 2018). To avoid ligand aggregation, artificial repulsive interactions among ligands may be used (Shan et al., 2011). Similar strategies are also extensively used in so-called mixed-solvent (or co-solvent) approaches, where high concentrations of fragments are used to identify and stabilize cryptic pockets (Guvench and MacKerell, 2009; Bakan et al., 2012; Schmidt et al., 2019; Kuzmanic et al., 2020). Another idea is to only use isolated beads as probes representing chemical groups or fragments, to predict the chemical topology in pockets and generate pharmacophore models (Michelarakis et al., 2018; Michelarakis, 2019). The combination of CG models, enhanced sampling, and ligand/fragment concentration strategies will allow simulations of competitive binding assays.
An advantage of Martini dynamic docking approach is the improved representation of protein flexibility via Gō-like potentials (Poma et al., 2017; Souza et al., 2019). Although some docking strategies can also include protein flexibility (Amaro et al., 2018; Evangelista Falcon et al., 2019), they usually depend on prior sampling of the protein conformational space, followed by docking in a specific chosen pocket. In the strategy proposed here, no a priori selection of the binding pocket is needed. Both induced-fit and conformational-selection mechanisms are included in MD simulations, as recently demonstrated (Souza et al., 2020); however, accuracy will depend on the quality of the protein CG model.
Another major advantage is the possibility to include complex environments, such as multicomponent membranes, crowded protein solutions, or other relevant in vivo-like conditions, allowing more realistic predictions. Competition with the environment may be relevant for proper interpretation of ligand biological activity. For instance, lipid membrane composition may affect kinetic rates and (un)binding constants in GPCRs (Vauquelin, 2010; Sykes et al., 2014, 2019; Yuan et al., 2018) by altering ligand partitioning to the membrane where the target protein is located. Atomistic MD simulations of such complex systems are computationally very costly, while they are already within reach with Martini (Marrink et al., 2019).
Combining “standard” docking algorithms with Martini provides a computationally cheap alternative to all-atom docking. As recently demonstrated by HADDOCK (Honorato et al., 2019; Roel-Touris et al., 2019), docking with Martini can be one order of magnitude faster than atomistic docking. This would allow to routinely explore very large ligand datasets (Lyu et al., 2019) or even to use massive docking with grids covering the whole the protein, or exploring multiple proteins/conformations at the same time. However, common problems of docking approaches (mentioned above) still would be present; probably Martini MD approaches represent a better compromise between accuracy and computational performance.
Lead Optimization: Backmapping and Coarse Graining in Chemical Space
Accurate predictions of ligand binding poses and affinities are key aspects for lead optimization (Jorgensen, 2009; Vivo et al., 2016). In atomistic pipelines, MD simulations can be used as a post-processing tool to validate and/or refine the binding poses from docking (Vivo et al., 2016). After this first check, more rigorous estimates of ligand binding affinities can be achieved by free energy perturbation (FEP) or thermodynamics integration (TI) (Jorgensen and Thomas, 2008; Jorgensen, 2009) – methods based on conversion of one ligand to another, allowing to add or replace substituents, in order to optimize ligand-protein interactions. In a Martini drug design pipeline (step 3A of Figure 1), one could simply convert the CG representation to all-atom (“backmapping” procedure) to verify and refine the CG docking poses. Currently, the most reliable approach for backmapping is the geometric projection implemented in Backward (Wassenaar et al., 2014). The main disadvantage is the need for mapping files for each ligand. After obtaining the atomistic structures, any MD-based simulations can be straightforwardly used. Careful equilibration is necessary to allow relaxation of the system, in particular, the water molecules may need to fill small cavities in pockets not accessible to CG water. One possibility is to model buried water molecules or ions using smaller beads, as previously showcased (Souza et al., 2020). Such difficulties are also common in standard docking approaches, as they usually do not include water molecules.
Figure 1. One of the possible pipelines for high-throughput dynamic docking based on Martini coarse-grained modeling. (1) The first step in the pipeline is the automatic conversion of input libraries of small compounds to Martini models. The library includes drug-like compounds and small-sized rigid molecules, useful for fragment-based drug discovery. (2) In the second step, thousands of parallel simulations are automatically set up, to sample small-molecule binding to pockets in the target protein. Competition in silico assays with endogenous ligands are possible in this step. Performance can be straightforwardly increased by optimizing the ligand concentration as well as by employing enhanced sampling techniques. At the end, automatic analysis and ranking of ligands is performed, to obtain estimates of binding affinity in relation to different pockets and environments (e.g., binding to protein in relation to water and/or bilayer). (3A) After defining the pocket and a set of candidates, the accuracy of the prediction can be improved in third step: backmapping to the atomistic models can be performed, providing high-resolution details of the binding modes. Additionally, free energy perturbation (FEP) or thermodynamic integration (TI) estimating the energetic cost of converting certain chemical groups into others can allow further optimization of the molecular structure. Here, coarse-graining in the chemical space is possible, as Martini CG moieties can represent more than one chemical fragment at the same time. (3B) An alternative or complementary third step, based on binding affinity and kinetics, is also considered here. Analysis of trajectories obtained in step 2 can help to identify the drug (un)binding pathways, which can be used in methods as Funnel-Metadynamics to provide lowest energy binding modes and dissociation rates koff (drug residence time) states determining. (4) The combined analysis of steps II and III can be used for predictions of activity, which in combination with ADMET predictions leads to the final rankings and selection of the lead compounds for in vitro assays. Part of the figure is adapted from Souza et al. (2020).
An alternative possibility for lead optimization in Martini would be to reverse the order of the steps, performing first a preliminary set of FEP/TI calculations at the Martini CG level. Such approach would allow to explore a broader portion of the chemical space. On top of the default computational efficiency of CG models, additional speed up could be obtained. First, given the smoother potential surface, the replacement of one bead for another (representing different chemical groups) could be performed in less FEP/TI windows. Additionally, as Martini CG beads generally represent more than one chemical fragment (Menichetti et al., 2019; Bereau, 2020), the exploration of chemical space increases computational efficiency by an additional factor 103–104 (Menichetti et al., 2019; Bereau, 2020) thanks to the reduction in the size of chemical space. Each bead of the CG model can be transformed into different chemical groups, for instance by using different mapping files for each bead in the Backward code (Wassenaar et al., 2014). With this alternative lead optimization approach, backmapping would be performed as the last step, to increase accuracy of the predictions.
Alternative Route: Ligand Binding Pathways, Binding Affinities, and Kinetic Rates
Drug discovery is historically focused on the elucidation and optimization of the ligand binding mode and binding affinity. However, in vivo drug activity is quantitatively correlated to the drug residence time – i.e., dissociation constant rate koff – more than binding affinity Kb (Copeland et al., 2006). The idea of integrating kinetic data in drug screening has been around since the beginning of 2000s (Limongelli, 2020; Nunes-Alves et al., 2020). However, ligand binding kinetics is determined by rare events, crossing ephemeral, high-energy states, elusive to both experiments, and computations (Copeland, 2016). The recent proof of concept with Martini 3 (Souza et al., 2020) opens the possibility of including information on ligand binding pathways in drug design pipelines (step 3B of Figure 1). The data coming from unbiased CG MD simulations should be integrated in a rigorous theoretical framework. One possibility would be to use Markov state models (Husic and Pande, 2018) based on Martini dynamic docking screening (step 2 of in Figure 1). The method has proven useful in atomistic ligand binding simulation studies (Buch et al., 2011) but it shows difficulties in defining the macrostates of the process, the choice of lag-time, and the sampling necessary to ensure statistical significance. An attractive strategy is to combine CG MD with Funnel-Metadynamics (Limongelli et al., 2013; Raniolo and Limongelli, 2020) that has emerged as a powerful method to reproduce binding mechanisms in ligand/protein and ligand/DNA complexes, identify crystallographic binding modes and predict binding free energies (Troussicot et al., 2015; Comitani et al., 2016; Moraca et al., 2017; Saleh et al., 2017; Yuan et al., 2018; D’Annessa et al., 2019). During FM simulations, the whole drug binding mechanism is reproduced, from the fully solvated state to the final binding mode, allowing to disclose important aspects of the binding process such as (i) the presence of alternative binding modes; (ii) the role of the solvent; and (iii) the kinetically relevant states (Tiwary et al., 2015; Brotzakis et al., 2019; Raniolo and Limongelli, 2020). CG-FM allows quantitative predictions of koff and Kb, ligand binding modes, and rate determining steps (Figure 1). This advance will represent a paradigm shift in drug design, as medicinal chemists would optimize the structure of drug candidates not only based on the static representation of the ligand binding mode, but also on the structures of kinetically relevant states. We point out that the reduction of friction from the missing atomistic degrees of freedom speeds up CG dynamics and affects kinetic estimates. However, estimating trends may be useful enough for ligand screening, while realistic kinetics rates might be recovered from estimates of the friction reduction (Español and Zúñiga, 2011).
Further Considerations and Discussion
We described a new vision of high-throughput drug screening based on Martini CG models. Although most of the recent efforts in new drug design approaches focused on artificial intelligence (AI), the development of new methods covering gaps in standard approaches is equally important. Machine learning and other AI approaches have great advantages when tackling problems with enough experimental data to be used as training dataset. In situations where this is not the case, physics-based approaches (such as CG molecular dynamics) can perform better. In particular, structural databases of transmembrane proteins are still limited. The same is also true for databases that include dynamic information, which can be important to elucidate hidden allosteric pockets, to properly model fit-induced ligand binding process or to determine ligand association/dissociation pathways. More than complementary, AI and physics-based approaches can be combined, with CG MD simulations being used for the training of AI models or for the further refinement of AI predictions.
The proposed Martini drug design workflow (Figure 1) could be applied in full, or specific modules could be adapted in more traditional virtual screening campaigns. Screening of drugs based on ligand binding pathways and dissociation rates is currently out of reach for all-atom descriptions, due to the prohibitively high computational cost. Flexible proteins in complex environments are also too costly for all-atom docking approaches. Martini greatly reduces the computational costs of MD, while offering reasonable accuracy and structural detail. Accuracy will be further improved with the implementation of polarizable models (Yesylevskyy et al., 2010; de Jong et al., 2013; Michalowsky et al., 2017, 2018; Khan et al., 2020). Additionally, protonation state changes and pH effects can be included with Titratable Martini approaches (Grünewald et al., 2020). Also within reach is the design of epitopes and nucleic acids, useful for rational vaccine development (Kulp and Schief, 2013; Hodgson, 2020; Norman et al., 2020). In this context, even CG MD simulations may be overly expensive, as the approach demands scanning of protein-protein and protein-nucleic acid interfaces. Here, combination with standard docking is already a reality, as recently implemented in HADDOCK (Honorato et al., 2019; Roel-Touris et al., 2019). Overall, we believe dynamic docking with CG models has great innovation potential, both in academic and private sectors, and we hope this Perspective will contribute to motivate the modeling community to expand the efforts in this area.
Data Availability Statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.
PCTS wrote the first draft of the manuscript and prepared the figure. All authors contributed to the conception of the perspective article, manuscript revision, read, and approved the submitted version.
VL acknowledges the support from the European Research Council (ERC Consolidator Grant “CoMMBi”), the Swiss National Science Foundation (Project No. 200021_163281), the Italian MIUR-PRIN 2017 (2017FJZZRC), and the Swiss National Supercomputing Centre (CSCS). SM acknowledges funding from the ERC through an Advanced grant “COMP-MICR-CROW-MEM.” LM acknowledges the Institut National de la Santé et de la Recherche Medicale (INSERM) and the Agence Nationale de la Recherche (ANR) for funding (Grant Nos. ANR-17-CE11-0003 and ANR-20-CE13-0030-03).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Ahalawat, N., and Mondal, J. (2021). An appraisal of computer simulation approaches in elucidating biomolecular recognition pathways. J. Phys. Chem. Lett. 12, 633–641. doi: 10.1021/acs.jpclett.0c02785
Alessandri, R., Sami, S., Barnoud, J., Vries, A. H., Marrink, S. J., and Havenith, R. W. A. (2020). Resolving donor–acceptor interfaces and charge carrier energy levels of organic semiconductors with polar side chains. Adv. Funct. Mater. 30:2004799. doi: 10.1002/adfm.202004799
Alessandri, R., Souza, P. C. T., Thallmair, S., Melo, M. N., de Vries, A. H., and Marrink, S. J. (2019). Pitfalls of the Martini model. J. Chem. Theory Comput. 15, 5448–5460. doi: 10.1021/acs.jctc.9b00473
Arnarez, C., Mazat, J.-P., Elezgaray, J., Marrink, S.-J., and Periole, X. (2013). Evidence for cardiolipin binding sites on the membrane-exposed surface of the cytochrome bc1. J. Am. Chem. Soc. 135, 3112–3120. doi: 10.1021/ja310577u
Bakan, A., Nevins, N., Lakdawala, A. S., and Bahar, I. (2012). Druggability assessment of allosteric proteins by dynamics simulations in the presence of probe molecules. J. Chem. Theory Comput. 8, 2435–2447. doi: 10.1021/ct300117j
Bereau, T. (2020). Computational Compound Screening of Biomolecules and Soft Materials by Molecular Simulations. Available online at: http://arxiv.org/abs/2010.03298 (accessed November 20, 2020).
Bereau, T., and Kremer, K. (2015). Automated parametrization of the coarse-grained Martini force field for small organic molecules. J. Chem. Theory Comput. 11, 2783–2791. doi: 10.1021/acs.jctc.5b00056
Brotzakis, Z. F., Faidon Brotzakis, Z., Limongelli, V., and Parrinello, M. (2019). Accelerating the calculation of protein–ligand binding free energy and residence times using dynamically optimized collective variables. J. Chem. Theory Comput. 15, 743–750. doi: 10.1021/acs.jctc.8b00934
Bruce, N. J., Ganotra, G. K., Kokh, D. B., Kashif Sadiq, S., and Wade, R. C. (2018). New approaches for computing ligand–receptor binding kinetics. Curr. Opin. Struct. Biol. 49, 1–10. doi: 10.1016/j.sbi.2017.10.001
Buch, I., Giorgino, T., and De Fabritiis, G. (2011). Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 108, 10184–10189. doi: 10.1073/pnas.1103547108
Casalini, T., Limongelli, V., Schmutz, M., Som, C., Jordan, O., Wick, P., et al. (2019). Molecular modeling for nanomaterial–biology interactions: opportunities, challenges, and perspectives. Front. Bioeng. Biotechnol. 7:268. doi: 10.3389/fbioe.2019.00268
Casasnovas, R., Limongelli, V., Tiwary, P., Carloni, P., and Parrinello, M. (2017). Unbinding kinetics of a p38 MAP kinase Type II inhibitor from metadynamics simulations. J. Am. Chem. Soc. 139, 4780–4788. doi: 10.1021/jacs.6b12950
Castillo, N., Monticelli, L., Barnoud, J., and Tieleman, D. P. (2013). Free energy of WALP23 dimer association in DMPC, DPPC, and DOPC bilayers. Chem. Phys. Lipids 169, 95–105. doi: 10.1016/j.chemphyslip.2013.02.001
Cheng, F., Li, W., Liu, G., and Tang, Y. (2013). In silico ADMET prediction: recent advances, current challenges and future trends. Curr. Top. Med. Chem. 13, 1273–1289. doi: 10.2174/15680266113139990033
Comitani, F., Limongelli, V., and Molteni, C. (2016). The free energy landscape of GABA binding to a pentameric ligand-gated ion channel and its disruption by mutations. J. Chem. Theory Comput. 12, 3398–3406. doi: 10.1021/acs.jctc.6b00303
Corradi, V., Sejdiu, B. I., Mesa-Galloso, H., Abdizadeh, H., Noskov, S. Y., Marrink, S. J., et al. (2019). Emerging diversity in lipid–protein interactions. Chem. Rev. 119, 5775–5848. doi: 10.1021/acs.chemrev.8b00451
D’Annessa, I., Di Leva, F. S., La Teana, A., Novellino, E., Limongelli, V., and Di Marino, D. (2020). Bioinformatics and biosimulations as toolbox for peptides and peptidomimetics design: where are we? Front. Mol. Biosci. 7:66. doi: 10.3389/fmolb.2020.00066
D’Annessa, I., Raniolo, S., Limongelli, V., Di Marino, D., and Colombo, G. (2019). Ligand binding, unbinding, and allosteric effects: deciphering small-molecule modulation of HSP90. J. Chem. Theory Comput. 15, 6368–6381. doi: 10.1021/acs.jctc.9b00319
de Jong, D. H., Liguori, N., van den Berg, T., Arnarez, C., Periole, X., and Marrink, S. J. (2015). Atomistic and coarse grain topologies for the cofactors associated with the photosystem II core complex. J. Phys. Chem. B 119, 7791–7803. doi: 10.1021/acs.jpcb.5b00809
de Jong, D. H., Singh, G., Bennett, W. F. D., Arnarez, C., Wassenaar, T. A., Schäfer, L. V., et al. (2013). Improved parameters for the Martini coarse-grained protein force field. J. Chem. Theory Comput. 9, 687–697. doi: 10.1021/ct300646g
Decherchi, S., Berteotti, A., Bottegoni, G., Rocchia, W., and Cavalli, A. (2015). The ligand binding mechanism to purine nucleoside phosphorylase elucidated via molecular dynamics and machine learning. Nat. Commun. 6:6155. doi: 10.1038/ncomms7155
Delort, B., Renault, P., Charlier, L., Raussin, F., Martinez, J., and Floquet, N. (2017). Coarse-grained prediction of peptide binding to G-protein coupled receptors. J. Chem. Inf. Model. 57, 562–571. doi: 10.1021/acs.jcim.6b00503
Dror, R. O., Pan, A. C., Arlow, D. H., Borhani, D. W., Maragakis, P., Shan, Y., et al. (2011). Pathway and mechanism of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. U. S. A. 108, 13118–13123. doi: 10.1073/pnas.1104614108
Duarte, Y., Márquez-Miranda, V., Miossec, M. J., and González-Nilo, F. (2019). Integration of target discovery, drug discovery and drug delivery: a review on computational strategies. Wiley Interdiscip. Rev. Nanomed. Nanobiotechnol. 11:e1554. doi: 10.1002/wnan.1554
Empereur-mot, C., Pesce, L., Bochicchio, D., Perego, C., and Pavan, G. M. (2020). Swarm-CG: automatic parametrization of bonded terms in coarse-grained models of simple to complex molecules via fuzzy self-tuning particle swarm optimization. ACS Omega 5, 32823–32843. doi: 10.26434/chemrxiv.12613427.v2
Evangelista Falcon, W., Ellingson, S. R., Smith, J. C., and Baudry, J. (2019). Ensemble docking in drug discovery: how many protein configurations from molecular dynamics simulations are needed to reproduce known ligand binding? J. Phys. Chem. B 123, 5189–5195. doi: 10.1021/acs.jpcb.8b11491
Ferré, G., Louet, M., Saurel, O., Delort, B., Czaplicki, G., M’Kadmi, C., et al. (2019). Structure and dynamics of G protein-coupled receptor–bound ghrelin reveal the critical role of the octanoyl chain. Proc. Natl. Acad. Sci. U. S. A. 116, 17525–17530. doi: 10.1073/pnas.1905105116
Graham, J. A., Essex, J. W., and Khalid, S. (2017). PyCGTOOL: automated generation of coarse-grained molecular dynamics models from atomistic trajectories. J. Chem. Inf. Model. 57, 650–656. doi: 10.1021/acs.jcim.7b00096
Grünewald, F., Souza, P. C. T., Abdizadeh, H., Barnoud, J., de Vries, A. H., and Marrink, S. J. (2020). Titratable Martini model for constant pH simulations. J. Chem. Phys. 153:024118. doi: 10.1063/5.0014258
Guvench, O., and MacKerell, A. D. Jr. (2009). Computational fragment-based binding site identification by ligand competitive saturation. PLoS Comput. Biol. 5:e1000435. doi: 10.1371/journal.pcbi.1000435
Ingólfsson, H. I., Lopez, C. A., Uusitalo, J. J., de Jong, D. H., Gopal, S. M., Periole, X., et al. (2014). The power of coarse graining in biomolecular simulations. Wiley Interdiscip. Rev. Comput. Mol. Sci. 4, 225–248. doi: 10.1002/wcms.1169
Khan, H. M., Souza, P. C. T., Thallmair, S., Barnoud, J., de Vries, A. H., Marrink, S. J., et al. (2020). Capturing choline-aromatics cation-π interactions in the Martini force field. J. Chem. Theory Comput. 16, 2550–2560. doi: 10.1021/acs.jctc.9b01194
Kmiecik, S., Gront, D., Kolinski, M., Wieteska, L., Dawid, A. E., and Kolinski, A. (2016). Coarse-grained protein models and their applications. Chem. Rev. 116, 7898–7936. doi: 10.1021/acs.chemrev.6b00163
Kuzmanic, A., Bowman, G. R., Juarez-Jimenez, J., Michel, J., and Gervasio, F. L. (2020). Investigating cryptic binding sites by molecular dynamics simulations. Acc. Chem. Res. 53, 654–661. doi: 10.1021/acs.accounts.9b00613
Lelimousin, M., Limongelli, V., and Sansom, M. S. P. (2016). Conformational changes in the epidermal growth factor receptor: role of the transmembrane domain investigated by coarse-grained metadynamics free energy calculations. J. Am. Chem. Soc. 138, 10611–10622. doi: 10.1021/jacs.6b05602
Li, C., Iscen, A., Sai, H., Sato, K., Sather, N. A., Chin, S. M., et al. (2020). Supramolecular–covalent hybrid polymers for light-activated mechanical actuation. Nat. Mater. 19, 900–909. doi: 10.1038/s41563-020-0707-7
López, C. A., Rzepiela, A. J., de Vries, A. H., Dijkhuizen, L., Hünenberger, P. H., and Marrink, S. J. (2009). Martini coarse-grained force field: extension to carbohydrates. J. Chem. Theory Comput. 5, 3195–3210. doi: 10.1021/ct900313w
Lucendo, E., Sancho, M., Lolicato, F., Javanainen, M., Kulig, W., Leiva, D., et al. (2020). Mcl-1 and Bok transmembrane domains: unexpected players in the modulation of apoptosis. Proc. Natl. Acad. Sci. U. S. A. 117, 27980–27988. doi: 10.1073/pnas.2008885117
Marrink, S. J., Corradi, V., Souza, P. C. T., Ingólfsson, H. I., Peter Tieleman, D., and Sansom, M. S. P. (2019). Computational modeling of realistic cell membranes. Chem. Rev. 119, 6184–6226. doi: 10.1021/acs.chemrev.8b00460
Marrink, S. J., Jelger Risselada, H., Yefimov, S., Peter Tieleman, D., and de Vries, A. H. (2007). The MARTINI force field: coarse grained model for biomolecular simulations. J. Phys. Chem. B 111, 7812–7824. doi: 10.1021/jp071097f
Melo, M. N., Arnarez, C., Sikkema, H., Kumar, N., Walko, M., Berendsen, H. J. C., et al. (2017). High-throughput simulations reveal membrane-mediated effects of alcohols on MscL gating. J. Am. Chem. Soc. 139, 2664–2671. doi: 10.1021/jacs.6b11091
Miao, Y., Feher, V. A., and McCammon, J. A. (2015). Gaussian accelerated molecular dynamics: unconstrained enhanced sampling and free energy calculation. J. Chem. Theory Comput. 11, 3584–3595. doi: 10.1021/acs.jctc.5b00436
Michalowsky, J., Schäfer, L. V., Holm, C., and Smiatek, J. (2017). A refined polarizable water model for the coarse-grained MARTINI force field with long-range electrostatic interactions. J. Chem. Phys. 146:054501. doi: 10.1063/1.4974833
Michelarakis, N., Sands, Z. A., Sansom, M. S. P., and Stansfeld, P. J. (2018). Towards dynamic pharmacophore models by coarse grained molecular dynamics. Biophys. J. 114:558a. doi: 10.1016/j.bpj.2017.11.3050
Miller, E., Murphy, R., Sindhikara, D., Borrelli, K., Grisewood, M., Ranalli, F., et al. (2020). A reliable and accurate solution to the induced fit docking problem for protein-ligand binding. ChemRxiv [Preprint]. doi: 10.26434/chemrxiv.11983845.v2
Mondal, J., Ahalawat, N., Pandit, S., Kay, L. E., and Vallurupalli, P. (2018). Atomic resolution mechanism of ligand binding to a solvent inaccessible cavity in T4 lysozyme. PLoS Comput. Biol. 14:e1006180. doi: 10.1371/journal.pcbi.1006180
Monticelli, L., Kandasamy, S. K., Periole, X., Larson, R. G., Peter Tieleman, D., and Marrink, S.-J. (2008). The MARTINI coarse-grained force field: extension to proteins. J. Chem. Theory Comput. 4, 819–834. doi: 10.1021/ct700324x
Moraca, F., Amato, J., Ortuso, F., Artese, A., Pagano, B., Novellino, E., et al. (2017). Ligand binding to telomeric G-quadruplex DNA investigated by funnel-metadynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 114, E2136–E2145. doi: 10.1073/pnas.1612627114
Negami, T., Shimizu, K., and Terada, T. (2020). Coarse-grained molecular dynamics simulation of protein conformational change coupled to ligand binding. Chem. Phys. Lett. 742:137144. doi: 10.1016/j.cplett.2020.137144
Norman, R. A., Ambrosetti, F., Bonvin, A. M. J. J., Colwell, L. J., Kelm, S., Kumar, S., et al. (2020). Computational approaches to therapeutic antibody design: established methods and emerging trends. Brief. Bioinform. 21, 1549–1567. doi: 10.1093/bib/bbz095
Periole, X., Cavalli, M., Marrink, S.-J., and Ceruso, M. A. (2009). Combining an elastic network with a coarse-grained molecular force field: structure, dynamics, and intermolecular recognition. J. Chem. Theory Comput. 5, 2531–2543. doi: 10.1021/ct9002114
Poma, A. B., Cieplak, M., and Theodorakis, P. E. (2017). Combining the MARTINI and structure-based coarse-grained approaches for the molecular dynamics studies of conformational transitions in proteins. J. Chem. Theory Comput. 13, 1366–1374. doi: 10.1021/acs.jctc.6b00986
Roel-Touris, J., Don, C. G., Honorato, R. V., João, P. G. L., and Bonyin, A. M. J. (2019). Less is more: coarse-grained integrative modeling of large biomolecular assemblies with HADDOCK. J. Chem. Theory Comput. 15, 6358–6367. doi: 10.1021/acs.jctc.9b00310
Saleh, N., Ibrahim, P., Saladino, G., Gervasio, F. L., and Clark, T. (2017). An efficient metadynamics-based protocol to model the binding affinity and the transition state ensemble of G-protein-coupled receptor ligands. J. Chem. Inf. Model. 57, 1210–1217. doi: 10.1021/acs.jcim.6b00772
Schmidt, D., Boehm, M., McClendon, C. L., Torella, R., and Gohlke, H. (2019). Cosolvent-enhanced sampling and unbiased identification of cryptic pockets suitable for structure-based drug design. J. Chem. Theory Comput. 15, 3331–3343. doi: 10.1021/acs.jctc.8b01295
Schneider, S., Provasi, D., and Filizola, M. (2016). How oliceridine (TRV-130) binds and stabilizes a μ-opioid receptor conformational state that selectively triggers G protein signaling pathways. Biochemistry 55, 6456–6466. doi: 10.1021/acs.biochem.6b00948
Schuetz, D. A., Bernetti, M., Bertazzo, M., Musil, D., Eggenweiler, H.-M., Recanatini, M., et al. (2019). Predicting residence time and drug unbinding pathway through scaled molecular dynamics. J. Chem. Inf. Model. 59, 535–549. doi: 10.1021/acs.jcim.8b00614
Shan, Y., Kim, E. T., Eastwood, M. P., Dror, R. O., Seeliger, M. A., and Shaw, D. E. (2011). How does a drug molecule find its target binding site? J. Am. Chem. Soc. 133, 9181–9183. doi: 10.1021/ja202726y
Souza, P. C. T., Alessandri, R., Barnoud, J., Thallmair, S., Faustino, I., Grunewald, F., et al. (2021). Martini 3: a general purpose force field for coarse-grain molecular dynamics. Nat. Methods. doi: 10.1038/s41592-021-01098-3
Souza, P. C. T., Thallmair, S., Conflitti, P., Ramírez-Palacios, C., Alessandri, R., Raniolo, S., et al. (2020). Protein-ligand binding with the coarse-grained Martini model. Nat. Commun. 11:3714. doi: 10.1038/s41467-020-17437-5
Souza, P. C. T., Thallmair, S., Marrink, S. J., and Mera-Adasme, R. (2019). An allosteric pathway in copper, zinc superoxide dismutase unravels the molecular mechanism of the G93A amyotrophic lateral sclerosis-linked mutation. J. Phys. Chem. Lett. 10, 7740–7744. doi: 10.1021/acs.jpclett.9b02868
Stark, A. C., Andrews, C. T., and Elcock, A. H. (2013). Toward optimized potential functions for protein–protein interactions in aqueous solutions: osmotic second virial coefficient calculations using the MARTINI coarse-grained force field. J. Chem. Theory Comput. 9, 4176–4185. doi: 10.1021/ct400008p
Sun, F., Schroer, C. F. E., Palacios, C. R., Xu, L., Luo, S.-Z., and Marrink, S. J. (2020). Molecular mechanism for bidirectional regulation of CD44 for lipid raft affiliation by palmitoylations and PIP2. PLoS Comput. Biol. 16:e1007777. doi: 10.1371/journal.pcbi.1007777
Sykes, D. A., Parry, C., Reilly, J., Wright, P., Fairhurst, R. A., and Charlton, S. J. (2014). Observed drug-receptor association rates are governed by membrane affinity: the importance of establishing “micro-pharmacokinetic/pharmacodynamic relationships” at the β2-adrenoceptor. Mol. Pharmacol. 85, 608–617. doi: 10.1124/mol.113.090209
Tiwary, P., Limongelli, V., Salvalaglio, M., and Parrinello, M. (2015). Kinetics of protein–ligand unbinding: predicting pathways, rates, and rate-limiting steps. Proc. Natl. Acad. Sci. U.S.A. 112, E386–E391. doi: 10.1073/pnas.1424461112
Troussicot, L., Guillière, F., Limongelli, V., Walker, O., and Lancelin, J.-M. (2015). Funnel-metadynamics and solution NMR to estimate protein–ligand affinities. J. Am. Chem. Soc. 137, 1273–1281. doi: 10.1021/ja511336z
Uusitalo, J. J., Ingólfsson, H. I., Akhshi, P., Peter Tieleman, D., and Marrink, S. J. (2015). Martini coarse-grained force field: extension to DNA. J. Chem. Theory Comput. 11, 3932–3945. doi: 10.1021/acs.jctc.5b00286
Van Eerden, F. J., Melo, M. N., Frederix, P. W. J. M., Periole, X., and Marrink, S. J. (2017). Exchange pathways of plastoquinone and plastoquinol in the photosystem II complex. Nat. Commun. 8:15214. doi: 10.1038/ncomms15214
Vazquez-Salazar, L. I., Selle, M., de Vries, A. H., Marrink, S. J., and Souza, P. C. T. (2020). Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid–liquid extraction. Green Chem. 22, 7376–7386. doi: 10.1039/d0gc01823f
Vivo, M. D., De Vivo, M., Masetti, M., Bottegoni, G., and Cavalli, A. (2016). Role of molecular dynamics and related methods in drug discovery. J. Med. Chem. 59, 4035–4061. doi: 10.1021/acs.jmedchem.5b01684
Wang, K., Chodera, J. D., Yang, Y., and Shirts, M. R. (2013). Identifying ligand binding sites and poses using GPU-accelerated Hamiltonian replica exchange molecular dynamics. J. Comput. Aided Mol. Des. 27, 989–1007. doi: 10.1007/s10822-013-9689-8
Wassenaar, T. A., Ingólfsson, H. I., Böckmann, R. A., Peter Tieleman, D., and Marrink, S. J. (2015). Computational lipidomics with insane: a versatile tool for generating custom membranes for molecular simulations. J. Chem. Theory Comput. 11, 2144–2155. doi: 10.1021/acs.jctc.5b00209
Wassenaar, T. A., Pluhackova, K., Böckmann, R. A., Marrink, S. J., and Tieleman, D. P. (2014). Going backward: a flexible geometric approach to reverse transformation from coarse grained to atomistic models. J. Chem. Theory Comput. 10, 676–690. doi: 10.1021/ct400617g
Yen, H.-Y., Hoi, K. K., Liko, I., Hedger, G., Horrell, M. R., Song, W., et al. (2018). PtdIns(4,5)P2 stabilizes active states of GPCRs and enhances selectivity of G-protein coupling. Nature 559, 423–427. doi: 10.1038/s41586-018-0325-6
Yesylevskyy, S. O., Schäfer, L. V., Sengupta, D., and Marrink, S. J. (2010). Polarizable water model for the coarse-grained Martini force field. PLoS Comput. Biol. 6:e1000810. doi: 10.1371/journal.pcbi.1000810
Keywords: molecular dynamics, coarse-grain, ligand-protein, protein-protein interaction, Martini, dynamic docking, high-throughput screening, drug design
Citation: Souza PCT, Limongelli V, Wu S, Marrink SJ and Monticelli L (2021) Perspectives on High-Throughput Ligand/Protein Docking With Martini MD Simulations. Front. Mol. Biosci. 8:657222. doi: 10.3389/fmolb.2021.657222
Received: 22 January 2021; Accepted: 05 March 2021;
Published: 29 March 2021.
Edited by:Agnel Praveen Joseph, Science and Technology Facilities Council, United Kingdom
Reviewed by:Valeria Losasso, United Kingdom Research and Innovation, United Kingdom
Sophie Sacquin-Mora, UPR 9080 Laboratoire de Biochimie Théorique (LBT), France
Copyright © 2021 Souza, Limongelli, Wu, Marrink and Monticelli. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Paulo C. T. Souza, firstname.lastname@example.org; Sangwook Wu, email@example.com; Vittorio Limongelli, firstname.lastname@example.org; Siewert J. Marrink, email@example.com; Luca Monticelli, firstname.lastname@example.org