Molecular discovery by optimal sequential search

Abstract

In the development of a new compound in chemistry and molecular biology, especially a new medicine in pharmaceutical industry, we often need to find candidate(s), a molecule or molecules, with the best desired property (e.g., binding affinity in medicine) from a large set of molecules with the same scaffold but m distinct functional substitutes at each of its n different sites. The total number $N_{\mathrm{lib}}$ of molecules in this library is $m^n$. In some cases, $N_{\mathrm{lib}}$ can be a very large number (e.g., millions). This is a challenging task because it is costly and often infeasible to synthesize and test all of these molecules. A new algorithm referred to as optimal sequential search is developed to overcome this difficulty. Especially, this algorithm is chemically intuitive which only uses the information of molecule composition, and accessible to practical chemists. The algorithm can be applied to small, medium and large size molecule libraries. With syntheses and property measurements for a limited number of molecules, the top best candidate molecules can be effectively captured from the whole library. Three examples with library size 64, 160,000 and 1,048,576, respectively, are used for illustration. For the first small library, syntheses and property measurements of 17 molecules are sufficient to capture the top 7 best candidate molecules; for the two medium and large libraries, syntheses and property measurements of about one thousand molecules can capture most or a large part of the top 500, especially the top 100 best candidate molecules. However, the algorithm needs to perform multiple (e.g., hundreds of) iterative syntheses and property measurements. The time cost may not be acceptable if the algorithm is performed manually. To make the algorithm practical, automation of the sequential searching process is the following task.

Appendix: Deduction of the formula of the expected improvement E(I) for maximization

Define improvement $I(x_*)$ for an one variable function $y=f(x)$ as

$$\begin{aligned} I(x_*)=\max (Y-f_{\mathrm{max}}, 0). \end{aligned}$$

(26)

The expected improvement $E(I(x_*))$ is then

$$\begin{aligned} E(I(x_*))={\mathbb {E}}[I(x_*)]\equiv {\mathbb {E}}[\max (Y-f_{\mathrm{max}},0)] \end{aligned}$$

(27)

i.e.,

$$\begin{aligned} E(I(x_*))= & {} \int _{-\infty }^{\infty } \max (y-f_{\mathrm{max}},0)\phi (y){\mathrm d}y \nonumber \\= & {} \int _{-\infty }^{f_{\mathrm{max}}} 0\phi (y){\mathrm d}y + \int _{f_{\mathrm{max}}}^{\infty } (y-f_{\mathrm{max}})\phi (y){\mathrm d}y \nonumber \\= & {} \int _{f_{\mathrm{max}}}^{\infty }(y-f_{\mathrm{max}})\frac{1}{\sqrt{2\pi }s} \exp \left[ -\frac{(y-f_*)^2}{2s^2}\right] {\mathrm d}y, \end{aligned}$$

(28)

where $\phi (y)$ is the probability density function (pdf) of normal distribution ${{\mathcal {N}}}(f_*, s^2)$. Setting

$$\begin{aligned} u=\frac{y-f_*}{s}, \qquad y=su+f_*, \qquad {\mathrm d}y = s{\mathrm d}u \end{aligned}$$

(29)

gives

$$\begin{aligned} E(I(x_*))= & {} \int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } (su+f_*-f_{\mathrm{max}}) \frac{1}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u \nonumber \\= & {} (f_*-f_{\mathrm{max}})\int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } \frac{1}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u + s\int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } \frac{u}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u \nonumber \\= & {} (f_*-f_{\mathrm{max}})\left[ \int _{-\infty }^{\infty } \frac{1}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u -\int _{-\infty }^{(f_{\mathrm{max}}-f_*)/s} \frac{1}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u\right] \nonumber \\+ & {} s\int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } \frac{u}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u \nonumber \\= & {} (f_*-f_{\mathrm{max}})\left[ 1-\Phi \left( \frac{f_\mathrm{max}-f_*}{s}\right) \right] + s\int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } \frac{u}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u, \nonumber \\ \end{aligned}$$

(30)

where $\Phi$ denotes the standard normal cumulative distribution. Note that

$$\begin{aligned} \phi (u)=\frac{1}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}, \qquad {\mathrm d}\phi (u)=-\frac{u}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u. \end{aligned}$$

(31)

The second term in Eq. (30) becomes

$$\begin{aligned} s\int _{(f_\mathrm{max}-f_*)/s}^{\infty } \!\!\!\!\!\!\!&\frac{u}{\sqrt{2\pi }}e^{-\frac{u^2}{2}}{\mathrm d}u =-s \int _{(f_{\mathrm{max}}-f_*)/s}^{\infty } {\mathrm d}\phi \nonumber \\= & {} s\left[ \phi (u)\right] _{\infty }^{(f_{\mathrm{max}}-f_*)/s} =s\phi \left( \frac{f_{\mathrm{max}}-f_*}{s}\right) . \end{aligned}$$

(32)

Substituting Eq. (32) into (30) yields

$$\begin{aligned} E(I(x_*))=(f_*-f_{\mathrm{max}})\left[ 1-\Phi \left( \frac{f_\mathrm{max}-f_*}{s}\right) \right] +s\phi \left( \frac{f_\mathrm{max}-f_*}{s}\right) . \end{aligned}$$

(33)

Note that in the whole process, only the output y is involved. For $f(\mathbf{x})$ of $n(>1)$-dimensional variables, the result is the same. Equation (33) is then also valid for $f(\mathbf{x})$, i.e.,

$$\begin{aligned} E(I(\mathbf{x}_*))=(f_*-f_{\mathrm{max}})\left[ 1-\Phi \left( \frac{f_\mathrm{max}-f_*}{s}\right) \right] +s\phi \left( \frac{f_\mathrm{max}-f_*}{s}\right) . \end{aligned}$$

(34)