Protein–ligand binding with the coarse-grained Martini model

The detailed understanding of the binding of small molecules to proteins is the key for the development of novel drugs or to increase the acceptance of substrates by enzymes. Nowadays, computer-aided design of protein–ligand binding is an important tool to accomplish this task. Current approaches typically rely on high-throughput docking essays or computationally expensive atomistic molecular dynamics simulations. Here, we present an approach to use the recently re-parametrized coarse-grained Martini model to perform unbiased millisecond sampling of protein–ligand interactions of small drug-like molecules. Remarkably, we achieve high accuracy without the need of any a priori knowledge of binding pockets or pathways. Our approach is applied to a range of systems from the well-characterized T4 lysozyme over members of the GPCR family and nuclear receptors to a variety of enzymes. The presented results open the way to high-throughput screening of ligand libraries or protein mutations using the coarse-grained Martini model.

P rotein-ligand binding interaction is fundamental to a large variety of cellular functions, including enzymatic reactions, catalysis, signal transduction, and regulation. The importance of protein-ligand interactions explains the ongoing interest to redesign these interactions, conferring novel functions or finding suitable drugs and molecular targets 1 . A typical binding pocket consisting of 10 residues gives rise to 20 10 possible mutations; together with an estimated amount of 10 60 potential drug-like compounds, this composes a vast chemical space. Therefore, the potential for rational design is enormous. To harness this potential, a lot of progress has been made in the development of high-throughput experimental screening techniques together with computational methods [2][3][4][5] . In this process, a significant increase of the drug design success rate arises out of the development of novel computational strategies.
Current approaches typically rely on docking assays to either predict or optimize the ligand-binding mode [6][7][8] . Although docking-based methods allow for high-throughput screening of large compound and/or protein mutant libraries, the accuracy of the predictions is limited. The main sources of limitation are the use of simplified energy ("scoring") functions as well as the limited treatment of protein and ligand flexibility and solvation models. In the case of protein engineering studies, only local mutations around the binding/catalytic site can benefit of docking-based methods. To alleviate these shortcomings, molecular dynamics (MD) simulation has become a popular tool in the field of drug design and discovery 9,10 . The MD technique is based on detailed interaction potentials (force fields) and, in principle, includes relevant dynamics of protein, ligand, and solvent provided that enough sampling can be achieved. In some cases, a few direct binding events of ligands can be simulated using brute force MD [11][12][13] , but typically sampling of the relevant degrees of freedom is a limiting factor. A variety of enhanced sampling methods, such as replica exchange, funnel-metadynamics, or adaptive sampling, exist to improve the sampling and to study the binding kinetics and pathways [14][15][16][17] . In addition, when the binding mode is known, calculations can be performed to obtain (relative) binding free energies [18][19][20][21] . In practice, nowadays docking and MD are often combined to increase the accuracy of the former method 4,10,22,23 . Fully atomistic MD simulations, however, are computationally too expensive to allow for high-throughput applications.
A potential solution is the use of coarse-grained (CG) force fields, which reduce the computational cost by uniting groups of atoms into effective interaction sites resulting in a substantial computational speed-up 24,25 . The Martini model 26 is among the most popular CG force fields, and has been applied to study a wide range of biomolecular processes including successful prediction of protein-lipid binding modes 27 . In some cases, binding of lipids to sites deeply buried inside the proteins were recovered by brute force MD 28,29 . Examples of protein-ligand binding with Martini are, however, still scarce [30][31][32] . The question remains, therefore, whether such CG models can be applied to capture protein-ligand binding in general, including small organic compounds, such as enzyme substrates, receptor ligands, drugs, or pesticides.
Here we show that the recently re-parameterized Martini model 33 can perform this task with high accuracy, based on unbiased sampling. For all the simulations presented here, the ligands (depicted together with the CG mapping in Supplementary Fig. 1) were initially positioned randomly in the solvent. We use three different classes of examples to illustrate the potential power of this approach. First, we show millisecond long sampling of the reversible binding of seven different ligands (binders as well as non-binders) to mutants of T4 lysozyme, a well-studied protein for which ample experimental and all-atom MD simulation data are available for comparison. Second, we show spontaneous binding and unbinding of both agonist and antagonist ligands to the adenosine A 2A receptor (A 2A R) and adrenergic β 2 receptor (β2AR), two different members of the membrane-embedded protein family of G protein-coupled receptors (GPCRs), as well as to the nuclear receptor farnesoid X receptor (FXR), all representing prominent pharmacological targets. Third, we demonstrate the versatility of our approach by simulating ligand binding to different enzymes, namely binding of a substrate to a complex catalytic site of an aminotransferase, which is an important representative of biocatalytic applications and binding of drugs to two members of the kinase family, proto-oncogene tyrosineprotein kinase (c-Src) and the AP2-associated protein kinase 1 (AAK1) involved in virus endocytosis and a potential target for COVID-19 treatment 34 .

