Method and computer program product for drug discovery using weighted grand canonical metropolis Monte Carlo sampling

Abstract

A method and computer program product for modeling a system that includes a protein and a plurality of different fragment types in order to identify drug leads is presented. The basis of the method is a weighted Metropolis Monte Carlo approach for sampling the grand canonical ensemble. This method distinguishes itself from an energy minimization approach in that it provides fragment distributions which are consistent with thermal fluctuations at physiologically relevant temperatures. The weighted Metropolis Monte Carlo scheme performs a quasi-uniform sampling of all regions of interest on the protein, and, in this way, enables to resolve the wide range in densities of the thermodynamic distribution which could not be achieved by a non-weighted Metropolis scheme. Making use of the properties of the grand canonical ensemble, the affinity of fragments for different regions on the protein surface can be efficiently computed, using a so-called “simulated annealing of the chemical potential” process. A protein binding site is then identified as a region with high affinity for multiple fragments with a diverse set of physico-chemical properties. Within a binding site, assembly of fragments into drug leads is finally carried out based on binding affinity of the different fragments, on geometric proximity, and a variety of rules by which organic fragments may bond together.

Description

BACKGROUND OF THE INVENTION

[0002] 1. Field of the Invention

[0003] The invention described herein relates to models for molecular interaction, and in particular the use of such models for drug discovery.

[0004] 2. Related Art

[0005] In determining drug leads, it is often desirable to model a system that includes a protein and a set of small molecular fragments. Given the three dimensional structure of a target protein, usually obtained experimentally from x-ray crystallography, the basic interactions between the protein and the small fragments (typical average molecular weight of 150) are computed. This computation can be carried out by Monte Carlo (MC)-type modeling and analysis (usually implemented in software) for a large collection of organic fragments with diverse physico-chemical properties. For such a fragment-based approach, the number of considered fragments can in practice be in the hundreds to thousands. What are needed, therefore, are an efficient method and computer program product for modeling such a system of fragments for purposes of determining drug leads.

SUMMARY OF THE INVENTION

[0006] The invention described herein includes a method and computer program product for modeling a system that comprises a protein and a plurality of fragments in order to identify drug leads. To analyze the interaction between a given type of fragment and a protein, the states of the fragment with respect to the protein are sampled from a thermodynamically relevant distribution. The underlying sampling algorithm is a Weighted Grand Canonical Metropolis Monte Carlo approach, referred to herein as WGCMMC.

[0007] The purpose of this weighted approach is to enable an essentially uniform numerical sampling of all states of interest of the fragment with respect to the protein, i.e. sampling deeper and shallower energy wells with the same thoroughness, while still avoiding the sampling of very unfavorable poses (typically resulting from steric clashes). The data is then appropriately re-weighted for the sampling to correctly represent the considered thermodynamic ensemble.

[0008] Saving a state of a system described by a grand canonical ensemble comprises saving the states of all fragments currently present in the system. In turn, saving a fragment state comprises storing its position, orientation, potential energy and weight. Note that in the framework of the grand canonical ensemble, the number of fragments in the system fluctuates from one system state to another.

[0009] By making use of this fragment data, binding modes can be identified and corresponding binding free energies estimated. The fact that the simulation system is considered in the framework of the grand canonical ensemble, instead of the canonical ensemble, enables through simulated annealing of the chemical potential an efficient estimation of the free energy of binding of a given fragment type for various binding modes on the protein surface. This binding data for the different fragment types can then in turn be used for identifying the relevant protein binding sites, and for assembling the different fragment types to obtain larger ligand molecules.

[0010] In practice, the weighting procedure is implemented by subdividing the sampling space with a grid. Each grid cell center x is assigned a local, numerical chemical potential field value Bnum,(x), which is adapted iteratively during the computation, based on preceding sampling statistics, so as to ultimately ensure an approximately uniform numerical sampling of fragment states at all regions of interest around the protein. Bnum is related to the energetic cost of inserting or removing a fragment from the numerical distribution in the cell centered at x, and the difference between its local value Bnum (x) and the actual physical chemical potential B of the system defines the weight w for each sampled fragment state.

[0011] Once the Bnum field has sufficiently converged, as a result of successive iterations, and the Markov chain associated with the Metropolis algorithm has equilibrated, the actual Monte Carlo sampling can be gathered. This is carried out by keeping the Bnum field fixed, and then periodically saving the state of the system along the Markov chain. The number of Markov steps interspacing the gathered states must naturally be sufficiently large to ensure proper decorrelation.

[0012] Further embodiments, features, and advantages of the present inventions, as well as the structure and operation of the various embodiments of the present invention, are described in detail below with reference to the accompanying drawings.

DESCRIPTION OF THE FIGURES

[0013] FIG. 1 is a flowchart illustrating overall processing of an embodiment of the invention.

[0014] FIG. 2 is a flowchart illustrating the initial step of preparing a molecular model for the system to be analyzed.

[0015] FIG. 3 is a flowchart illustrating the modeling process at the systemic level for computing the fragment-protein interactions using a Weighted Grand Canonical Metropolis Monte Carlo (WGCMMC) approach, according to an embodiment of the invention.

[0016] FIG. 4 is a flowchart illustrating the convergence phase of the simulation system, according to an embodiment of the invention.

[0017] FIG. 5 is a flowchart illustrating the sampling phase of the simulation system, according to an embodiment of the invention.

[0018] FIG. 6 is a flowchart illustrating the process of identifying potential binding sites, according to an embodiment of the invention.

[0019] FIG. 7 is a flowchart illustrating the process of clumping fragments before assembly into drug leads, according to an embodiment of the invention.

[0020] FIG. 8 is a block diagram illustrating a computing platform on which a software embodiment of the invention can be stored and executed.

DETAILED DESCRIPTION OF THE INVENTION

[0021] A preferred embodiment of the present invention is now described with reference to the figures, where like reference numbers indicate identical or functionally similar elements. Also in the figures, the left-most digit of each reference number corresponds to the figure in which the reference number is first used. While specific configurations and arrangements are discussed, it should be understood that this is done for illustrative purposes only. A person skilled in the relevant art will recognize that other configurations and arrangements can be used without departing from the spirit and scope of the invention. It will be apparent to a person skilled in the relevant art that this invention can also be employed in a variety of other devices and applications.

[0022] I. Overview

