Flexible Topology: A Dynamic Model of a Continuous Chemical Space

Abstract

Ligand design problems involve searching chemical space for a molecule with a set of desired properties. As chemical space is discrete, this search must be conducted in a pointwise manner, separately investigating one molecule at a time, which can be inefficient. We propose a method called "Flexible Topology", where a ligand is composed of a set of shapeshifting "ghost" atoms, whose atomic identities and connectivity can dynamically change over the course of a simulation. Ghost atoms are guided toward their target positions using a translation-, rotation-, and index-invariant restraint potential. This is the first step toward a continuous model of chemical space, where a dynamic simulation can move from one molecule to another by following gradients of a potential energy function. This builds on a substantial history of alchemy in the field of molecular dynamics simulation, including the Lambda dynamics method developed by Brooks and co-workers [X. Kong and C.L. Brooks III, J. Chem. Phys. 105, 2414 (1996)], but takes it to an extreme by associating a set of four dynamical attributes with each shapeshifting ghost atom that control not only its presence but also its atomic identity. Here, we outline the theoretical details of this method, its implementation using the OpenMM simulation package, and some preliminary studies of ghost particle assembly simulations in vacuum. We examine a set of 10 small molecules, ranging in size from 6 to 50 atoms, and show that Flexible Topology is able to consistently assemble all of these molecules to high accuracy, beginning from randomly initialized positions and attributes.

Figure 1.

Framework of the Flexible Topology method. The gray and blue spheres show the ghost particles. MD simulations start with an initial structure of the protein of interest and randomly initialized ghost particles. The MLForce plugin computes external forces on the ghost particle positions and attributes using a loss function, $L.X_{\mathrm{G}\mathrm{P}}$ denotes the ghost particle positions, $A_{\mathrm{G}\mathrm{P}}$ represents the atomic attributes of ghost particles, $M$ shows the mapping function, and $F_{\mathrm{G}\mathrm{P}}$ and $F_{\mathrm{T}}$ show ghost particle and target ligand features, respectively. Final output shows the ghost particles assembled into a target structure in the binding pocket of the protein.

Figure 2.

Assembly simulation target molecules randomly chosen from the OpenChem dataset are shown in their 2D structures. The number of atoms for each molecule is written under it. The 2D structures are generated and drawn by MarvinSketch using implicit hydrogens.

Figure 3.

Time series analysis of the per particle loss (blue) and the assembly RMSD (gray) for successful assembly trajectories for molecule 8 (top), 20 (middle), and 40 (bottom). 2D structures are shown in Figure 2. A dashed line showing a loss-per-particle value of 0.002 is shown on each graph, which can be used across systems as an indicator of excellent assembly accuracy. The structures on the right show snapshots of the molecule at intermediate stages of assembly and correspond to the three black circles on the graphs on the left-hand side. The structures are shown in size-charge representation, where the radius of each particle roughly corresponds to its $\sigma$ value and the color shows its charge $(q)$, with blue being positive, red being negative, and white being neutral. The target structures show the conformers that were used to generate the target set of features that are almost indistinguishable from the 125 ps structures in each case.

Figure 4.

Performance of the Flexible Topology algorithm upon parameter perturbation. Each panel shows a set of box plots that show the minimum loss per particle achieved in a set of 10 independent assembly trajectories for each target. The boxes show the interquartile range, whiskers show 150% of the IQR, and outliers are shown as circles. The orange lines in each box mark the median datapoint of each set. A single parameter is explored in each panel: the number of steps per cycle (A), the periodicity of atom assignment (B), the value of the MLForce Scale (C), and the strength of the second-order bond restraints (D). All parameter values were perturbations from the central point of: 500 steps per cycle, 5000 step assignment period, 2000 MLForce Scale, and 0.5 bond restraint strength.