Results
Binding of ligands to T4 lysozyme. T4 lysozyme serves as a benchmark system to investigate ligand binding 18 . Notably, the L99A mutant is known to bind hydrophobic ligands into a well-defined cavity in the core of the C-terminal domain 35 . To investigate the ligand binding to the L99A mutant of T4 lysozyme, the protein is embedded in a cubic box with 10 nm edge length containing~8850 CG water beads (corresponding to 35,400 water molecules) and one ligand molecule, as shown in Fig. 1a. Thus, the ligand concentration is around 1.6 mM. Figure 1b-e summarizes the results obtained for the seven examined ligands benzene, phenol, indole, thieno-pyridine, toluene, ethylbenzene, and n-propylbenzene with a sampling time of 0.9 ms each (30 independent trajectories of 30 µs per ligand). In addition, the binding of benzene and toluene to the L99A/M102Q double mutant of T4 lysozyme was studied. The additional M102Q mutation reduces the hydrophobic character of the binding pocket 36 which is expected to reduce the affinity to benzene and toluene.
In case of benzene and the L99A T4 lysozyme, each individual simulation of 30 µs shows between 2-9 binding and unbinding events. In total, 156 binding and 147 unbinding events are observed ( Table 1). The crystal structure of L99A T4 lysozyme including a bound benzene molecule (green) is depicted in Fig. 1b. In addition, the figure contains several CG snapshots of benzene taken from our binding simulations (red) and the benzene density in the binding pocket (transparent red isosurface). A comparison of the benzene position in the crystal structure with the CG snapshots as well as the averaged benzene density clearly shows that the CG force field Martini properly locates the buried hydrophobic binding pocket. The experimentally observed binding pose of the ligand in the binding pocket is excellently reproduced with a root mean square deviation (RMSD) of 1.4 ± 0.2 Å, which is similar to the one from atomistic studies 37 . Fig. 1c visualizes the benzene density around L99A T4 lysozyme. The blue, cyan, red, and violet isosurfaces represent occupancies which are 10, 100, 1000, and 10,000 fold higher than in water. A binding channel indicated by the blue isosurface leads from the surface of the protein toward the buried binding pocket. The density inside the binding pocket (red isosurface) is about 1000 fold higher compared to the bulk. Figure 1d unravels the existence of more than just one binding/ unbinding path which could be already easily spotted by eye based on the overall ligand densities shown in Fig. 1c. The benzene density in the vicinity of the binding pocket as well as snapshots of the benzene molecules show four binding/unbinding paths in total. Note that compared to Fig. 1b, c the protein is rotated to enable the distinction of the different paths. The high density channel from Fig. 1c corresponds to path 1 which also has   Fig. 1 Unbiased simulations of ligand binding to T4 lysozyme at the CG Martini level. a Simulation box containing the L99A mutant of T4 lysozyme and one benzene molecule (red) solvated in water (transparent blue surface). b Crystal structure of the L99A T4 lysozyme with benzene (green) in its binding pocket (pdb code: 181L 43 ). In addition, several CG snapshots of benzene (red) and the benzene density in the binding pocket (transparent red isosurface) are shown. The histogram of the RMSD of benzene and the contact protein beads is depicted on the lower right. c Benzene density around L99A T4 lysozyme obtained from averaging 0.9 ms of CG simulations. The blue, cyan, red, and violet isosurfaces correspond to a 10, 100, 1,000, and 10,000 fold higher benzene density than in water. These densities translate to the free energy values shown at the color map. The experimental binding free energy of benzene is between −21.7 kJ/mol 98 and −17.7 kJ/mol 99 . d The benzene density and corresponding snapshots in the vicinity of the binding pocket reveal four binding/unbinding pathways. e Binding free energies calculated from the 0.9 ms of unbiased simulations for all ligands simulated here (blue) in comparison with experimental data (red). The re-parametrized Martini force field can separate non-binders from binders (dashed line). The error bars of the simulated data are <0.7 kJ/mol and not depicted here.  18,44,97,98 . Where a second value is given, this is taken from NMR experiments 99 . d The RMSD is calculated for the binding pocket residues and ligand after aligning the binding pocket residues to the respective crystal structure (for details see Supplementary Methods). e Crystal structures taken from ref. 43 . (benzene, indole) and 45 (toluene, ethylbenzene, and n-propylbenzene). Because no crystal structures are available for non-binders, the experimental binding mode of phenol M102Q 44 , indole 43 , and benzene 43 were used as references the structure to compute the RMSD for phenol, thieno-pyridine, and benzene M102Q, respectively. been found to be the dominant dissociation path in a weighted ensemble approach 38 and the second dominant dissociation pathway in a biased MD study 39 , both at atomistic resolution. The other three paths have been reported previously as well in biased [38][39][40][41] and unbiased 37,42 atomistic MD studies. An example of a binding event is provided in Supplementary Movie 1.
The ligand densities for the other ligands as well as the L99A/ M102Q double mutant are depicted in Supplementary Fig. 2. Also in these eight cases, the hydrophobic binding pocket shows the highest occupancy, in agreement with the crystal structures of T4 lysozyme with the respective ligands [43][44][45] . The binding channel density representing path 1 is very similar to the one of benzene for all ligands except phenol, for which this region is less populated. The experimentally observed binding poses are well reproduced with an average RMSD ≤ 2.1 Å for each of the eight examples (see Table 1, Supplementary Fig. 3, and Supplementary Table 1). In the case of thieno-pyridine, indole, toluene, ethylbenzene, and n-propylbenzene, the RMSD distributions show two binding modes ( Supplementary Fig. 3). Figure 1e compares the binding free energies obtained from our simulations ΔG CG bind with the available experimental binding free energies ΔG exp bind (see also Table 1). The ΔG CG bind values are calculated by integrating the one-dimensional potentials of mean force depicted in Supplementary Fig. 4. Based on the simulated results, it is possible to distinguish between non-binding (three) and binding protein-ligand combinations (six). Moreover, ΔG CG bind is in very good agreement with the experimental values with a mean absolute error of 1 kJ/mol and a maximum error of 2 kJ/mol. To the best of our knowledge, no other study is available to date with such a good agreement of theoretical and experimental binding free energies for this range of T4 lysozyme ligands.
Our results for T4 lysozyme demonstrate the capability to accurately predict the ligand-binding pocket as well as the diverse binding pathways with the recently parametrized Martini 3 force field. Remarkably, as the results were obtained using unbiased simulations, no a priori knowledge of the binding pocket is required. Moreover, a nearly quantitative agreement of the binding free energy for all nine examined systems including different ligands and protein mutants with experimental data is achieved.
Binding to membrane and nuclear receptors. As a second showcase, we present the binding of two small organic molecules to a pharmacologically relevant membrane protein, namely the adenosine A 2A receptor (A 2A R). It belongs to the class A of the GPCRs. The GPCR superfamily is a major therapeutic target whose functioning regulates several physiological processes such as vision, smell, taste, cardiovascular, neurological, and reproductive mechanisms, and A 2A R in particular has been recognized as a drug target for the treatment of many diseases, including cancer, Parkinson's and Alzheimer's disease 46,47 . Here, we examine the ability of the Martini 3 force field to identify the crystallographic binding poses of the endogenous agonist adenosine and the natural antagonist caffeine. The recognition of the A 2A R binding site by the ligands is challenging, since the access to the binding pocket is narrow and regulated by the extracellular loops like extracellular loop 2 (ECL2). In A 2A R, ECL2 contains a small helix that is partially folded over the ligand entry site 48 . Furthermore, interaction of the ligands with either ECL1, ECL2, or ECL3 is relevant to the binding kinetics as proved in different experimental and computational studies 13,49 . In addition to the molecular details of the pocket entrance, the protein environment is more complex, because A 2A R is a transmembrane protein.
Thus, the ligand can also be sorted to the lipid bilayer.
For each system (A 2A R-adenosine and A 2A R-caffeine) we performed 12 simulations with variable number of ligands, spanning from 7 to 13, in a cubic box with 12 nm edge length (Fig. 2a). The binding poses of adenosine and caffeine, after backmapping to all-atom resolution 50 , are compared with the ones in the crystal structures 2YDO 51 and 3RFM 52 , respectively. Figure 2b shows all binding poses of adenosine with an RMSD < 2.6 Å. The RMSD distribution of all observed binding poses is depicted in Fig. 2c (top) and has an average of 3.3 ± 0.5 Å. The best binding pose obtained from the simulations has an RMSD of 2.2 Å (Fig. 2d). Comparing this pose to the crystal structure 51 , we note that adenosine is shifted toward the transmembrane helix 7 (TM7), with a slight tilt of 14°. Notably, the best binding mode closely resembles the experimental one, engaging all the interactions observed in the crystallographic structure (for details see Supplementary Discussion) 51 . In the case of the A 2A Rcaffeine system, we observe 142 binding events, which lead to a well-characterized binding mode. Figure 2c (bottom) shows the RMSD distribution of all binding poses, which has an average of 3.4 ± 0.5 Å. Figure 2e depicts the best binding mode which has a RMSD of 1.9 Å and is remarkably similar to the crystallographic one, being shifted by only 1.4 Å toward TM3 and tilted by 4°( for details see Supplementary Discussion) 52 . The larger number of binding events observed for caffeine with respect to adenosine can be attributed to the planar geometry of the former that better inserts into the gorge forming the entry of the A 2A R binding site. Moreover, the two ligands show a different way to approach to A 2A R from the environment. Caffeine demonstrates a preference to interact with the glycerol and phosphate beads of the lipid bilayer, which is in agreement with experiments and atomistic simulations 53 . On the other hand, this behavior was not observed for adenosine. As a result, caffeine interacts with lipids prior to binding in 87% of the cases, whereas in the case of adenosine this happens in only 40% of the binding events.
A more detailed analysis of our simulations provides structural insight into the ligand-binding mechanism. In the majority of cases-all 15 for adenosine and 89 out of 142 for caffeine-the ligand firstly interacts with either ECL2 or ECL3 before reaching the binding pocket through the passage made by TM3, TM6, and TM7 (Fig. 2d, e). While adenosine directly interacts with ECL2 or ECL3 to approach A 2A R from the water phase, caffeine is also able to establish contacts with the transmembrane helices in the initial binding phase. In such a case, caffeine interacts with ECL3 and TM6 (Fig. 2e, green solid line, 57/142 binding events) or with TM1 and TM7 entering into the cleft formed by the latter two helices (Fig. 2e, black solid line, 51/142 binding events). In 32/142 cases the ligand binds A 2A R through ECL2 (Fig. 2e, red solid line). More details about the adenosine and caffeine binding mechanism to A 2A R are provided in the Supplementary Discussion. An example of an adenosine binding/unbinding event to A 2A R is provided in the Supplementary Movie 2.
As a third showcase, we studied another pharmacologically prominent member of the GPCR superfamily, the adrenergic β 2 receptor (β2AR) 54 . We successfully simulated the binding to β2AR of one natural agonist, adrenaline, and one inverse agonist, propranolol. More details about the β2AR results are given in Supplementary Fig. 5 and Supplementary Discussion.
The fourth case study is the farnesoid X receptor (FXR), a type II nuclear receptor involved in the control of bile acids, triglyceride, cholesterol, and glucose metabolism [55][56][57] . After activation, FXR in complex with the retinoid X receptor (RXR) binds to DNA regulating the expression of proteins involved in bile acid synthesis, triglyceride clearance, cholesterol reduction and modulating insulin sensitivity 56,57 . Therefore, in recent years FXR has become a prominent target for treatment of metabolic disorders, primary biliary cirrhosis, and non-alcoholic steatohepatitis syndrome 58,59 . Here, we investigate the binding of the potent agonist 6-ethyl-chenodeoxycholic acid (also known as obeticholic acid) to FXR. Obeticholic acid is a steroid acid with a specifically bent shape, containing a hydroxyl group and a charged carboxylate in the opposite extremities of the molecule (Fig. 3a). These characteristics make obeticholic acid a challenging ligand for binding studies.
In detail, we simulated 20 µs for 72 replicas of FXR in presence of 4 ligand molecules (Fig. 3a). Overall, we observe 622 binding events. The best bound pose has an RMSD of 1.2 Å and it is remarkably similar to the crystallographic one, being shifted toward helix 12 by only 1.0 Å (Fig. 3b). Figure 3c shows the RMSD values distribution of the identified binding poses computed relative to the crystal ligand-binding mode (pdb code: 1OSV) 60 . Three peaks corresponding to different conformational families are detected at 3.6 Å, 2.7 Å, and 2.0 Å, hereinafter defined as peaks 1, 2, and 3, respectively. The centroid of each family is depicted in Fig. 3d. The presence of three peaks indicates a multistep ligand-binding mechanism. We observe that obeticholic acid first approaches FXR through the loop connecting helices 5 and 6 assuming an external binding conformation (ligand at peak 1 represented as yellow sticks in Fig. 3c), and then upon the receptor rearrangement it reaches an inner binding pose (ligand at peak 2 and 3 represented as ice-blue and orange sticks, respectively, in Fig. 3c). In all the observed binding events, obeticholic acid approaches the binding pocket contacting the loop connecting helices 5 and 6 either from the bottom, passing over helices 2 and 6 ( Fig. 3b black and red arrow), or from the top through helix 1 (Fig. 3b green arrow). Such pathways resemble the so-called paths III, observed for other type II nuclear receptors 61 . In all cases, the ligand points the inner steroidal scaffold toward the binding pocket and the outer carboxylate side chain toward the solvent. This orientation is favored by the salt bridge interaction engaged by the negatively charged carboxylic group of obeticholic acid with several positively charged residues placed around the cleft, namely Arg261, Lys272, Arg328, Lys336, and Arg348 (numbering from pdb code: 1OSV), which are known to be involved in the binding of FXR agonists 62,63 . It is also interesting to note that the loop connecting helices 5 and 6 changes its conformation passing from the apo (unbound) to the holo form (ligand bound) (Fig. 3e). This leads to a reduced solvent accessible surface area (SASA) of the binding pocket computed for the apo FXR if compared with that of the holo form ( Supplementary Fig. 6, blue solid line). The conformational change of the protein is further quantified by the change in the RMSD values of its backbone beads if compared with those of the crystal structure. In particular, the RMSD changes from 4.1 ± 1.0 Å in the unbound conformation to 6.9 ± 0.9 Å during ligand binding and to 5.9 ± 0.8 Å in the final bound state (blue line in  and other type II nuclear receptors 66,67 also indicate stabilization upon ligand binding of the region around the binding site, including helices 2, 5, and 6. Together, our results demonstrate the ability of the reparametrized Martini CG model to reproduce experimental binding modes for a variety of receptors, even in challenging cases as GPCRs where ligand binding occurs in environments of different polarity and the nuclear receptor FXR where moderate conformational rearrangement of the binding pocket is induced by ligand binding. Binding to enzymes. As a fifth showcase, we consider the human proto-oncogene tyrosine-protein kinase (c-Src) which regulates signal transduction in cells 68 . Tyrosine kinase dysfunction is linked to many diseases, and is of particular interest for cancer treatment. Regulation of tyrosine kinases is known to be governed by the presence of multiple allosteric binding sites, which have become the target of many pharmaceuticals 69 . Here, we investigate the binding of the antileukemia drug dasatinib (Fig. 4a) of which the experimental binding mode is known and its binding has been previously simulated at atomistic level 11 . Furthermore, dasatinib as a larger organic molecule compared to the previously investigated ligands, is endowed with four rings with high relative flexibility. It represents a differently challenging system where the ligand not only needs to find the binding site but also orient itself correctly in a specific conformation. The protein was solvated iñ 14,500 CG water beads, and one dasatinib molecule was added in the box, resulting in a substrate concentration of 0.96 mM (Fig. 4a). The chlorobenzene part of dasatinib binds deeply inside the highly conserved ATP pocket, while the piperazine ring remains in the outer part in contact with the solvent. Inverse binding poses, i.e. poses where the piperazine ring is buried in the active site whereas the chlorobenzene ring stays in the solvent, were rarely observed. Using a higher dasatinib concentration (i.e., five dasatinib molecules per box) lead to strong aggregation around the hinge region between the two kinase domain lobules ( Supplementary Fig. 7), although binding to the ATP pocket still occurred. We found an excellent agreement in the RMSD of the bound ligand compared to the crystallographic binding pose (pdb code: 1Y57 70 ), with an average RMSD of 3.0 Å (Fig. 4b), while the RMSD of the atomistic simulation was 2.0 Å 11 . 5 out of 10 independent 30 µs replicas resulted in dasatinib binding with a total of 11 individual binding events observed in the aggregate 300 µs simulation time. The binding events lasted anywhere from a few hundreds of nanoseconds to 5 µs, adding up to 8% of time spent in the bound state. Figure 4b depicts the high-density regions or dasatinib around the protein. One binding pathway was found to be coming from the region between the two protein lobules (Fig. 4b green arrow), in accordance with the sole dasatinib binding path identified in an atomistic simulation 11 . However, in 7 out of 11 binding events dasatinib reached the catalytic site directly from the solvent (Fig. 4b blue arrow). We estimate a k on of 40 s −1 µM −1 , which can be compared to the k on of 1.9 s −1 µM −1 calculated from atomistic simulations 11 and 5 s −1 µM −1 from experiments 71 . Considering the temporal speedup of the Martini model due to the smoother energy landscape (see discussion), our estimate is well in line with the expected values.
As a sixth showcase, we studied the S-selective aminotransferase of Vibrio fluvialis which stereoselectively catalyzes the transfer of an amino group from a donor to an acceptor ketone 72 . Using unbiased simulations we obtain good agreement of the binding pose of the substrate acetophenone. More details are provided in the Supplementary Fig. 8 and Supplementary Discussion. The seventh and last case study presented here targets a serine/threonine kinase involved in the regulation of cellular virus endocytosis, namely the AP2-associated protein kinase 1 (AAK1) 73,74 . Targeting AAK1 has been shown to successfully prevent viral infection of cells 74 and the compound baricitinib was recently suggested as potential drug to prevent endocytosis of SARS-CoV-2 by inhibiting AAK1 34 . Baricitinib is a known inhibitor of the Janus kinases 1 and 2 JAK1/2 and is applied e.g. against rheumatoid arthritis 75 . Due to the lack of a crystal structure of AAK1 in complex with baricitinib, we used the crystal structure of AAK1 bound to the inhibitor LKB1 (pdb code: 5L4Q 76 ) to generate the CG protein model, while the structure of the homolog BMP-2-inducible kinase bound to baricitinib (pdb code: 4W9X 77 ) was employed as reference model for the baricitinib binding mode. Overall, the BMP-2-inducible kinase has 72% of identity with AAK1. Moreover, the catalytic site is fully conserved (100% identity), thus representing a good reference for the baricitinb binding mode to AAK1. Figure 4c shows the simulation box containing AAK1 and one molecule of baricitinib. We performed 30 simulations of 30 µs each and observed 4 binding events to the buried catalytic site (red surface, Fig. 4d). During one of these binding events, the ligand stayed for the last 15 µs of the replica in the catalytic site. Despite the complexity of the system, we obtain an RMSD for the main binding mode of 2.6 Å with respect to the crystal pose (Fig. 4d). Two additional binding conformations are also found at farther distance from the binding site. Moreover, the ligand exhibits a tendency to bind to a more external site next to the catalytic pocket that seems to be involved in the association pathway. This pathway seems to resemble one of the pathways observed for c-Src, which may be a common feature for kinases. Several additional low-affinity pockets (gray surfaces) on the protein surface are identified, which might be relevant for the activity of baricitinib.
Together, the results presented on binding of drugs and substrates to different classes of enzymes further show the potential of our approach to predict accurate binding modes and pathways for large, flexible, and complex ligands.

Discussion
Applications such as structure-based drug design are particularly challenging for CG modeling because of several requirements: (1) high chemical specificity of pocket-ligand interactions; (2) capability to represent all possible components of the system (as proteins, cofactors, drug candidates, solvent, lipids, etc.) in a coherent way; (3) realistic representation of conformational flexibility of each molecule in the system; (4) accurate thermodynamics and kinetics of binding. Currently, none of the CG force fields available fulfills all the requirements listed above. The examples showed in this work indicate that the current state of the Martini model seems to finally achieve most of these requirements with reasonable accuracy in relation to atomistic models. The key improvement in Martini 3 33 to enable such applications is the enhanced packing of the CG beads, achieved by re-balancing of the cross-interactions of different bead sizes 78 , as well as by the re-parametrization of bonded distances based on molecular volume and shape. As a result, protein cavities are represented more realistically and ligands can fit better. Besides, the expansion of bead chemical type options in Martini 3 allows for a better coverage of the chemical space, facilitating CG modeling of rather complex small molecules such as baricitinib and dasatinib.
Given the improvements in accuracy in Martini 3, another key advantage of CG models in general is their computational performance. For instance, Martini CG based docking of biomolecular complexes can be around one order of magnitude faster compared to atomistic models, as recently demonstrated for the Haddock program 79,80 . Benchmarks tests performed with the program package Gromacs (version 2018) 81 showed that Martini based MD simulations of protein-ligand systems can be 110-350 times faster than all-atom simulations, with the performance gain increasing with growing system size (see Supplementary  Discussion and Supplementary Table 2). Considering diffusioncontrolled processes such as ligand binding, the smoother potential landscape of Martini can also provide 2-3 times faster association/dissociation to the proteins (Supplementary Table 3). In addition, faster protein dynamics can play an important role in induced-fit processes, as exemplified by the binding of the obeticholic acid to FXR and the dasatinib binding to c-Src kinase. The latter shows a~8 times speed-up in association rate constant compared to atomistic simulations. Another performance advantage can be achieved by coarse-graining in the chemical space, as recently demonstrated by Menichetti et al. 82 . Because certain Martini CG moieties can represent more than one chemical fragment at the same time, virtual screening using fragment-based strategies could lead to an improvement in performance of three to four orders of magnitude. The combination of increased performance, smother potential surfaces, and coarsegraining of the chemical space can bring virtual screening protocols to speed-ups in the order of 10 5 -10 7 in comparison to approaches based on atomistic models.
The main limitation to achieve such a high-throughput screening pipeline is the currently limited set of available ligand parameters for Martini. The development of a curated and validated database for Martini CG ligands is therefore of paramount importance. Such ligand databases in combination with automated tools to generate CG models 83,84 (https://github.com/ marrink-lab/cartographer) can expand the accessible chemical space of Martini to millions of compounds. Other aspects that can limit the accuracy of the model in certain applications, and that should be kept in mind, are: (1) the poor representation of protein flexibility by elastic network models; (2) the hydration of pockets by water molecules that cannot be represented by CG water models; (3) limited accuracy at the CG level to differentiate enantiomers or to fully represent directionality in binding poses; (4) the approximate nature of binding kinetics. Most of these issues are inherent to the process of coarse-graining, but some of these problems can be alleviated at least to some extent. For instance, protein flexibility can be greatly improved by the combination of Martini and Gō-like models 85,86 . Hydration of pockets can be modeled by the usage of smaller water beads, as it was already applied here in the A 2A R and β2AR case (see Supplementary Methods). Polarizable Martini models 87,88 can be applied to cases for which directionally of hydrogen bonds or an improved description of the electrostatic interactions are necessary. Differentiation of enantiomers is a clear challenge for CG models, and might require two-state CG models 89 or the use of multi-resolution tools to couple Martini to all-atom models 90,91 . Further refinement of the ligand pose and the binding pocket can be achieved by back-mapping to atomistic resolution 50 as it was already used in the present study in the case of the GPCRs and FXR. Concerning kinetics, the reduction of friction from the missing atomistic degrees of freedom causes a natural speed-up of the dynamics (Supplementary Table 3). This implies that only order of magnitude estimates can be provided when using a CG model such as Martini. However, trends for ligands with the same level of resolutions can be expected to be captured well. In the case of binding, diffusional encounter between the ligand and the binding pocket entrance largely determines the binding kinetics. Because the molecular shape of the ligand, the protein, and the binding pocket entrance are represented well in Martini 3, trends in binding kinetics are expected to be represented reasonably. Moreover, realizing that in particular unbinding rates of ligands from proteins often involve large free energy barriers, the overall accuracy of kinetic estimates for unbinding will mostly rely on a careful representation of the barrier energetics. Here, Martini is expected to perform well as the model heavily relies on reproducing free energy data and the correct kinetics rates might be retrieved from a rigorous estimate of the friction reduction 92 .
In summary, we have demonstrated that the Martini 3 force field can be used to simulate protein-ligand binding in a brute force approach. We have illustrated its capability by spanning a range of systems from the well-characterized model system T4 lysozyme over pharmacologically relevant receptors such as the A 2A R, β2AR, and FXR, to a number of different enzymes, namely the aminotransferase Vf-ATA and the kinases c-Src and AAK1. In the future, the computational performance can be straightforwardly increased by optimizing the ligand concentration as well as by employing enhanced sampling techniques. Moreover, our results pave the way to an efficient computational approach to quantitatively predict binding thermodynamics and potentially capture trends in kinetics. Because no a priori knowledge of the binding pocket is required and a known pocket does not influence our approach in any way, it entails the possibility of finding new additional pockets. In this view, we could envision computational competitive binding assays. Together with ongoing developments of an automated topology builder for Martini, the presented results open the way to a high-throughput screening pipeline, potentially screening millions of drugs and protein mutants. Protein CG models. The bonded parameters of the protein models were slightly adapted from the standard Martini 2.2 settings including the recently suggested side-chain corrections 93 , applied not only for β-strands but to all loops and secondary structure elements. An elastic network comparable to the one of the Martini 2.2 protein model 88 was used to maintain the secondary and ternary protein structure without exclusions of the non-bonded interactions between the backbone beads connected by the elastic network 78 . See more details about the protein models in the Supplementary Methods.
Generic ligand parametrization. Models for the ligands were parametrized according to the Martini (3.0) procedure: First, mappings were designed based on the following principles: (i) minimize number of CG beads used and use regular-(R-), small-(S-), and tiny-(T-) beads for 4-to-1, 3-to-1, and 2-to-1 (non-hydrogen) atoms-to-CG-site mappings; (ii) describe aromatic rings by T-beads; (iii) take into account the symmetry of the molecule. Secondly, Martini bead types were assigned based on the chemical building block they are taken to represent. Bonded interactions are then obtained based on atomistic models. Note that we used center of geometry (COG)-based mapping taking into account also the hydrogen atoms. COG-based mapping leads to better molecular (e.g., solvent accessible surface area, SASA) and bulk (e.g., mass densities) properties for the models 33,94 . Finally, models are validated using partitioning free energy data and comparison to atomistic SASA values. More details about the ligands models are given in the Supplementary Methods.

RMSD calculation.
To calculate the RMSD between the simulated binding poses and the crystal structure, the binding pocket was aligned to the CG crystal structure. The list of residues used for the alignment of the various systems is given in the Supplementary Methods. Finally, before the RMSD calculation, the CG ligand in the crystal structure, which was obtained by transforming the bound atomistic ligand to its CG resolution, was minimized for one step to account for slight changes in the bonded parameters. Note that in case of ligands with high symmetry like benzene, all possible orientations have to be taken into account. The lowest RMSD of all possible orientations is the correct value because all other values are too high due to flipping or rotation of the structure which does not change the chemical structure of the binding pose.
Binding free energy calculation. Supplementary Fig. 4 shows the radial ligandreceptor potentials of mean force (PMFs) obtained from unbiased MD simulations. One dimensional PMFs were computed based on normalized distance distributions (p(r)) between the center of geometry of the pocket and of the ligand, including volumetric (also called entropic) correction 2 × ln(r) 95,96 .
PMFðrÞ ¼ ÀRTln pðrÞ ð Þþ2 lnðrÞ ð 1Þ The distance r was shifted 0.25 nm to take into account the average radial volume of the ligands, which was estimated based on their average SASA values. Besides, the entire PMF was shifted to zero at 4.5 nm distance. The binding free energies (displayed in Table 1) were estimated by integrating the PMF over the distance r in terms of Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data supporting the findings of this paper are available from the corresponding authors upon reasonable request.

Code availability
Force-field parameters and procedures (e.g. tutorials) are publicly available at http:// cgmartini.nl.