[0023] The invention described herein is a fragment-based approach for designing drug leads. For this purpose, Locus Pharmaceuticals, Inc., Blue Bell, Pa., developed the Locus Monte Carlo (LMC) code. The approach described herein makes use of a Weighted Grand Canonical Metropolis Monte Carlo (WGCMMC) algorithm for sampling fragment states around the target protein, of a given fragment type. This sampling data can then be directly used for estimating the free energy of binding for different binding modes of the given fragment type on the protein surface. This computation can be carried out simultaneously for hundreds to thousands of different fragment types on a computing platform consisting of multiple processors, such as a PC cluster. This approach mainly distinguishes itself from a similar process implemented by Mezei and Guamieri in their Metropolis Monte Carlo (MMC) code [Guamieri, F. and Mezei, M., J. Am. Chem. Soc. 118:8493-8494 (1996)], in that it removes fragment-fragment interactions.

[0024] During the Monte Carlo sampling, a set of attributes are saved for each rigid fragment instance, including the coordinates (x,y,z) of the fragment's center of mass, the unit quatemion q=(q1, q2, q3, q4) characterizing its orientation, and the potential energy of interaction E between the fragment and the protein.

[0025] This LMC data for the different fragment types can be analyzed for identifying potential binding sites using diagnostic tools such as the Locus Cluster Analysis (LCA) code and the Locus Binding Analysis (LBA) code (Locus Pharmaceuticals, Inc., Blue Bell, Pa.). These tools are based on the postulate that a protein binding site must be a localized high affinity region for a diverse collection of fragments, i.e. fragments with different physico-chemical properties. It is indeed assumed, that diverse interactions in a localized region are the necessary condition for ensuring the specificity of a binding site. If available, one naturally also makes use of experimental binding site data (e.g., co-crystal X-ray data and residue mutational analysis) in determining the final site within which the leads are designed.

[0026] Within the chosen binding site, fragments can be assembled into the actual candidate drug leads, usually composed of four to five fragments and thus having a molecular weight of the order of 600-800, using a software package such as the Locus Chemistry Design (LCD) software (Locus Pharmaceuticals, Inc., Blue Bell Pa.). Here again, use is made of the LMC fragment data in providing preferred fragment states—positions and orientations—with respect to the protein (also called fragment poses). Assembly of fragments is carried out based on geometric proximity, and using a variety of rules by which organic fragments may bond together. In somewhat more detail, two fragments can be assembled, if the relative positions of their atoms enable, within given tolerances, to establish a certain type of covalent bond, with specific bond lengths and angles. The most elementary example of bonding rule is of the form

—CH2—H+H—CH2——CH2—CH2

[0027] Other bonding rules, such as the fusing of methyl groups or merging of cyclic rings, for example, may also be considered.

[0028] Fragment-based computational approaches are well-known. One example is the Multiple Copy Simultaneous Search (MCSS) numerical tool presently commercialized by Accelrys, of San Diego, Calif., which derives from an original version developed by the group of Karplus, Harvard University, MA, [Miranker, A. and Kaprlus, M., Proteins: Struc. Func. Gen. 11:29-34 (1991); Caflish, A., et al., J. Med. Chem. 36:2142-2167 (1993); Joseph-McCarthy, D., et al., J. Am. Chem. Soc. 123:12758-12769 (2001)]. (These references are incorporated herein by reference in their entirety.)

[0029] What distinguishes the LMC approach from previous fragment-based methods is its ability to compute the actual thermodynamic fragment distributions around the protein, i.e. distributions consistent with thermal fluctuations at physiological temperatures. Information on the thermodynamic distribution is essential for computing free energies of binding, which, as presented further on, is the basic biologically relevant quantity for quantifying the binding affinity of a ligand.

[0030] Indeed, the MCSS approach for example is essentially based on an energy minimization procedure, providing fragment states corresponding to various local minima of the potential energy field representing the fragment-protein interaction. Such a procedure is computationally more expeditious than computing a thermodynamic ensemble of states, but is unable to provide information on entropic effects, essential for free energy estimates.

[0031] For computing the thermodynamic distributions, the LMC code package makes use of a Metropolis Monte Carlo approach [Metropolis, N., et al., J. Chem. Physics 21:1087-1092 (1953)] for sampling from a grand canonical ensemble of states [Adams, D. J., Molecular Physics 29:307-311 (1975); Mezei, M., Molecular Physics 61:565-582 (1987)]. (These references are incorporated herein by reference in their entirety.) In addition to exchanging just energy with a surrounding thermal bath, as in the case of a canonical ensemble, the system described by a grand canonical ensemble exchanges particles (or fragments in the case of LMC) with its surroundings as well. The energy cost associated with inserting/deleting a fragment from the system is determined by its chemical potential. By varying this chemical potential, so-called simulated annealing of the chemical potential, one may vary the average number of fragments in the simulation system. It is shown further on, that measuring the values of the chemical potential at which fragments leave various sites on the protein provides an estimate of the free energy of binding for the different binding modes over the protein surface.

[0032] The practicality of the simulated annealing procedure for estimating binding affinities was demonstrated by Guarnieri and Mezei for differentiating hydration propensities of different DNA grooves [Guarnieri, F. and Mezei, M., J. Am. Chem. Soc. 118:8493-8494 (1996)]. (This reference is incorporated herein by reference in its entirety.) These results were obtained with the Metropolis Monte Carlo (MMC) code developed by the group of Mezei, Mount Sinai School of Medicine, NY. For these simulations, the system was composed of a molecule fraction of DNA surrounded by a varying number of interacting water molecules.

[0033] In its original form, the LMC algorithm carried out a series of calculations similar to the MMC approach for each fragment-type of interest, i.e. simulations in which both the fragment—protein as well as all fragment—fragment interactions were considered. However, it has been acknowledged, that considering fragment-fragment interactions is actually detrimental to the physical interpretation of the simulation results for all fragments but water. Indeed, due to the high dilution of the solute molecules in actual biochemical relevant conditions, considering interactions between non-water fragments is not realistic. Furthermore, the drug leads assembled by LCD usually are composed of only one fragment of each type. Fragment-fragment interactions in the LMC simulation thus lead to undesirable correlation effects.

[0034] Finally, in the original MMC code, carrying out the simulated annealing of the chemical potential for computing the free energies of binding for a single fragment type required the data from multiple ensemble samplings at various B values, i.e. data from multiple simulations. In the absence of fragment-fragment interactions however, the data required for applying the simulated annealing procedure can be directly derived from the sampling of a single simulation. As will be shown further on, this simplification results from the ability of establishing the analytical dependence in B of the fragment density when fragment interactions are omitted. This fact naturally provides an opportunity for significant computational speedup.

