MolAlign: an algorithm for aligning multiple small molecules

Abstract

In small molecule drug discovery projects, the receptor structure is not always available. In such cases it is enormously useful to be able to align known ligands in the way they bind in the receptor. Here we shall present an algorithm for the alignment of multiple small molecule ligands. This algorithm takes pre-generated conformers as input, and proposes aligned assemblies of the ligands. The algorithm consists of two stages: the first stage is to perform alignments for each pair of ligands, the second stage makes use of the results from the first stage to build up multiple ligand alignment assemblies using a novel iterative procedure. The scoring functions are improved versions of the one mentioned in our previous work. We have compared our results with some recent publications. While an exact comparison is impossible, it is clear that our algorithm is fast and produces very competitive results.

Appendix

Termination criterion for simultaneous multiple pose generation

After each cycle (each query triplet), the pair pose registers may be updated. The raw score improvement is defined to be the average (over all registers) of the improvement of the score at the fifth slot (i.e., the fifth best score in the register) over this cycle. The first slot is not used because it does not get updated very often (for the first slot to be updated, a new pose with a score better than any of the existing ones must be found). Position 5 is used because each register is guaranteed to have more than five slots, as per Eq. (2). In principle, we would terminate the iteration when the score improvement, SI, becomes too small. We noticed that some cycles produced very small SIs although there was still significant room for improvement. In order to smooth out the fluctuation, we average out the SI of the current cycle with the SIs of a few past cycles to create an effective SI, and use this to determine termination. Let ${n_{curr}}$ be the number of compatible triplets (which equals the number of conformer poses in the grand assembly) encountered in the current cycle. Let ${n_{avg}}$ be the hitherto average number of compatible triplets per cycle (i.e., ${n_{avg}}$ is the sum of all ${n_{curr}}$′s encountered so far, divided by the number of cycles). The effective SI of the current cycle, ${s_{curr}}$, is given by

$${s_{curr}}=~\frac{{{n_{avg~}}{s_{prev}}~+~{n_{curr~}}{s_0}}}{{{n_{avg}}~+~{n_{curr}}}}$$

(3)

where ${s_{prev}}$ is the effective SI of the previous cycle, and ${s_0}$ is the raw score improvement, as described above. In other words, ${s_{curr}}$ is a weighted average of ${s_{prev}}$ and ${s_0}$. ${s_0}$ is weighted less if there are few compatible triplets in the current cycle. This should make the effective SI, ${s_{curr}}$, less prone to statistical fluctuations than the raw score improvement, $~{s_0}$.

Iteration terminates when the effective SI drops below a cut-off value, ${s_{cut}}$. ${s_{cut}}$ is set to a value of 0.01. To put this in perspective, the pairwise score (according to any scoring scheme) of one hydrophobic atom overlapping on another is 2.0.

For the first cycle, the value of ${s_{prev}}$ in Eq. (3) is set to ${2^5} \times {s_{cut}}$. 5 is the index of the slot used for determining the raw score improvement, ${s_0}$, as mentioned above. This way, if ${s_0}$ is 0 and ${n_{curr}}$ remains unchanged over the cycles (implying that ${n_{avg}}$ always equals ${n_{curr}}$), the iteration would go on for 5 cycles before terminating. Put in another way, if only one slot of each pair pose register is filled per iteration cycle, there is time for the register to be filled to the fifth slot, where ${s_0}$ is determined.

Choosing random query triplets for simultaneous multiple pose generation

Early in our research, we simply chose a random triplet from a random conformer of a random ligand. We found this inefficient. Since we would like each pair pose register to be at least partially filled, we can bias the choice of the ligands to encourage ligands with empty registers to be chosen. First, ligands having fewer than three representative points are not involved, because they do not have representative triplets. An array is constructed out of the indices of all ligands having at least three representative points. Each index appears once. We then go through all ligand pairs. For each ligand pair whose pair pose register is still empty, three copies of the index of each ligand in the pair are inserted into the array. (Experiments showed that using five instead of three has little effect on the results.) Finally, a random element from this array is selected to give a random ligand. Then a random conformer of this ligand is chosen, and a random triplet of representative points of this conformer is chosen. Similar to the systematic pose generation method, triplets that are geometrically too small are never chosen. The diameter (see “Systematic pose generation” section) of a chosen triplet has to be at least half that of the smallest diameter amongst the representative point sets of all conformers of all ligands.

Starting points for multiple ligand alignment

As described in the section “Generating multiple ligand assemblies”, one base template conformer gives rise to one starting point assembly, and subsequently one solution assembly. We discovered that using a maximum of 250 starting points strikes a balance between good results and speed. Before the iteration begins, each conformer of each ligand is used as a base template to generate a starting assembly. Given one base template conformer, for each of the other ligands, the best scoring pairwise alignment pose between the ligand and the base template is used to build the starting assembly. It is possible that no pairwise alignment pose exists in the pair pose register between the base template conformer and the ligand. If this is the case, the missing ligand count of this base template will be increased by one. After all the non-missing ligands have been positioned in the assembly, each missing ligand will be positioned using the non-missing ligands as helpers, one at a time. For each missing ligand, amongst all poses generated using the helpers, the one resulting in the best score for the assembly (all non-missing ligands plus the current ligand) is used. Any ligands still missing after this step will be positioned randomly. The multiple ligand alignment score of the starting point assembly will then be calculated, using fast scoring. The full list of starting point assemblies will then be sorted, first by the missing ligand count, then by the multiple ligand alignment score. If there are more than 250 assemblies, the 250 with the smallest missing ligand counts and the best scores are used as starting points. In a real run, there is of course no need to patch the missing ligands if there are already enough starting points without any missing ligands.