Large-scale integrated super-computing platform for next generation virtual drug discovery.

Traditional drug discovery starts by experimentally screening chemical libraries to find hit compounds that bind to protein targets, modulating their activity. Subsequent rounds of iterative chemical derivitization and rescreening are conducted to enhance the potency, selectivity, and pharmacological properties of hit compounds. Although computational docking of ligands to targets has been used to augment the empirical discovery process, its historical effectiveness has been limited because of the poor correlation of ligand dock scores and experimentally determined binding constants. Recent progress in super-computing, coupled to theoretical insights, allows the calculation of the Gibbs free energy, and therefore accurate binding constants, for usually large ligand-receptor systems. This advance extends the potential of virtual drug discovery. A specific embodiment of the technology, integrating de novo, abstract fragment based drug design, sophisticated molecular simulation, and the ability to calculate thermodynamic binding constants with unprecedented accuracy, are discussed.

Introduction

Drug discovery is distressed. The number of approved new molecular entities has declined steadily for 15 years [1], the cost per new approved compound has breached the one billion United States Dollar (USD) benchmark [2] and, by one informed estimate, the financial return provided by all therapeutic product categories does not even recover the capital costs of their development. Moreover, when biologics are removed from this model, the net present value of pharmaceutical Research and Development (R&D) investment is actually negative [3]. Disturbingly, poor research performance occurs in the context of increasing social need. The rising incidence of drug resistant bacteria is well documented [4], viral pathogens such as Severe Acute Respiratory Syndrome (SARS), Dengue and H1N1 pose pandemic threats [5], and an aging and more affluent global population drives up the prevalence of chronic diseases that are if anything more difficult than infectious diseases to drug. In the United States, where by 2023 the Census Bureau projects ∼15% of males and ∼19% of females will be 65 years old or older, rates of cancer are expected to rise from 33% to 62%, of diabetes from 33% to 53%, of cardiovascular complaints from 6% to 39%, of mental disorders from 35% to 54%, for a total increase in chronic morbidity in the U.S. population from 17% to 42% [6]. An increase in successful drug discovery, especially against difficult chronic targets, is clearly desirable. However, progress must materialize in the framework of reduced per compound cost and enhanced efficiency, since capital will not continue to flow into a sector that offers net negative returns. Computational or ‘virtual’ drug discovery strategies, potentially cheaper and faster, offer attractive alternative, or at least complimentary, routes to improved R&D performance in the therapeutics sector [7].

Computational drug discovery

The physiological effect of a drug is mediated by electrostatic and geometrical interactions of the atoms of the ligand with the atoms of its corresponding receptor, interactions which conform to the laws of physics and quantum chemistry, and which can therefore be described by predictive mathematical models [8]. Although these models are complex (and in the quantum case inherently non-exact), researchers active in the computer intensive field of molecular graphics realized thirty years ago that in silico assessment of drug-receptor binding could be deployed to accelerate drug discovery [9, 10]. They also realized that continuing operation of Moore's Law, with 18 month doublings in computing efficiency and economy, implied an on-going improvement and refinement in computational techniques. The crux of the computational drug discovery paradigm is this coupling of the fundamental laws of biophysics with the accelerating technological performance of the semi-conductor industry. The former inspires faith that the approach can work in theory, the latter that it will work in practice.

Virtual compound screening