[0035] It turns out however that the standard Metropolis Monte Carlo algorithm has difficulty in handling simulations where fragment-fragment interactions are removed. Indeed, the absence of fragment-fragment interactions leads to the possible overlap of fragments and thus to a broad range (typically orders of magnitude) of fragment densities between the higher and lower affinity binding sites on the protein, which the standard Metropolis Monte Carlo scheme has trouble in resolving. This problem has been overcome in the current implementation of LMC by developing a weighted Metropolis Monte Carlo approach.

[0036] The system in which fragment-fragment interactions have been removed can be referred to as being linear by reference to the linear properties of the partial differential equation (Liouville-type) that describes the time-evolution of the fragment density away from thermodynamic equilibrium.

[0037] II. Process

[0038] A. Formulation

[0039] First, the derivation of the single fragment density in the framework of the grand canonical ensemble is presented.

[0040] The potential energy of the system composed of N equivalent, rigid fragments is denoted U(Γ, N). In general, U includes both contributions from fragment-protein and fragment-fragment interactions. The configuration of the system is characterized by

Γ=(Y1,Y2, . . . , YN),  (1)

[0041] where Yi=(xi, Ωi) stands for the position xi and orientation Ωi of the rigid fragment i. In practice, the orientation Ωi is conveniently represented by a unit quaternion q.

[0042] In the grand canonical ensemble, the probability that the system has N fragments in configuration Γ is given by 1$\begin{array}{cc}f\left(\Gamma ,N\right)=\frac{1}{Q}\frac{1}{V^N{\sigma }^N}\frac{1}{N!}\exp \left[\mathrm{BN}-\beta U\left(\Gamma ,N\right)\right],& \left(2\right)\end{array}$

[0043] with the normalization factor given by the grand partition function 2$\begin{array}{cc}Q=\sum _{N=0}^{\infty }\frac{1}{N!}\exp \left(\mathrm{BN}\right)\int \frac{d{Y}^N}{V^N{\sigma }^N}\exp \left[-\beta U\left(\Gamma ,N\right)\right].& \left(3\right)\end{array}$

[0044] Here V is the volume of the system, σ is the volume of orientation space, β=1/KBT, T is the temperature, KB the Boltzmann constant, and B is related to the excess chemical potential μex, i.e. the energy cost in units of β−1 for a particle to leave the system, according to the following relation:

B =βμex+log<N>,  (4)

[0045] where <N> is the average number of fragments in the system. The integral in Eq. (3) is taken over the whole configuration space (Vσ)N.

[0046] Assuming no fragment-fragment interactions, the potential energy U of the system becomes: 3$\begin{array}{cc}U\left(\Gamma ,N\right)=\sum _{i=1}^{N}E\left(Y_i\right),& \left(5\right)\end{array}$

[0047] where E(Y) is the energy of interaction of a single fragment of the considered type with the protein.

[0048] The grand partition function can then be written as 4$\begin{array}{cc}Q=\sum _{N=0}^{\infty }\frac{1}{N!}{\left(\exp \left(B\right)\int \frac{dY}{V\sigma }\exp \left[-\beta E\left(Y\right)\right]\right)}^N=\exp Z,\mathrm{with}& \left(6\right)\\ Z=\exp \left(B\right)\int \frac{dY}{V\sigma }\exp \left[-BE\left(Y\right)\right].& \left(7\right)\end{array}$

[0049] In this case, the probability P(N) for having N fragments in the system is given by 5$\begin{array}{cc}P\left(N\right)=\int d{Y}^{N}f\left(\Gamma ,N\right)=\exp \left(-Z\right)\frac{Z^N}{N!}.& \left(8\right)\end{array}$

[0050] This is simply the Poisson distribution with parameter Z. In particular, the average number of fragments in the system is given by 6$\begin{array}{cc}⟨N⟩=\sum _{N=1}^{\infty }\mathrm{NP}\left(N\right)=Z,& \left(9\right)\end{array}$

[0051] which, according to Eq.(7), thus scales exponentially with B.

[0052] In fact, more generally, the probability P(n, ΔV) of finding n fragments in any given sub-volume ΔV of configuration space is given by a Poisson distribution: 7$\begin{array}{cc}\begin{array}{c}P\left(n,\Delta V\right)\approx \sum _{N=n}^{\infty }\frac{N!}{\left(N-n\right)!n!}{\int }_{\Delta V}d{Y}_1\dots \\ d{Y}_n{\int }_{V}d{Y}_{n+1}\dots d{Y}_{N}f\left(\Gamma ,N\right)\\ =\frac{1}{n!}{\left(\exp \left(B\right){\int }_{\Delta V}\frac{dY}{V\sigma }\exp \left[-\beta E\left(Y\right)\right]\right)}^n\frac{1}{Q}\sum _{N=n}^{\infty }\frac{Z^{N-n}}{\left(N-n\right)!}\\ =z^n/n!,\end{array}\mathrm{with}& \left(10\right)\\ z=\exp \left(B\right){\int }_{\Delta V}\frac{dY}{V\sigma }\exp \left[-\beta E\left(Y\right)\right].& \left(11\right)\end{array}$

[0053] Finally, the single fragment density is given by 8$\begin{array}{cc}\begin{array}{c}{f}_{\mathrm{gc}}\left(Y\right)=\sum _{N=1}^{\infty }N\int d{Y}_2\dots \int d{Y}_{N}f\left(\Gamma =\left(Y,Y_2,\dots ,Y_N\right),N\right)\\ =\exp \left(-Z\right)\frac{1}{V\sigma }\exp \left[B-\beta E\left(Y\right)\right]\sum _{N=1}^{\infty }\frac{1}{\left(N-1\right)!}{Z}^{\left(N-1\right)}\\ =\frac{1}{V\sigma }\exp \left[B-\beta E\left(Y\right)\right],\end{array}& \left(12\right)\end{array}$

[0054] which again scales exponentially with respect to B. Here the subscript “gc” stands for Grand Canonical.

[0055] As expected, note that one recovers Eq. (9) for the average number of fragments in the system by integrating fgc over all configurations:

∫dY fgc(Y)=Z.  (13)

[0056] B. Numerical Method

[0057] Equation (12) for the physical single fragment density shows the large dynamical range that may result from the exponential dependence of this quantity with respect to the single fragment-protein potential energy E(Y). This dependence comes from the possible overlap of the non-interacting fragments. This is not an issue in the presence of fragment-fragment interactions, as an upper bound to the fragment density is set by the tightest possible packing of the molecules.

