Crossover operators for molecular graphs with an application to virtual drug screening

Abstract

Genetic algorithms are a powerful method to solve optimization problems with complex cost functions over vast search spaces that rely in particular on recombining parts of previous solutions. Crossover operators play a crucial role in this context. Here, we describe a large class of these operators designed for searching over spaces of graphs. These operators are based on introducing small cuts into graphs and rejoining the resulting induced subgraphs of two parents. This form of cut-and-join crossover can be restricted in a consistent way to preserve local properties such as vertex-degrees (valency), or bond-orders, as well as global properties such as graph-theoretic planarity. In contrast to crossover on strings, cut-and-join crossover on graphs is powerful enough to ergodically explore chemical space even in the absence of mutation operators. Extensive benchmarking shows that the offspring of molecular graphs are again plausible molecules with high probability, while at the same time crossover drastically increases the diversity compared to initial molecule libraries. Moreover, desirable properties such as favorable indices of synthesizability are preserved with sufficient frequency that candidate offsprings can be filtered efficiently for such properties. As an application we utilized the cut-and-join crossover in REvoLd, a GA-based system for computer-aided drug design. In optimization runs searching for ligands binding to four different target proteins we consistently found candidate molecules with binding constants exceeding the best known binders as well as candidates found in make-on-demand libraries.Scientific contributionWe define cut-and-join crossover operators on a variety of graph classes including molecular graphs. This constitutes a mathematically simple and well-characterized approach to recombination of molecules that performed very well in real-life CADD tasks.

Keywords: Crossover; Genetic algorithm; Graph theory; Virtual screening.

Fig. 1

An example of a 2-cut on naphthalene. The application of the cut set $C(V^{\prime },V^{\prime \prime })$ leads to the formation of the shown fragments, i.e the induced subgraphs $G[V^{\prime }]$ and $G[V^{\prime \prime }]$

Fig. 2

Lewis structure (left) and molecule multigraph (right) of allyl alcohol

Fig. 3

Increasingly restrictive join operators. While a join allows any nodes of two subgraphs to be linked, for a cut node preserving join only connections between nodes involved in the cut sets of their respective parental graphs are permitted. A degree preserving join further requires all nodes of the offspring graph to adhere to the degree of their respective parental vertices. The most restrictive join presented is the natural join, which introduces edge weights, e.g. bond orders, that must be conserved

Fig. 4

Choosing $V^{\prime }$ and $W^{\prime }$ as the top vertices of the two triangles yields two cuts consisting of two edges, each with weight 1. The only possible join in this case is a multigraph. This join is degree-preserving, as each of the two vertices still has two incident edges, but is not natural since it converted (parts of) a simple graph into a multigraph

Fig. 5

Cut-and-join crossover of two planar graphs can be restricted to yield again a planar graph. Since every minimal cut of a planar graph is cycle in its complement, there is a circular order to the severed edges and thus of the vertices incident with the cutset. It sufficies to preserve the circular order of these vertices (in arbitrary orientation). The example here is in addition degree-preserving

Fig. 6

Left: Search-tree of the cut enumeration algorithm. Right: The different steps of the cut enumeration algorithm shown on camphor. Enumeration of a subtree can be halted, if a cut results in > 2 connected components or the cut node set has already been visited (colored red)

Fig. 7

Left: Functional groups can be interconverted by a simple cut-and-join crossover. The general position (primary, secondary, tertary, quaternary) of a functional group can be changed by means of crossover operations that replace a hydrogen atom by a larger fragment (dotted arrows). Right: Cut-and-join crossover involving hydrogens enable access to simple reactions

Fig. 8

Within a minimal polycyclic set, bridged polycycles sharing multiple edges are accessible as well as cubanes

Fig. 9

A minimal set of base graphs $\mathcal{B}$ needed to access all graphs of degree at most 4 using a sequence of natural join operators

Fig. 10

Diversity of crossover products taken from unrestricted cut-and-join operation on a sample of 10 molecules from the USPTO-10k. The homogeneous distribution of the recombinants in the t-SNE embedding shows that offspring molecules cover a high molecular diversity

Fig. 11

Synthesizability scores calculated by SAScore and SCScore on the 10-sample molecule recombinant set in comparison to the USPTO-10k dataset. Both tools show an increase in the difficulty of synthesizablility. The average $\mathrm{\Delta }$SA Score is listed in Additional file 4

Fig. 12

Frequencies of molecules with a better SAscore than the nth percentile of USPTO-10k molecules are depicted. A significant fraction of the recombinant offsprings has good SAscores, and thus it is feasible in each step of the GA to reject candidate offsprings that are difficult to produce according to their SA score

Fig. 13

Molecular weight (MW) distribution of the 10 molecule sample set. The cut-and-join crossover operator shows an increase of molecular weight. The $\mathrm{\Delta }$MW average of the different sets compared to the USPTO-10k can be found in Additional file 4. For better readability 56 molecules with a MW >1100 in the USPTO-10k were removed

Fig. 14

Distribution of ring topologies of parent and recombinant molecules. Molecules are assigned a category based on the priority: heteroaromatic > aromatic > heterocyclic > carbocyclic > acyclic. The cut-and-join crossover operator tends to generate more polycycles, especially with heteroaromatic and aromatic rings. This is also due to the sample, which is also largely composed of heteroaromatics

Fig. 15

Distribution of normalized docking scores for four protein targets using either only known actives, random samples from the Enamine REAL space, REvoLd exploring the Enamine REAL space, or our custom chemical space. Dotted lines are docking scores for the co-crystallized ligands from the used protein structure. Next to each violin plot are first the highest scoring known active and second the highest scoring molecule from running REvoLd on our custom chemical space. Lid_root2 is the an interface score normalized by the squared number of non-hydrogen atoms. All scoring results except EFR+REvoLd were taken from the original REvoLd publication

Fig. 16

Illustration of $a_{\mathit{ef}}$ in matrix representation for different joins to demonstrate the necessity of the embedded constraints. The first two constraints only allow one connection per row and column, which is infringed in the first example by $e_3$ and $f_3$, which both give rise to two newly formed edges. The third constraint aligns the join to impede multi-graph formation, as it is depicted in the second example. The fourth constraint ensures that the edge labels of introduced edges adhere to the definition of a weight preserving join, violated in the forth example