Virtual screening, whether of compounds or molecular fragments, has two stages. First, the algorithms attempt to find the correct conformation and position the ligand in the active site of the receptor, and then they try to quantify the quality of particular atomic arrangements by assigning a score. Several technically different approaches to predicting ligand–receptor interactions have been developed, but all are known as ‘docking’ algorithms after the suggestively named primogenitor program, ‘DOCK’ [9]. The modeling of ligand–receptor atomic interactions presupposes an accurate three-dimensional molecular structure of the receptor so that inter-atomic forces can be calculated. Since protein folding cannot yet be modeled, this means having an X-ray crystal or NMR structure of the receptor, or a homology model which maps a related protein sequence onto a known structure. A priori one might suppose the experimentally determined crystal structure to be inevitably superior. Surprisingly, a meta-analysis of DOCKing studies concluded that in some cases virtual screening was more successful on homology models compared to experimental structures [11]. This seems counter-intuitive, but may indicate that the relaxed precision of the homology models indirectly capture conformational flexibility that is lost in ‘frozen’, possibly subtly distorted crystal structures. In any case, an homology model must start from a closely related experimental structure, so an important contributing factor in the increased utility of computational drug discovery is the rapid growth in the number of available protein structures (currently approaching 75,000 structures in the Protein Data Base http://www.pdb.org/pdb/statistics/contentGrowthChart.do?content=total&seqid=100), a number which in turn reflects improvements in protein production, robotic crystallization regimens, and the wide availability of sophisticated advanced light sources [12]. Another positive development in virtual screening infrastructure is the creation of curated virtual compound databases that provide large prebuilt sets of virtual representations of commercially available molecules suitable for input to virtual screens. ZINC at the University of California San Francisco [13], and EDULISS at Edinburgh University [14], are two examples. The much larger Chemical Universe Database GDB-13 takes a different approach, attempting to construct the universe of ‘synthetically plausible’, rather than ‘available’ compounds [15].

Molecular docking: successes and limitations

At a high level the performance of dock programs can be measured by two criteria: ‘DOCKing power’ (the ability to identify the correct experimental ligand binding pose in a collection of incorrect, computer generated ‘decoy’ poses, i.e. the ability to correctly position ligands in the active site, or to ‘pose’ them); and ‘scoring power’ (the ability to produce dock binding scores that correlate with experimentally determined binding affinities). In the past decade a large number of comparative studies of the performance of various dock programs have been undertaken, in both academic and pharmaceutical settings [11••, 16, 17••, 18•, 19, 20, 21, 22]. Despite the diverse backgrounds of the investigators, and although these studies differ in methodology and are not directly comparable, they nevertheless unanimously agree on two points. One, dock algorithms fairly accurately pose ligands in the active site, and two, the same dock algorithms poorly score those ligands’ affinity. In other words, dock programs correctly identify the geometry of ligand–receptor systems, but, do not in general accurately predict the binding energy, and therefore cannot predict ligand potency. To make this concrete, a typical dock screen might produce 1000 ‘hit’ compounds, but, the most potent compounds are as likely to be ranked at the bottom of that list by the scoring function as they are to appear near the top. This is a significant deficiency since the expected potency of a compound often will be the operational feature of interest, for example in prioritizing compounds for medicinal chemistry. In sum, Dock algorithms can ‘pose’ ligands well but they ‘score’ them poorly.

Beyond dock scores: accurate binding affinity from thermodynamic calculation with MAPLE CAFEE

ΔG  =  RT  ln  K d exactly relates the computed Gibbs free energy difference ΔG and experimentally measured dissociative constant K d under temperature T (where, R is the gas constant). Free energy differences between bound and unbound equilibrium states of a protein–ligand–water system (ΔG) gives the binding affinity of the ligand, which in general translates into drug efficacy. In other words, correct computation of ΔG values for a series of ligands leads immediately to a ligand list accurately ranked by potency.

Computational methods to perform the ΔG calculation have been studied enthusiastically since the late 1990s when it was proved that a nonequilibrium process in finite-time can derive the binding free energy exactly [23, 24, 25]. This theoretical insight was followed up in 2005 when ΔG was shown to be approachable by massively parallel computation. Scaling up to thousands of concurrent CPUs reduced the computational requirement from years to days, and allowed binding free energy of real molecular systems to be computed. Access to this high performance computational resource made possible an important series of proof of concept experiments that in turn produced calculated binding affinities in excellent agreement with corresponding experimental values [26, 27]. Subsequently, the computational methodology has been improved, a better force field refining method has been implemented, and the platform, christened Massively Parallel Computation of Absolute binding Free Energy with well Equilibrated system (MAPLE CAFEE), has been further validated in follow-on projects [28, 29, 30]. Useful thermodynamic calculation comes with a heavy computational burden (an example of the trade-offs between speed and accuracy that turn up in many computational problems), but recent advances in the construction of large-scale computing environments have shrunk the envelope of computational time and brought these methods into the realm of practical application ([31], http://www.nsc.riken.jp/index-eng.html, http://www.fujitsu.com/global/news/pr/archives/month/2010/20100928-01.html). Cumulatively, improved computational methodology, enhanced infrastructure, and theoretical advances have combined to achieve chemical thermodynamic calculations that are sufficiently accurate to provide reliably ordered and prioritized lists of hit compounds, either from virtual libraries or de novo design.

Comparison of binding affinity calculation by four methods

In 2008, three separate organizations compared binding affinity prediction methods on a defined data set using identical force field parameters [32]. The system of interest was FK Binding Protein (FKBP)–ligand–water of about 17,000 atoms, with water explicitly modeled with TIP3P. A co-crystal structure of the FKBP-FK506 complex was obtained from the Protein Data Bank (PDB). Then other ligands replaced FK506 one by one. The four methods tested were Fragment Molecular Orbital method (FMO), Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA), the hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) and MAPLE CAFEE.