[0058] The underlying method developed for the WGCMMC approach to enable the accurate resolution of the above-mentioned dynamical range in densities is presented here.

[0059] For numerical purposes, instead of considering a constant B value over the whole system, one may consider a field Bnum(Y) in the single fragment configuration space Y (the subscript “num” standing for numerical). This field can be interpreted as the energy cost for a particle to leave the system specifically from position Y. Instead of Eq. (2), the density of states in this generalized grand canonical ensemble is now given by 9$\begin{array}{cc}{f}_{\mathrm{num}}\left(\Gamma ,N\right)=\frac{1}{Q_{\mathrm{num}}}\frac{1}{V^N{\sigma }^N}\frac{1}{N!}\exp \left[\sum _{i=1}^N{B}_{\mathrm{num}}\left(Y_i\right)-\beta U\left(\Gamma ,N\right)\right],& \left(14\right)\end{array}$

[0060] with the normalization factor (grand partition function) now given by 10$\begin{array}{cc}{Q}_{\mathrm{num}}=\sum _{N=0}^{\infty }\frac{1}{N!}\int \frac{d{Y}^N}{V^N{\sigma }^N}\exp \left[\sum _{i=1}^N{B}_{\mathrm{num}}\left(Y_i\right)-\beta U\left(\Gamma ,N\right)\right].& \left(15\right)\end{array}$

[0061] An analogous derivation as the one used for obtaining Eq. (12) leads to the corresponding single fragment density: 11$\begin{array}{cc}{f}_{\mathrm{gc},\mathrm{num}}\left(Y\right)=\frac{1}{V\sigma }\exp \left[B_{\mathrm{num}}\left(Y\right)-\beta E\left(Y\right)\right].& \left(16\right)\end{array}$

