SCORING FUNCTION FOR MOLECULAR ALIGNMENT When creating a hypothesis for an alignment of molecules, it is desirable for atoms (from different molecules that are being aligned) with similar properties to be near to each other. Intermolecular atom pairs that are far away should have little inﬂuence on the goodness of an alignment. This can be achieved by associating each intermolecular atom pair with similar properties with a Gaussian scoring term: w exp(-Rr2), where r is the distance between the two atoms, R is a parameter controlling the range of the interaction, and w is the weight of this term. It is also desirable for the alignment hypothesis for atoms with incompatible properties (e.g., hydrophobic vs hydrophilic) to not superpose. This can be achieved using the same functional form but with a negative value for w. Besides atoms, hydrogen-bond pro- jected features can also be considered. These features represent the expected positions of the hydrogen-bonding partner on the binding pocket based on the ligand geometry. The overlap of hydrogen-bond projected features of the same type should be favorable for an alignment hypothesis. Again, the same Gaussian score can work, with r being the distance between the intermolecular features. While the Gaussian function is widely used in molecular alignment methods (e.g., FAP,4 Pharao,5 MOE,16 SQ,19 and SEAL20), there are other functional forms that can provide a similar effect. These include piecewise linear functions (e.g., FLAME),12 and piecewise linear and quadratic func- tions (e.g., TFIT).21 It can be construed that indeed any bell- shaped function can work in a similar fashion. In their pioneering work 20 years ago, Kearsley and Smith20 noticed that the Lorentzian function (1/(1 + Rr2)) yielded similar results to the Gaussian but was signiﬁcantly faster to compute. However, the advance in computational technology to calculate exponents has eliminated this advantage of the Lorentzian function. Simple tests showed that the calculation speeds are now similar. For this work we chose the Gaussian because it involves a minimal number of tunable parameters and has a simple, continuous shape. Another consideration is that function optimization procedures often require the ﬁrst and second derivatives of the function. In the case of the Gaussian, once the function has been calculated, its ﬁrst and second derivatives can be obtained without much extra work. A scoring function S could thus be constructed based on a sum of Gaussian terms: S ) ∑ Tk ) k ∑ wk ∑ ∑ exp(-Rri2j) k i∈Ak j∈Bk In the formula, Tk is a term encouraging/discouraging the overlap of a particular pair of types (Ak, Bk) of atoms/ features; wk is the weight of the term; i and j are atoms/ features from the two molecules; rij is the distance between the intermolecular atom/feature pair; and R is a parameter controlling the range of this type of interaction. Details of the various terms, Tk’s, are given in Table 1. Ligands binding in the same mode are expected to have substantial volume overlap. The ﬁrst term in the scoring function, T0, rewards overlap between the heavy atoms of the different ligands. Next, hydrogen bonding is an important interaction deter- mining ligand binding. The terms T1 and T2 encourage the overlap of pairs of intermolecular donors and acceptors. Following this is the hydrophobic term T3 rewarding the overlap of hydrophobic atoms. After these “attractive” terms J. Chem. Inf. Model., Vol. 50, No. 9, 2010 1725 Table 1. Details of the Scoring Function Terms, Tk’sa atom/feature type Tk Ak (molecule 1) Bk (molecule 2) wk R T0 heavy heavy 1 Ra T1 donor donor wDA Ra T2 acceptor acceptor wDA Ra T3 hydrophobic hydrophobic wHH Ra T4 hydrophobic acceptor/donor wR1 Ra T5 acceptor/donor hydrophobic wR1 Ra T6 acceptor-not-donor donor-not-acceptor wR2 Ra T7 donor-not-acceptor acceptor-not-donor wR2 Ra T8 donor projected feature donor projected feature wPF Rp T9 acceptor projected feature acceptor projected feature wPF Rp a The scoring function consists of a summation of terms, Tk’s, each of which has the form: Tk ) wk ∑ ∑ exp(-Rr2ij) i∈Ak j∈Bk Each row of the table corresponds to one term, Tk, and i and j denote atoms (for T0-T7) or features (for T8 and T9) coming from the ﬁrst and second molecule, respectively. These atoms/features have to be of certain types, Ak and Bk, as given in the second and third columns of the table, and rij is the distance between the atoms/features. The names of the various constants for the different terms, wk and R, are given in the fourth and ﬁfth columns of the table. All the atom-based terms have the same value for R, namely Ra, while all hydrogen-bond projected feature-based terms have another value for R, namely Rp. come the “repulsive” terms that discourage the overlap of intermolecular atom pairs of incompatible types. First, hydrophobic and hydrophilic atoms (atoms capable of hydrogen bonding) are not expected to occupy the same part of the binding pocket, as reﬂected by T4 and T5. Second, hydrogen-bond acceptors that are not also donors are expected to avoid positions of donors that are not also acceptors. This is reﬂected by terms T6 and T7. The last two terms, T8 and T9, are based on hydrogen-bond projected features. They encourage the overlap of hydrogen-bond projected features of the same type. To keep the number of parameters low, the same value of R, Ra, was used for all atom-based terms. The projected feature terms used a different value of R, Rp, allowing these terms to be more diffuse than the atom-based terms if necessary. Note that we have ordered the terms in the scoring function. Atom-based terms were followed by projected feature terms. For the atom-based terms, attractive terms were followed by repulsive terms. And within the attractive terms, the ordering was by the expected importance and intensity of that type of interaction. Based on this order, the scoring function was built up term by term. At each stage, the one or two parameter(s) of the newly added term was/were optimized. A series of scoring functions were created using different values for the newly added parameter(s). The goodness of a scoring function was measured, as described below. The parameter value(s) corresponding to the best scoring function in the series would be picked. Then we proceeded to the next stage to optimize the next parameter(s). Test Sets. After the scoring function was optimized, it was incorporated into the ﬂexible alignment functionality in molecular operating environment (MOE).22 Molecular align- ment results were compared with two recent publications.6,23 Proschak et al.6 performed rigid-body alignments for eight thermolysin ligands using their Shapelets method and