Figure 1 a–d, reprinted from [32], shows the correlation between observed binding affinity and computed binding free energy for the four methods. While the computational time required for each method conformed to expectation, and each computational method produced results that correlate somewhat with experimental values, there were significant differences in the scale and shifts from the origin of the vertical axis. MAPLE CAFEE produced the best agreement with a difference from observed affinities within 0.5 [kcal/mol] (an error of 1.4 [kcal/mol] is a good benchmark for computational accuracy since the error of experimentally measured affinity is in this range).

Figure 1

Controlled comparison of four binding affinity prediction methods on the same 17,000 atom FK506-FKBP data set. Each graph shows the correlation between observed (experimental) binding affinity (on the horizontal axis) and computed binding free energy (on the vertical axis). The dotted lines show a slope of 1 (perfect correlation, with possible offset). The yellow box shows the area with no X–Y offset. (a) Massively Parallel Computation of Absolute binding Free Energy with Well Equilibrated system (MAPLE CAFEE). (b) Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA). (c) Quantum Mechanics/Molecular Mechanics (QM/MM). (d) Fragment Molecular Orbital Method (FMO).

A blind challenge of MAPLE CAFEE

In 2009, MAPLE CAFFE was challenged by a pharmaceutical industry collaborator to a blind-test of binding affinity prediction (unpublished results). The algorithm was given the structure of a cancer related protein and five ligands, one a co-crystal structure, three without co-crystal data but where the small molecule was similar to the first, and a final ligand without co-crystal data and structurally dissimilar to the others. The experimentally determined binding affinities were hidden. As seen in Figure 2 the results showed very good agreement between binding free energies computed by MAPLE CAFEE, and experimentally measured dissociation constants. There was one outlier, however, the pharmaceutical partner suggested this particular experimental measurement was unreliable, and indeed, when repeated, the experimental value fell within 1.4 kcal/mol of the computed value.

Figure 2

Blind-test of MAPLE CAFEE. The binding energies of five ligands on a cancer related protein were calculated by MAPLE CAFÉ, then plotted against the experimentally determined binding energies. The dotted lines define a band within 1.5 kcal of perfect correlation. The binding constant of the outlier was measured a second time in a repeat experiment. The recalculated value is shown by the red arrow.

A coherent structure based drug-design platform

MAPLE CAFEE is part of a larger Structure and Simulation Based Drug-Design (SSBDD) platform (Figure 3 ), that also includes the Optimum Packing of Molecular Fragments (OPMF) module, an abstract fragment based, de novo, drug-design tool (Figure 3a–e). Each component of SSBDD can be applied to specific purposes independently, but the platform is designed as an integrated system for creating novel and active chemical entities.

Figure 3