[0062] Thanks to the field Bnum((Y), one now has a direct handle on the value of the density at each position Y of the single particle configuration space. Thus, by appropriately adapting Bnum(Y) during the convergence phase of the Metropolis Monte Carlo simulation, typically through an iterative process, one may obtain good sampling in all regions of interest. For a Bnum field continuous over Y, this would be achieved by taking

B num (Y)≈min(βE(Y)+const, Bmax),  (17)

[0063] leading to similar numerical densities of fragment instances in various regions of space. An upper bound Bmax is set on Bnum to avoid unnecessary sampling in strongly unfavorable positions, i.e., essentially for configurations leading to steric clashes. Such an upper bound ensures to preserve the advantages of the Metropolis Monte Carlo scheme over standard Monte Carlo integration algorithms. In practice, the field Bnum(Y) is typically chosen to be independent of the fragment orientation, and to be piece-wise constant on a 3-D grid in x-space (translational-space). Eq. (16) and (17) also show how the purpose of the Bnum(Y) field could have equivalently been achieved by rescaling the single fragment potential energy field E(Y).

[0064] Making use of the exponential dependence in B of the density, one can infer the physical fragment density fgc(Y) at any B=B0=constant value from the simulation results for a given numerical Bnum(Y) field. Assume that one has a sampling {Γi=(Y1, . . . , YNi)}i=1, . . . , nsnap of nsnap snapshots from the numerical distribution fnum(Γ,N). The average of any single fragment quantity A(Y) over the distribution fgc(Y) is then given by 12$\begin{array}{cc}\begin{array}{c}⟨A⟩=\int dY{f}_{\mathrm{gc}}\left(Y\right)A\left(Y\right)\\ =\int dY{f}_{\mathrm{gc},\mathrm{num}}\left(Y\right)\frac{f_{\mathrm{gc}}\left(Y\right)}{f_{\mathrm{gc},\mathrm{num}}\left(Y\right)}A\left(Y\right)\\ \simeq \frac{1}{n_{\mathrm{snap}}}\sum _{i=1}^{n_{\mathrm{snap}}}\sum _{j=1}^{N_i}{w}_{j}A\left(Y_j\right),\end{array}& \left(18\right)\end{array}$

[0065] where wj is the weight assigned to the fragment state Yj, and defined by 13$\begin{array}{cc}{w}_j=\frac{f_{\mathrm{gc}}\left(Y_j\right)}{f_{\mathrm{gc},\mathrm{num}}\left(Y_j\right)}=\exp \left(B_0-B_{\mathrm{num}}\left(Y_j\right)\right).& \left(19\right)\end{array}$

[0066] Results for any B value can thus be inferred from Eqs. (18)-(19). In this way, thanks to the absence of fragment-fragment interactions, simulated annealing of the chemical potential (i.e. variation of B) can be derived analytically given the sampling data for a single Bnum(Y) field.

[0067] C. Handling WGCMMC data

[0068] The following addresses how the WGCMMC data is to be handled and analyzed.

[0069] The starting point for the data interpretation is the relation linking the WGCMMC data to the association constant Kα characterizing the binding of the considered fragment to a given region on the protein. This relation for Kα is rederived here.

[0070] The association constant Kα characterizes the equilibrium of the binding process

F+P⇄FP,  (20)

[0071] and is defined by 14$\begin{array}{cc}{K}_a=\frac{\left[\mathrm{FP}\right]}{\left[F\right]\left[P\right]},& \left(21\right)\end{array}$

[0072] where [P], [F], and [FP] are respectively the concentrations of protein P alone, fragment F alone, and of a particular protein-fragment complex FP (binding mode). The association constant is the basic biologically relevant quantity.

[0073] Let us consider a single protein in a volume V. For the sake of the following discussion, take V to be large, although for the actual LMC simulation this need not be the case. The protein concentration is thus given by [P]=1/V. Furthermore, let us note n the average number of fragments in the binding volume ΔVb (in general a volume with limits both in translational and orientational space), and N the average total number of fragments in the system, so that [F]=(N−n)/V and [FP]=n/V. The association constant can thus be written 15$\begin{array}{cc}{K}_a=\frac{n/V}{\left(N-n\right)/V1/V}\simeq V\frac{n}{N}& \left(22\right)\end{array}$

[0074] having invoked the thermodynamic limit of large volume V, so that n<<N (N/V→const, for V→∞). The values n and N can be obtained from the fragment density (12): 16$\begin{array}{cc}n={\int }_{\Delta V_b}^{}dY{f}_{\mathrm{gc}}\left(Y\right)=\frac{e^B}{V\sigma }{\int }_{\Delta V_b}^{}dY\exp \left[-\beta E\left(Y\right)\right],& \left(23\right)\\ N={\int }_V^{}dY{f}_{\mathrm{gc}}\left(Y\right)=\frac{e^B}{V\sigma }{\int }_V^{}dY\exp \left[-\beta E\left(Y\right)\right]\simeq e^B,& \left(24\right)\end{array}$

[0075] having again invoked the assumption of high protein dilution, so that the total system volume V is much larger than the effective region of interaction between the fragment and the protein, and thus one may consider E(Y)≅0 in deriving the last approximate equality in (24). The association constant now becomes: 17$\begin{array}{cc}{K}_a=\frac{1}{\sigma }{\int }_{\Delta V_b}^{}dY\exp \left[-\beta E\left(Y\right)\right].& \left(25\right)\end{array}$

[0076] On the basis of Eq. (25) one can also write the association constant in terms of the free energy of binding ΔA:

K α =V exp(−βΔA).  (26)

[0077] where ΔA=AFP−AF, with AFP and AF respectively the free energies of the fragment-protein complex FP and of the fragment F alone: 18$\begin{array}{cc}{A}_{\mathrm{FP}}=-\frac{1}{\beta }\log \left({\int }_{\Delta V_b}^{}dY\exp \left[-\beta E\left(Y\right)\right]\right),& \left(27\right)\\ A_F=-\frac{1}{\beta }\log {\int }_V^{}dY=-\frac{1}{\beta }\log \left(V\sigma \right).& \left(28\right)\end{array}$

[0078] The critical value Bc that is associated to the binding volume ΔVb is defined as the value for which the average number of fragments in the binding site is one. From Eq. (23) follows: 19$\begin{array}{cc}\begin{array}{ccc}n\left(B_c\right)=1& \iff & e^{-B_c}=\frac{1}{V\sigma }{\int }_{\Delta V_b}^{}dY\exp \left[-\beta E\left(Y\right)\right],\end{array}& \left(29\right)\end{array}$

[0079] and from (25), (26) and (29) one sees that Bc is directly related to Ka and ΔA as follows: 20$\begin{array}{cc}{K}_a=V{e}^{-B_c},& \left(30\right)\\ \Delta A=\frac{1}{\beta }{B}_c.& \left(31\right)\end{array}$

[0080] Thus, a low Bc value reflects a high affinity binding mode, and inversely a high Bc value reflects a low affinity mode.

[0081] The critical value Bc can be computed from the WGCMMC data using definition (29), as well as Eqs (18) and (19): 21$\begin{array}{cc}\begin{array}{c}1=n\left(B_c\right)\\ ={\int }_{\Delta V_b}dY{f}_{\mathrm{gc}}\left(Y\right)\cong \frac{1}{n_{\mathrm{snap}}}\sum _{i=1}^{n_{\mathrm{snap}}}\sum _{\mathrm{frag}j\in \Delta V_b}\exp \left[B_c-B_{\mathrm{num}}\left(Y_j\right)\right]\iff \\ B_c=-\log \left(\frac{1}{n_{\mathrm{snap}}}\sum _{i=1}^{n_{\mathrm{snap}}}\sum _{\mathrm{frag}j\in \Delta V_b}\exp \left[-B_{\mathrm{num}}\left(Y_j\right)\right]\right).\end{array}& \left(32\right)\end{array}$

[0082] Equations (30), (31) and (32) provide the basic relations for interpreting the WGCMMC data.

[0083] Binding Analysis

[0084] A first estimate of the binding affinity of a given fragment for different regions on the protein surface can be obtained by assigning a critical Bc to each fragment-residue pair. These Bc values are obtained from the WGCMMC data by applying relation (32), and by assigning a binding volume ΔVb to each residue based of the following proximity criteria: A fragment state is considered to be in proximity of a given residue if at least one fragment-protein atom pair (a, b) is such that

rab<a (RVdW,a+Rvdw,b),  (33)

[0085] where rab is the distance between the two atoms, RVdW is the Van der Walls radii, and α is a numerical parameter (typically α=1.2). The Van der Walls radii are typically defined as half the Lennard-Jones parameter from the considered molecular-mechanics force-field (e.g. AMBER) used for the Monte Carlo simulation.

[0086] The volume defined on the basis of the proximity criteria is in general only a crude estimate of a binding mode volume. The corresponding Bc values must therefore be interpreted accordingly. Nonetheless, comparing sets of Bc values obtained in this way for different fragment types has proven valuable to help identify protein binding sites as follows: A binding site is identified as a set of neighboring residues with low Bc values (high affinity) for multiple fragments with different physico-chemical properties. This approach is based on the assumption that diverse interactions in a localized region are the necessary condition for ensuring the specificity of a binding site. This numerical identification of binding sites is preferably complemented by experimental binding information, such as co-crystal X-ray data and mutational analysis.

[0087] Compared to the above described residue-based proximity criteria, more detailed calculations of the binding mode volumes ΔVb are necessary to provide more accurate estimates of the free energy of binding using Eq. (32). Such improved binding mode volume estimates are determined by identifying “humps” in the fragment distribution. This can be achieved by clustering sampled fragment states belonging to a same potential energy well. For this purpose one makes use of the potential energies saved for the sampled fragment states.

[0088] Chemistry Design

[0089] With the purpose of data reduction, the LCD chemistry design software clumps the sampled fragment instances together. Clumping in LCD is usually carried out at a relatively fine-grained level, so that the clumping volume AVC (limited both in translational and orientational space) is different from a true binding mode volume ΔVb of the fragment. In fact, a binding mode volume is usually composed of many clump volumes. Each clump is thus assigned the Bc value of the binding mode volume to which it belongs.

[0090] Using the WGCMMC-type data, average clump positions xc and quaternion representation qc of average clump orientation can be computed by the following weighted averages: 22$\begin{array}{cc}⟨x_c⟩=\frac{\sum _i{w}_i{x}_i}{\sum _i{w}_i},& \left(37\right)\\ ⟨q_c⟩=\sum _i{w}_i{q}_i->\mathrm{Normalize}{q}_c,& \left(38\right)\end{array}$

[0091] where the sums are over all fragments i in the clump.

[0092] Within the chosen protein binding site, clumps of different fragment types can then be assembled into actual candidate drug leads, usually composed of four to five fragments. Assembly of fragments is carried out based on binding affinity of the different fragments (Bc values), and on geometric proximity using a variety of rules by which organic fragments may bond together, as is well known in the art.

[0093] D. Process Implementation

[0094] In light of the above analytical description of WGCMMC processing, its logic can be implemented in the broader simulation context as illustrated in FIG. 1, according to an embodiment of the invention. The overall process starts at step 110. In step 120, a model is constructed for the molecules to be simulated, i.e., a protein as well as different types of rigid molecular fragments whose interaction with the protein will be analyzed. In step 130, the thermodynamic equilibrium of the system is modeled so that the interactions between a given fragment type and the protein at thermodynamic equilibrium can be understood. This step results in simulation data that includes, for each fragment state, the fragment's position, orientation, weight, and fragment-protein energy. Step 130 is carried out for each fragment type of interest. These simulations for the different fragment types are performed independently from one another and can therefore be conveniently carried out simultaneously on a cluster of computers, without any communication required between processes. In step 140, potential binding sites are identified on the protein. In step 150, fragments are assembled into drug leads. The overall process concludes at step 160. Each of these steps is described in greater detail below.

[0095] Molecule Preparation

[0096] Step 120, the preparation of the molecular model, is illustrated in FIG. 2. This process starts at step 210. Protein preparation takes place in step 220. A protein can be viewed as a biological macro-molecule to which a prospective ligand binds. The basic protein structure is provided by experimental X-ray crystallography data, typically downloaded from a data base [e.g. the from the Protein Data Bank (PDB), Research Collaboratory for Structural Bioinformatics (RCSB), Rutgers Univ., NJ]. If required, the protein structure is completed for missing substructures, which in some cases may be a limited number of heavy atoms or, in other cases, entire segments of an amino-acid chain. Hydrogen atoms, not resolved by X-ray crystallography, are added as well. Conformer and protonation state issues for the amino-acids HIS, ASP, GLU, CYS, TYR, THR, and SER are also resolved at this stage. Such a process for protein preparation is disclosed and claimed in a co-pending U.S. patent application Ser. No. 60/450,711, filed on Mar. 3, 2003, and incorporated herein by reference in its entirety.

[0097] Fragment preparation takes place in step 230. The structure and partial charges of the small organic fragments are completed with an ab initio, i.e. quantum mechanical based, code. This calculation is typically carried out in the framework of the Density Functional Theory (DFT) approximation using the code Gaussian (M. J. Fish et.al., “Gaussian 98, revision A.9,” 1998. Gaussian Inc., Pittsburgh, Pa.). This step also assigns the atom types from the molecular mechanics force-field (e.g. AMBER) applied in the subsequent Monte Carlo simulation. The process concludes at step 240.

[0098] The step of modeling the thermodynamic system of the protein-fragment interaction is illustrated in greater detail in FIG. 3, according to an embodiment of the invention. The process starts at step 310. In step 320, a convergence phase of the weighted grand canonical Metropolis Monte Carlo simulation is executed. This is followed by a sampling phase in step 330. Steps 320 and 330 are described in greater detail below. The resulting simulation data is saved in step 340. The process concludes in step 350.

[0099] Convergence Phase of LMC Simulation

[0100] Step 320, the convergence phase of the LMC simulation, is illustrated in FIG. 4. In this first stage of the LMC simulation, the numerical B-field, Bnum, and the Markov chain generated by the LMC stepping are converged.

[0101] The process starts with step 410. In step 420, the simulation space is subdivided with a grid. Typically, the 3-dimensional translational space of the simulation system is subdivided by an orthogonal, equidistant grid, with centers xi. Grid size is based on the variation scale of the interaction force-field, typically of the order of one Angstrom. The Bnum field is then initialized on this grid to a constant value Bo (Bnum≡Bo), as shown in step 430.

[0102] Stepping of the system state is then carried out using the Metropolis Monte Carlo scheme for grand canonical simulations [Adams, D. J., Molecular Physics 29:307-311 (1975); Mezei, M., Molecular Physics 61:565-582 (1987)]. (These references are incorporated herein by reference in their entireties.) At regular intervals in the stepping of the convergence phase, sufficiently long to ensure decorrelation of states, the fragment distributions are saved, as shown in 440.

[0103] With progressively improved statistics, the fragment distributions are then monitored periodically as given in step 450, where the weighted number of sampled fragments in each grid cell xi is computed as follows: 23$\begin{array}{cc}{n}_{B=0}\left(x_i\right)=\frac{1}{n_{\mathrm{samples}}}\sum _{\mathrm{samples}}\sum _{\mathrm{frag}jin\mathrm{cell}i}\exp \left[-B_{\mathrm{num}}\left(Y_j\right)\right],& \left(40\right)\end{array}$

[0104] where nsamples is the number of samples saved up to the current point in the convergence phase. Equation (40) is an application of Eqs. (18)-(19) for B0=0.

[0105] Based on these statistics, the field Bnum(x) is then adapted in step 460, by making use of the exponential dependence in Bnum(x) of the number of fragments in each grid cell i. In this way, each cell is assigned a constant value Bnum(xi) as follows: 24$\begin{array}{cc}{B}_{\mathrm{num}}\left(x_i\right)=\log \left(\frac{n_{\mathrm{target}}}{n_{B=0}\left(x_i\right)}\right),& \left(41\right)\end{array}$

[0106] the goal being to achieve a similar average number of sampled fragments ntarget within all cells. An upper bound Bmax is set on Bnum to avoid spending too much computing time on sampling very unfavorable positions, i.e., mainly for configurations leading to steric clashes or for fragment states far away from the protein surface where the binding interaction is low. In this way one still ensures the numerical advantages of the Metropolis Monte Carlo scheme over basic Monte Carlo integration algorithms.

[0107] Adapting the field Bnum(X) is an iterative process of steps 440 to 460. Indeed, the first Bnum updates are based on some very non-uniform sampling, thorough in deep energy pockets, but poor in shallow ones. As the Bnum(x) field is adapted, the sampling is globally improved and the adjustment of Bnum(x) can be further refined.

[0108] In step 470 of the convergence phase, the Bnum(x) field is finally kept fixed, which enables the Markov chain to fully equilibrate.

[0109] The acceptance probabilities for the various types of Monte Carlo steps in the framework of the grand canonical ensemble with spatially varying Bnum(x) field are as follows:—

[0110] Moving a fragment within the simulation system: Assuming symmetric attempts, moving a fragment from position Ya=(xa, Ωa) to position Yb=(xb, Ωb) is accepted with probability:

acc(Ya→Yb)=min(1, α),  (42)

with α=exp([B(xb)−B(xa)]−β[E(Yb)−E(Ya)]).  (43)

[0111] Inserting a fragment into the simulation system: Assuming no biased sampling, and considering that N fragments are already present in the system, the probability of accepting the insertion of a fragment at position Y=(x, Ω) is given by: 25$\begin{array}{cc}\mathrm{acc}\left(N->N+1\right)=\min \left(1,\alpha \right),\mathrm{with}& \left(44\right)\\ \alpha =\frac{1}{N+1}\exp \left(B\left(x\right)-\beta E\left(Y\right)\right).& \left(45\right)\end{array}$

[0112] Deleting a fragment from the simulation system: The probability of deleting a fragment at position Y=(x, Ω), assuming that N+1 fragments are initially in the system, is given by:

acc(N+1→N)=min(1, α),  (46)

with α=(N+1)exp(−B(x)+βE(Y)).  (47)

[0113] Equations (42) to (47) can naturally be generalized to various types of biased sampling, such as preferential sampling or cavity bias.

[0114] Sampling Phase of MC Simulation

[0115] The numerical B-field, Bnum, is kept fixed throughout the second stage, the so-called sampling phase, of the MC simulation. This phase, step 330 of FIG. 3, is illustrated in greater detail in FIG. 5. The process starts with step 510. In step 520, Bnum(x) is further kept fixed. In step 530, the equilibrated Markov chain is sampled periodically at sufficiently decorrelated states until a statistically appropriate amount of sampling data is acquired. As shown in step 540, saving the state of the system consists of storing the positions x, orientations Ω, weights w=exp(−Bnum(x)), and fragment-protein potential energies E(Y) of all fragments currently present in the system. The sampling process concludes at step 550.

[0116] Identifying Binding Modes

[0117] FIG. 6 illustrates the process of identifying potential binding sites, according to an embodiment of the invention. The process starts with step 610. In step 620, logic such as the Locus Binding Analysis (LBA) software package begins execution. In step 630, a value Bc is assigned to each fragment-residue pair. In step 640, potential binding sites are identified on the basis of the Bc values. As discussed above, these Bc values are obtained from the WGCMMC data by applying relation (32), where the volume ΔVb is defined for each residue on the basis of the proximity criteria. Recall from Eq. (33) above that a fragment is considered to be in proximity of a given residue if at least one fragment-protein atom pair (a, b) is such that

rab<α(RVdW,a+RVdW,b),  (33)

[0118] where rab is the distance between the two atoms, RVdW is the Van der Walls radii, and α is a numerical parameter (typically α=1.2). The Van der Walls radii are typically defined as half the Lennard-Jones parameter from the considered molecular-mechanics force-field (e.g. AMBER) used for the Monte Carlo simulation. A binding site is then identified as a set of residues with low Bc values (high affinity) for multiple fragment types with diverse physico-chemical properties. The process concludes at step 650.

[0119] Assembling Fragments in the Binding Site

[0120] Step 150 of FIG. 1, the step of assembling fragments into drug leads, is illustrated in greater detail in FIG. 7, according to an embodiment of the invention. The process starts with step 710. With the purpose of data reduction, fragment instances are clumped together in step 720. Clumping is carried out at a relatively fine-grained level (both in translational and orientational space), so that the clumping volume ΔVc is different from a true binding volume. In fact, a binding mode volume is usually composed of many clump volumes. The purpose of this clumping is to achieve some level of data reduction before carrying on with the fragment assembly into drug leads. From a combinatorial point of view, this assembly indeed becomes increasingly complex and therefore computationally intensive with increasing number of considered fragment poses.

[0121] In steps 730 through 750, weighted average clump positions xc, and quaternion representation of weighted average clump orientation qc is computed as described earlier: 26$\begin{array}{cc}⟨x_c⟩=\frac{\sum _i{w}_i{x}_i}{\sum _i{w}_i},& \left(37\right)\\ ⟨q_c⟩=\sum _i{w}_i{q}_i->\mathrm{Normalize}{q}_c,& \left(38\right)\end{array}$

[0122] where the sums are over all fragments i in the clump.

[0123] In the same way, as appears in step 760, one may also compute the average potential energy of the clump: 27$\begin{array}{cc}⟨E_c⟩=\frac{\sum _i{w}_i{E}_i}{\sum _i{w}_i},& \left(39\right)\end{array}$

[0124] where Ei is the potential energy of interaction of fragment i with the protein.

[0125] In step 770, each clump is assigned the Bc value of the binding mode volume to which it belongs.

[0126] In step 780, within the chosen protein binding site, clumps of different fragment types are then assembled into actual candidate drug leads, usually (though not always) composed of four to five fragments. Assembly of fragments is carried out based on binding affinity of the different fragments (Bc values), and on geometric proximity, using a variety of rules by which organic fragments may bond together as is well known in the art.

[0127] III. Computing Environment

[0128] The present invention may be implemented using software and may be implemented in conjunction with a computing system or other processing system. An example of such a computer system 800 is shown in FIG. 8. The computer system 800 includes one or more processors, such as processor 804. It is to be noted that the here-described fragment-based computation is particularly well suited for being carried out on a computer cluster, each cluster node computing the interaction of a given fragment type with the target protein. The processor 804 is connected to a communication infrastructure 806, such as a bus or network. Various software implementations are described in terms of this exemplary computer system. After reading this description, it will become apparent to a person skilled in the relevant art how to implement the invention using other computer systems and/or computer architectures.

[0129] Computer system 800 also includes a main memory 808, preferably random access memory (RAM), and may also include a secondary memory 810. The secondary memory 810 may include, for example, a hard disk drive 812 and/or a removable storage drive 814, representing a magnetic tape drive, an optical disk drive, etc. The removable storage drive 814 reads from and/or writes to a removable storage unit 818 in a well-known manner. Removable storage unit 818 represents a magnetic tape, optical disk, or other storage medium that is read by and written to by removable storage drive 814. As will be appreciated, the removable storage unit 818 can include a computer usable storage medium having stored therein computer software and/or data.

[0130] In alternative implementations, secondary memory 810 may include other means for allowing computer programs or other instructions to be loaded into computer system 800. Such means may include, for example, a removable storage unit 822 and an interface 820. An example of such means may include a removable memory chip (such as an EPROM, or PROM) and associated socket, or other removable storage units 822 and interfaces 820 which allow software and data to be transferred from the removable storage unit 822 to computer system 800.

[0131] Computer system 800 may also include one or more communications interfaces, such as network interface 824. Network interface 824 allows software and data to be transferred between computer system 800 and external devices. Examples of network interface 824 may include a modem, a network interface (such as an Ethernet card), a communications port, a PCMCIA slot and card, etc. Software and data transferred via network interface 824 are in the form of signals 828 which may be electronic, electromagnetic, optical or other signals capable of being received by network interface 824. These signals 828 are provided to network interface 824 via a communications path (i.e., channel) 826. This channel 826 carries signals 828 and may be implemented using wire or cable, fiber optics, an RF link and other communications channels.

[0132] In this document, the terms “computer program medium” and “computer usable medium” are used to generally refer to media such as removable storage units 818 and 822, a hard disk installed in hard disk drive 812, and signals 828. These computer program products are means for providing software to computer system 800.

[0133] Computer programs (also called computer control logic) are stored in main memory 808 and/or secondary memory 810. Computer programs may also be received via communications interface 824. Such computer programs, when executed, enable the computer system 800 to implement the present invention as discussed herein. In particular, the computer programs, when executed, enable the processor 804 to implement the present invention. Accordingly, such computer programs represent controllers of the computer system 800. Where the invention is implemented using software, the software may be stored in a computer program product and loaded into computer system 800 using removable storage drive 814, hard drive 812 or communications interface 824.

IV. CONCLUSION

[0134] While various embodiments of the present invention have been described above, it should be understood that they have been presented by way of example, and not limitation. It will be apparent to persons skilled in the relevant art that various changes in detail can be made therein without departing from the spirit and scope of the invention. Thus the present invention should not be limited by any of the above-described exemplary embodiments.

Claims

1. A method for modeling a system that includes a protein and a plurality of fragment types in order to identify drug leads, the method comprising: initiating a weighted grand canonical Metropolis Monte Carlo simulation of the system; subdividing the space of the simulation system with a grid, with xi the centers of the grid cells; initializing a numerical chemical potential field Bnum=B0 on the grid; periodically sampling the Markov chain associated with the Metropolis Monte Carlo simulation, so as to compute the weighted number of sampled fragments per cell: 28nB=0⁡(xi)=1nsamples⁢∑samples⁢∑frag⁢ ⁢j⁢ ⁢i⁢ ⁢n⁢ ⁢cell⁢ ⁢i⁢exp⁡[-Bnum⁡(Yj)];iteratively adapting the field Bnum(x) according to 29Bnum⁡(xi)=log⁡(ntargetnB=0),fixing the field Bnum(X) such that the Markov chain associated with the Metropolis Monte Carlo simulation equilibrates; and outputting samples from the equilibrated Markov chain.

2. The method of claim 1, further comprising: sampling the Markov chain periodically, with sufficiently long interspacing to ensure decorrelated states of the system; and saving positions, orientations, fragment-protein potential energies, and statistical weights for all fragments present in a current state of the system.

3. The method of claim 2, further comprising: performing binding analysis of the system, based on the positions, orientations, fragment-protein potential energies, and statistical weights for all fragment states provided by the sampling.

4. The method of claim 3, wherein said performing step comprises: i) making use of the properties of the grand canonical ensemble to estimate the binding affinity of the fragment for different regions of the protein surface by assigning a critical value Bc to each fragment-residue pair, using the positions, orientations, and statistical weights for all fragment states provided by the sampling; and ii) identifying potential binding sites on the protein based on the Bc values.

5. The method of claim 2, further comprising: assembling the fragments into drug leads for a considered binding site, based on binding affinity of the fragment types (Bc values) for the considered binding site, and on geometric proximity using rules by which organic fragments may bond together.

6. A computer program product comprising a computer usable medium having computer readable program code that enables a computer to model a system that comprises a protein and a plurality of fragments in order to identify drug leads, the computer program product comprising: first computer readable program code that initiates a weighted grand canonical Metropolis Monte Carlo simulation; second computer readable program code that causes the computer to subdivide the space of the simulation system with a grid, with xi the centers of the grid cells; third computer readable program code that causes the computer to initialize a field Bnum(xi)=B0; fourth computer readable program code that causes the computer to compute the weighted number of sampled fragments per cell, 30nB=0⁡(xi)=1nsamples⁢∑samples⁢∑frag⁢ ⁢j⁢ ⁢i⁢ ⁢n⁢ ⁢cell⁢ ⁢i⁢exp⁡[-Bnum⁡(Yj)],fifth computer readable program code that causes the computer to iteratively adapt the field Bnum(x) according to 31Bnum⁡(xi)=log⁡(ntargetnB=0),sixth computer readable program code that causes the computer to keep the field Bnum(x) fixed, so that the Markov chain associated with the Metropolis Monte Carlo scheme can equilibrate; and seventh computer readable program code that causes the computer to output samples from the equilibrated Markov chain.