The OPMF and MAPLE CAFEE workflows. The standard drug-design process of OPMF comprises five steps. (a) Search energetically stable positions for each abstract fragment in/on a targeted protein with molecular mechanics simulation (MM). At the discretion of the investigator, physically discovered molecular fragments, for example from NMR or crystal soaks, can be substituted for the virtual fragment inputs. (b) Select a stably positioned set of abstract fragments following drug-design strategies. (c) Exhaustively connect the selected fragments with constraints to generate as many abstract molecular skeletons as possible, then stabilize them with MM. (d) Assign possible real-atoms to the abstract skeletons considering biophysical constraints and synthetic feasibility. (e) Screen out unfavorable structures using heuristics such as the occurrence of pharmacophores that raise toxicological ‘red flags’. The standard MAPLE CAFEE workflow comprises five steps (f)–(j). (f) Setup molecular models for protein, ligand, and water and assign the force field parameters with Force Field Formulator for Organic Molecules (FF-FOM) [28]. When no co-crystal is available, reasonable initial biding poses should be obtained with docking or other methods. It is notable that having realistic molecular models with refined force field parameters is essential as well as efficient algorithms with reliable computational method for estimation. (g) Equilibrate the system in the initial pose with a MD job. (h) Compute the nonequilibrium work term between bound and unbound equilibrium states. Actually, the two end states are divided into microstates with decoupling parameters, and the trajectories for each microstate are computed with MD jobs. (i) Estimate the most likely representative binding free energies between microstates with the BAR (Bennett Acceptance Ratio) statistical method [34]. (j). Finally, sum up all the micro binding free energies to get the total binding free energy. Molecules can be submitted to MAPLE CAFEE for free energy calculation from any source. For example, they might originate as proposed molecules in a medicinal chemistry series (k), or from a DOCK type virtual screen (l).

OPMF is a virtual fragment based de novo drug-design tool that generates chemical structures predicted to bind and modulate the activity of target proteins whose 3D structures are known. In order to examine the vast virtual chemical space efficiently, OPMF generalizes real-atom fragments to abstract conceptual fragments at the first step. This strategy allows a small number of abstract fragments to represent thousands of real fragments. For example, benzene, pyridine, pyrimidine, triazine are generalized to an abstract six-membered aromatic ring built of abstract atoms having representative diameter and properties. The approximation is reversed later when abstract structures are replaced by precisely assigning real-atoms that rigorously capture exact physical features including electrostatics, orbitals and van der Waals radii.

MAPLE CAFEE is a binding affinity prediction system, as already mentioned, that computes ΔG as precisely as possible using a massively parallel sampling method reinforced by many molecular dynamic (MD) simulations [26] (Figure 2f–j). Parallelization, via the Message Passing Interface (MPI) plays a critical role in the system equilibration step of MAPLE CAFEE (Figure 3g), for example by allowing information such as hydrogen bond frequency to be derived and fed back into the drug-design process. Computing the nonequilibrium work term (Figure 3h) invokes numerous MD jobs, and this is made practical by adopting a massively parallel approach. Each MD job generally takes a couple of days on a single CPU core, but MPI parallelism for a microstate MD job would reduce the computational time almost in proportion to the number of cores. Therefore, if 10,000 CPU cores were available, it would be possible to examine hundreds of novel structures in a week without any actual synthesis. This is the biggest advantage of in silico technologies. As an alternative to generic parallelization on commodity clusters, it has recently been shown that specialized, purposed built computer architecture can accelerate MD calculations by two orders of magnitude [33].

Conclusion and future works

An integrated SSBDD platform has been built to create active novel chemical structures for chemical synthesis. A novel feature of the platform is its ability to calculate ligand affinity constants with great accuracy. Despite the substantial computational resources required to achieve this accuracy, the approach is cost effective because it dramatically reduces the time and expense associated with unnecessary experimental synthesis and assays, enabling chemists to focus their synthetic efforts chemical structures that are most likely to lead to therapeutic success. Currently, binding affinity prediction with MAPLE CAFEE tends to be unstable when the actual affinity is very weak or ligands are protonated. We are now improving the sampling algorithms for weak binding state, and charge correction methods for protonation. We are also reducing the requisite CPU resources as much as possible without sacrificing accuracy.

The availability of next generation super-computers puts a practical thermodynamic guided drug-design platform, both accurate and swift, within reach. Moving forward, accelerated computation, and the virtual drug discovery platforms it will support, will bolster successful therapeutic R&D outcomes and contribute to better human health.

Conflict of interest statement

SM is General Manager of the in silico Drug Discovery Research Division, Bio-IT Business Development Unit, Fujitsu Limited, which has a commercial interest in MAPLE CAFEE, OPMF, and the SSBDD PLATFORM. WM collaborates with the Fujitsu group but has no financial interest. WM is supported by a block grant to the Experimental Therapeutics Centre from the Biomedical Research Council of the Agency for Science and Technology Research (A*STAR), Ministry of Trade and Industry, Singapore.

References and recommended reading

Papers of particular interest, published within the period of review, have been highlighted as:

• of special interest

•• of outstanding interest