7. The computer program product of claim 6, further comprising: seventh computer readable program code that causes the computer to sample the Markhov chain periodically at sufficiently decorrelated states of the system; and eighth computer readable program code that causes the computer to obtain positions, orientations, fragment-protein potential energies, and statistical weights for all fragments present in a current state of the system.

8. The computer program product of claim 7, further comprising: ninth computer readable program code that causes the computer to perform binding analysis based on the positions, orientations, and statistical weights for all fragments at each sampled state of the system.

9. The computer program product of claim 8, wherein said ninth computer readable program code comprises: computer readable program code that causes the computer to assign a critical value Bc to each fragment-residue pair based on the positions, orientations, and statistical weights for all fragments at each state; and computer readable program code that causes the computer to identify potential binding sites on the protein based on the Bc values.

10. The computer program product of claim 8, further comprising: tenth computer readable program code that causes the computer to assemble the fragments into drug leads for a considered binding site based on binding affinity of the fragment types (Bc values), and on geometric proximity using rules by which organic fragments may bond together.

11. A system for modeling a system that includes a protein and a plurality of different fragment types in order to identify drug leads, the system comprising: A. means for initiating a weighted grand canonical Metropolis Monte Carlo simulation of the system; B. means for subdividing the space of the simulation system with a grid, with xi the centers of the grid cells; C. means for initializing a numerical chemical potential field Bnum=B0 on the grid; D. means for computing the weighted number of sampled fragments per cell, 32nB=0⁡(xi)=1nsamples⁢∑samples⁢∑frag⁢ ⁢j⁢ ⁢in⁢ ⁢cell⁢ ⁢i⁢exp⁡[-Bnum⁡(Yj)];E. means for iteratively adapting the field Bnum(X) such that 33Bnum⁡(xi)=log⁡(ntargetnB=0),F. means for fixing the field Bnum(x) such that the associated Markhov chain equilibrates; and G. means for outputting samples from an equilibrated Markov chain.