EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2021-058 LHCb-PAPER-2020-047 May 4, 2021

Measurement of $C\\!P$ asymmetry

in ${D}^{0}$ $\rightarrow$ ${K}^{0}_{\mathrm{S}}$ ${K}^{0}_{\mathrm{S}}$
decays

LHCb collaboration†††Authors are listed at the end of this Letter.

A measurement of the $C\\!P$ asymmetry in
${{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decays is
reported, based on a data sample of proton-proton collisions collected by the
LHCb experiment from 2015 to 2018, corresponding to an integrated luminosity
of $6\text{\,fb}^{-1}$. The flavor of the ${D}^{0}$ candidate is determined
using the charge of the ${D}^{*\pm}$ meson, from which the decay is required
to originate. The ${{D}^{0}}\rightarrow{{K}^{+}}{{K}^{-}}$ decay is used as a
calibration channel. The time-integrated $C\\!P$ asymmetry for the
${{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ mode is
measured to be:

${{\mathcal{A}}^{{C\\!P}}}({{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=(-3.1\pm
1.2\pm 0.4\pm 0.2)\%,$

where the first uncertainty is statistical, the second is systematic, and the
third is due to the uncertainty on the $C\\!P$ asymmetry of the calibration
channel. This is the most precise determination of this quantity to date.

Submitted to Phys. Rev. D Lett.

© 2024 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence.

The existence of charge-parity violation ($C\\!PV$) effects in charm hadrons
has recently been established [1], which constitutes the only evidence of
$C\\!PV$ in the up-quark ($u$, $c$ or $t$) sector. However, current
experimental evidence is limited to a single observable,
$\Delta{{\mathcal{A}}^{{C\\!P}}}$, the difference between the $C\\!P$
asymmetry in ${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays. Given the
uncertainties in the theoretical predictions on $C\\!PV$ in the charm-quark
sector within the Standard Model (SM), it is currently not possible to test
these [2, 3]. Further measurements in charm-hadron decays are crucial to shed
light on $C\\!PV$ phenomenology. This could involve dynamics beyond the SM,
which is not constrained to be the same as in the down-quark ($d$, $s$ or $b$)
sector and could enter the amplitudes via loop contributions, affecting
observables in a detectable way.

Amongst many charm-hadron decay modes in which $C\\!PV$ could manifest, the
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ mode is
a promising target because of the expected size of the effect. Its $C\\!P$
asymmetry is defined by

${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=\dfrac{\Gamma({{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})-\Gamma({{\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})}{\Gamma({{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})+\Gamma({{\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})},$
(1)

where $\Gamma$ is the decay width of the ${D}^{0}$ (${\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{0}$) meson, and it can be larger than
in other channels, up to the percent level in the SM [4, 5, 6, 7, 8]. In fact,
only amplitudes proceeding via loop-suppressed and tree-level exchange
diagrams, which vanish in the flavor–SU(3) limit, contribute to this decay,
and they are similar in size. Their interference could therefore result in a
detectable $C\\!P$ asymmetry. In addition, the $C\\!P$ asymmetry in the
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decay is
sensitive to a different mix of amplitudes compared to
${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ and
${{D}^{0}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays. Therefore, measuring
${\mathcal{A}}^{{C\\!P}}$ (${K}^{0}_{\mathrm{S}}$ ${K}^{0}_{\mathrm{S}}$)
provides independent information which can help to elucidate the mechanisms of
$C\\!PV$ in charm hadron decays.

The current world average for the time-integrated $C\\!P$ asymmetry is
${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=(0.4\pm
1.4)\%$[9], the precision of which is still insufficient for detecting an
effect. In this work, a new measurement of this quantity, performed with
proton-proton ($pp$) collisions collected from 2015 to 2018 (Run 2) by the
LHCb experiment at the LHC at CERN, is reported. Data collected at a center-
of-mass energy of 13 TeV, corresponding to an integrated luminosity of
$6\,\text{\,fb}^{-1}$, are used. The subsample of data collected during 2015
and 2016 has been analyzed previously [10]. In the present work a number of
improvements in the analysis leads to a more efficient selection and a
sizeable improvement in sensitivity.

The measurement of
${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})$
requires knowledge of the ${D}^{0}$ flavor at production. A sample of flavor-
tagged ${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$
decays is obtained by selecting only ${D}^{0}$ mesons that originate from
${{D}^{*+}}\rightarrow{{D}^{0}}{{\pi}^{+}}$ decays.111Inclusion of the charge-
conjugate process is implied throughout this document unless explicitly
specified. The charge of the pion (tagging pion) in this decay identifies the
flavor of the accompanying ${D}^{0}$ meson. While ${D}^{0}$ oscillations may
cause some of them to change flavor before decaying, this is a small effect in
comparison to the resolution of the current measurement and is not considered
further.

The decay widths $\Gamma$ in Eq. 1 are related to the number of observed
candidates $N$ by

$N(\kern 1.79993pt\shortstack{{(\rule[2.15277pt]{12.50002pt}{0.51212pt})} \\\
[-.7ex] $\kern-1.79993ptD$}^{\rm
0}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})\propto\sigma({{D}^{*\pm}}){\varepsilon}^{\pm}\Gamma(\kern
1.79993pt\shortstack{{(\rule[2.15277pt]{12.50002pt}{0.51212pt})} \\\ [-.7ex]
$\kern-1.79993ptD$}^{\rm
0}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}),$ (2)

where $\sigma({{D}^{*\pm}})$ are the production cross sections and
${\varepsilon}^{\pm}$ the detection probabilities for ${D}^{*\pm}$ decays.
Both factors are charge-asymmetric, due to the ${D}^{*\pm}$ production
asymmetry arising from the hadronization of charm quarks in $pp$ collisions,
and to asymmetries in the geometry and response of the detector. Because of
the charge-symmetric final state of the ${D}^{0}$ decay, the only source of
detection asymmetry comes from the tagging pion. All these quantities are
calibrated using a large sample of ${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$
decays, for which ${\mathcal{A}}^{{C\\!P}}$ is known with much greater
precision than in the
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decay
mode [11].

The LHCb detector is a single-arm forward spectrometer designed for the study
of particles containing $b$ or $c$ quarks, as described in detail in Refs.
[12, 13]. The LHCb tracking system exploits a dipole magnet to measure the
momentum of charged particles, and it consists of a silicon-strip vertex
detector surrounding the $pp$ interaction region [14] and tracking stations
placed upstream and downstream of the magnet. The magnetic-field polarity is
reversed periodically during data taking, alternating between pointing upwards
(MagUp) and downwards (MagDown), to mitigate the differences of reconstruction
efficiencies of particles with opposite charges. The online event selection is
performed by a trigger [15], which consists of a hardware stage, based on
information from the calorimeter and muon systems, followed by a software
stage, which applies a full event reconstruction and selection in two steps.
The simulated data used in this analysis are produced using the software
described in Refs. [16, 17, 18, 19, 20].

For an event to be considered in this analysis, the hardware-trigger decision
is required to have been based either on the transverse energy deposited in
the hadronic calorimeter by a particle associated with the ${{D}^{0}}$
candidate, or on signatures not associated with the ${{D}^{*+}}$ candidate,
such as a muon with high transverse momentum, or a large transverse-energy
deposit in the electromagnetic or hadronic calorimeter from the rest of the
event. Similarly, the first software stage decision is required to have been
based either on the presence of a single track associated with the ${D}^{0}$
decay, with sufficient transverse momentum and impact parameter relative to
the primary $p$ $p$ collision vertex (PV), or on track-related signatures
independent of the ${{D}^{*+}}$ candidate. This ensures that the trigger
decision is independent of the tagging pion, preventing a possible charge-
dependent bias which would be difficult to measure and correct accurately.

At the second software stage, a complete reconstruction of the decay chain is
performed, requiring a pair of ${K}^{0}_{\mathrm{S}}$ mesons compatible with
the decay of a ${D}^{0}$ particle, and the presence of an additional pion,
compatible with originating from a
${{D}^{*+}}\\!\rightarrow{{D}^{0}}{{\pi}^{+}}$ decay. Candidate
${{K}^{0}_{\mathrm{S}}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays are
classified in two categories: those in which the decay occurred early enough
for the pions to be reconstructed in the vertex detector (“long”, indicated in
the following by an “L” label); and those decaying later, such that segments
of pion tracks can only be formed in the silicon-strip tracker upstream of the
magnet, and the tracking stations after the magnet (“downstream”, indicated in
the following by a “D” label). The geometric acceptance of the LHCb detector
for typical ${K}^{0}_{\mathrm{S}}$ momenta is about a factor of two larger for
D-type versus L-type decays, but the mass, momentum and vertex resolution of
the latter category are better than those of the former.

Candidates are separated into three categories: LL, LD, or DD, according to
the types of the two ${K}^{0}_{\mathrm{S}}$ mesons associated to the ${D}^{0}$
meson. These are analyzed separately, and the results are combined only at the
end. To prevent any experimenters’ bias, the flavor of ${D}^{0}$ candidates
was not examined until the analysis was finalized.

Some of the ${D}^{*+}$ mesons are not promptly produced in the primary $p$ $p$
interaction, but rather come from the decay of a beauty-flavored hadron. These
secondary decays are affected by different production asymmetries with respect
to prompt decays, which is estimated at the 1% level. No attempt is made in
this analysis to reject secondary decays, as selection requirements based on
pointing observables are less effective in ${D}^{0}$ mesons reconstructed with
displaced vertices, such as
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$,
compared to other charm hadron decays. The spurious asymmetry introduced by
this contribution is instead canceled by use of a dedicated calibration
sample. For this purpose, a sample of
${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ decays produced in ${D}^{*+}$
decays is selected in a way to ensure it contains prompt and secondary decays
in the same proportions as the
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ sample.
Its selection, starting from the first trigger level, refrains from any
requirements for which the effect on the secondary-to-prompt ratio cannot
reliably reproduce that in the signal sample. This requires, amongst other
things, to avoid selections based on observables sensitive to the location of
the ${D}^{*}$ or ${D}^{0}$ vertices, which have different resolutions in the
two samples [21].

The ${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$
signal peaks in the difference of invariant masses ${\mbox{$\Delta
m$}}=m({{D}^{*+}})-m({{D}^{0}})$, where $m({{D}^{0}})$ is the invariant mass
of the ${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$
candidate and $m({{D}^{*+}})$ is the invariant mass of the
${{D}^{*+}}\\!\rightarrow{{D}^{0}}{{\pi}^{+}}$ candidate. This sample may be
contaminated by other physics processes. Amongst them, partially reconstructed
${D}^{0}$ decays coming from ${D}^{*+}$ mesons, such as
${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}{{\pi}^{0}}$,
can also peak in $\Delta m$. Some non-${D}^{0}$ decays can also contaminate
the sample, the main contributor being
${{D}^{+}_{s}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}{{\pi}^{+}}$
decays, where the ${{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ pair is
incorrectly associated with a ${{D}^{0}}$ decay. These background components
are effectively suppressed by selecting only candidates with $m({{D}^{0}})$
within $\pm 20(30)\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ of the known ${D}^{0}$
mass[9] in LL/LD (DD) samples.

A background source which is difficult to eliminate completely is the abundant
${D}^{0}$ $\rightarrow$ ${K}^{0}_{\mathrm{S}}$ ${\pi}^{+}$ ${\pi}^{-}$ decay,
with the ${{\pi}^{+}}{{\pi}^{-}}$ pair mimicking a ${{K}^{0}_{\mathrm{S}}}$.
After a loose selection based on the flight distance of the
${K}^{0}_{\mathrm{S}}$, it is statistically separated from the signal through
a simultaneous fit of three observables: the $\Delta m$ and the invariant
masses, $m({{K}^{0}_{\mathrm{S}}})$, of the two ${{K}^{0}_{\mathrm{S}}}$
candidates.

To maximize the sensitivity of the measurement, the LL and LD subsamples are
further split according to whether the ${D}^{*+}$ candidate is compatible or
incompatible with having originated from the PV. Compatibility is defined by a
fixed threshold on the goodness of fit of the complete decay chain [22], when
constrained to be originating from the PV. For PV-compatible candidates,
independent of their true nature, a factor of two better mass resolution is
obtained by constraining the origin vertex to coincide with the PV. Therefore,
separate analyses of the two subsamples, followed by the combination of the
results provides optimal exploitation of the available data. This is not done
for the DD sample due to its limited signal yield. While this requirement
might in principle have slightly different effects on the secondary fraction
of signal and calibration sample, differences are negligible, as it has been
explicitly verified on data.

The data are further split according to a classifier sensitive to the signal-
to-background ratio. The classifier combines a number of track-related
observables, including track and vertex quality, transverse momenta of
${K}^{0}_{\mathrm{S}}$ and ${D}^{0}$ candidates, helicity angles of the
${K}^{0}_{\mathrm{S}}$ and ${D}^{0}$ decays, and particle identification
parameters of the ${D}^{0}$ final state particles. The kinematics of the
tagging pion is excluded to avoid introducing possible charge asymmetry
biases. All these observables are combined in a k-nearest-neighbors (kNN)
classifier [23], trained using a simulated signal sample, and data from the
${D}^{0}$ mass sidebands for background. The resulting discriminant is used to
split each sample in three categories (low-, medium- and high-purity); the
low-purity class contains very little signal and is removed. The thresholds
separating the three categories are numerically optimized in a way to achieve
the best combined resolution on ${\mathcal{A}}^{{C\\!P}}$ [21].

Finally, 2015–2016 data are analyzed separately from 2017–2018 data, due to
different trigger conditions. In addition, DD candidates were not collected in
2015–2016 data taking. In about 10% of cases, multiple ${D}^{*+}$ candidates
are found in the same event. This is mainly due to ${D}^{0}$ candidates that
are associated with multiple tagging pions. In these cases, one ${D}^{*+}$
candidate is randomly selected for the analysis.

Distributions of the $\Delta m$ variable for some representative subsamples
are shown in Fig. 1.

Figure 1: Distributions and fit projections of the $\Delta m$ observable for
representative candidate categories (2017–2018 data). See text for the
definition of purities.

An unbinned maximum-likelihood fit is performed to the joint distribution of
$\Delta m$ and the two $m({{K}^{0}_{\mathrm{S}}})$ observables, simultaneously
to candidates of both flavors to obtain the value of
${\mathcal{A}}^{{C\\!P}}$. The total probability density function is
parameterized by the sum of eight components: the signal component, peaking in
the three observables, and seven additional components, each describing a
specific background source. This includes ${{D}^{0}}\rightarrow\
{{K}^{0}_{\mathrm{S}}}{{\pi}^{+}}{{\pi}^{-}}$ decays, which peak in $\Delta m$
and in one $m({{K}^{0}_{\mathrm{S}}})$ distribution but not in the other, and
all possible combinations of unrelated particles. The peaking component in the
$\Delta m$ distribution is described by a Gaussian function in PV-incompatible
samples, or a Johnson $S_{U}$ distribution [24] in PV-compatible samples. The
peaking component in the $m({{K}^{0}_{\mathrm{S}}})$ distribution is described
by the sum of two Gaussian functions with different widths and a common mean,
for both L- and D-category ${K}^{0}_{\mathrm{S}}$ candidates. The non-peaking
component in the $\Delta m$ distribution is described by an empirical
threshold function, while in the $m({{K}^{0}_{\mathrm{S}}})$ distribution it
is described by a first- and second-order Chebyshev polynomial for L- and
D-category ${K}^{0}_{\mathrm{S}}$ candidates, respectively. In each subsample,
the parameters defining the signal and background probability density
functions are shared between flavors, while the normalization of each
component is allowed to differ.

Each candidate participating in the fit is appropriately weighted with the aim
of correcting for all spurious asymmetries, with the help of the calibration
${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ sample, selected in a way to
contain the same proportions of primary and secondary decays. The calibration
sample is similarly split between PV-compatible and PV-incompatible
categories, and is required to have $m({{K}^{+}}{{K}^{-}})$ within $\pm
20\text{\,Me\kern-1.00006ptV\\!/}c^{2}$ around the known ${{D}^{0}}$ mass. In
addition, only candidates in a $\pm 1.5\text{\,Me\kern-1.00006ptV\\!/}c^{2}$
range around the $\Delta m$ peak are used.

To compute weights, each ${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ subsample
is classified in categories of observed charge asymmetry, by a kNN classifier
based on the space of kinematic parameters of the ${D}^{0}$ candidate. The
detector and production asymmetries observed in the calibration sample are
independent of the ${D}^{0}$ decay mode, being charge-symmetric, and can be
used to correct the signal sample for the same effects. This is achieved by
weighting each signal candidate in the global fit by a charge-dependent factor

$w^{\pm}(\vec{p}_{0})=\frac{n_{C}^{+}(\vec{p}_{0})+n_{C}^{-}(\vec{p}_{0})}{2n_{C}^{\pm}(\vec{p}_{0})}\left[1\pm{{\mathcal{A}}^{{C\\!P}}}({{K}^{+}}{{K}^{-}})\right],$
(3)

where $\vec{p}_{0}$ is the ${D}^{0}$ three-momentum, and
$n_{C}^{\pm}(\vec{p}_{0})$ is the density of calibration ${D}^{*\pm}$ decays
in the $\vec{p}_{0}$ space[21].

To account for a possible dependence of the detection asymmetry on magnet
polarity, weights are separately calculated for MagUp and MagDown
configurations. Their distributions are shown in Fig. 2, where their
difference is clearly visible. To avoid weights affected by large
uncertainties, a negligible fraction of candidates having very large weights
($>10$) are dropped at this stage. The weighting procedure described above is
a better alternative to the procedure of binning the ${D}^{0}$ kinematic space
and weighting the two samples to have the same distribution. The advantage is
a reduction of dimensionality to the single variable of actual relevance, that
is the charge asymmetry to be corrected. This reduces the loss in statistical
power that occurs when weighting bins of very different populations in the two
samples under comparison. In addition, it holds exactly even when the involved
asymmetries are not small. This is helpful in the case of the LHCb detector,
where some kinematic regions are characterized by sizeable detection
asymmetries for charged pions [10]. This affects the tagging pion in 30% of
our events.

Figure 2: Distributions of weights (2017–2018 data) applied to signal
candidates in the global fit to correct for detector, production, and physics
of the calibration channel. Differences between MagUp (left) and MagDown
(right) are a consequence of different ${D}^{0}$ and ${\kern
1.79993pt\overline{\kern-1.79993ptD}}{}^{0}$ acceptances, and of detector
asymmetries. A small fraction of entries with large weights fall outside the
range of the plot.

This weighting method has been extensively checked with both data and
simulation. One of the checks is to apply the weighting procedure to half of
the ${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ sample, using the second half
to calculate weights. This has the expected effect of canceling the asymmetry
of the sample, that was initially highly significant, $(1.3\pm 0.1)\%$, due to
a combination of physics and detector effects.

The systematic uncertainty associated to the limited size of the calibration
mode is determined by a bootstrap sampling of the data, and is found to be
negligible. The uncertainty due to the presence of residual background in the
${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ sample has also been evaluated and
found to be a minor effect.

Another important source of systematic uncertainty is due to the limited
knowledge of the shape of the mass distributions, as the fit function used is
purely empirical. Any potential bias in the procedure is evaluated by fitting
the model to simulated samples of pseudo-experiments from alternative models.
The number and type of peaking components and the behavior of the background
at the kinematic threshold are varied, and the largest observed variations are
taken as the systematic uncertainty, ranging between 0.4% and 5.6%, depending
on the subsample.

A systematic uncertainty is also assessed for possible residual differences in
secondary decay fractions between signal and calibration samples. The largest
discrepancy between the two samples is the presence of slightly different
trigger requirements on the proper decay time of the ${D}^{0}$ candidate. This
has been corrected for by increasing the weight of candidates close to the
threshold, to emulate the effect of those that have been lost. An uncertainty
of this correction is conservatively assessed, that is between 0.1% and 0.2%.
Finally, an uncertainty of $0.15\%$ in the input value of
${{\mathcal{A}}^{{C\\!P}}}({{K}^{+}}{{K}^{-}})$ is separately taken into
account, measured by the LHCb experiment as
${{\mathcal{A}}^{{C\\!P}}}({{K}^{+}}{{K}^{-}})=(0.04\pm 0.12\pm 0.10)\%$ [11].

The results obtained in each subsample are summarized in Table 1. The
subsamples corresponding to different kNN classifier ranges are fitted
simultaneously with common parameters, so they do not produce separate
results. As a check of goodness of fit, a $\chi^{2}$ has been calculated for
each one-dimensional projection of the fit (Fig. 1). All $p$-values are found
to be greater than $0.2$.

Table 1: Measurements of yields and ${{\mathcal{A}}^{{C\\!P}}}({{D}^{0}}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})$ in individual subsamples. For asymmetries, the first uncertainty is statistical and the second is systematic. Sample | 2015 + 2016 (2 $\text{\,fb}^{-1}$) | 2017 + 2018 (4 $\text{\,fb}^{-1}$)
---|---|---
| Yield |  ${\mathcal{A}}^{{C\\!P}}$ [%] | Yield |  ${\mathcal{A}}^{{C\\!P}}$ [%]
LL PV-comp. | $1388\pm 41$ | $\phantom{(}\phantom{0}\phantom{-}0.3\pm\phantom{0}2.5\pm 0.6$ | $4056\pm 77$ | $\phantom{0}-4.3\pm\phantom{0}1.6\pm 0.4$
LL PV-incomp. | $\phantom{0}178\pm 31$ | $\phantom{}-11\phantom{.0}\pm 17\phantom{.0}\pm 2\phantom{.0}$ | $\phantom{0}430\pm 41$ | $\phantom{0}-3.0\pm\phantom{0}7.9\pm 1.1$
LD PV-comp. | $\phantom{0}411\pm 25$ | $\phantom{0}-7.2\pm\phantom{0}5.8\pm 1.1$ | $1145\pm 49$ | $\phantom{0}-2.9\pm\phantom{0}3.8\pm 0.7$
LD PV-incomp. | $\phantom{00}58\pm 18$ | $\phantom{}-10\phantom{.0}\pm 31\phantom{.0}\pm 4\phantom{.0}$ | $\phantom{0}349\pm 64$ | $\phantom{0}-5\phantom{.0}\pm 17\phantom{.0}\pm 2\phantom{.0}$
DD | $-$ | $-$ | $\phantom{00}87\pm 28$ | $\phantom{}-35\phantom{.0}\pm 47\phantom{.0}\pm 6\phantom{.0}$

All partial results in Table 1 are statistically compatible with each other.
The weighted average of all measurements using 2015–2016 data, is

${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=(-1.1\pm
2.3\pm 0.5\pm 0.2)\%,$

where the first uncertainty is statistical, the second is systematic, obtained
by taking the individual contributions as uncorrelated, and the third is from
the uncertainty on ${{\mathcal{A}}^{{C\\!P}}}({{K}^{+}}{{K}^{-}})$. This
result is compatible with the previously published value based on the same
data sample [10], but has a better precision by about $30$%, corresponding to
an effective doubling of the yields. The sensitivity increase is due to the
combined effect of several improvements in the analysis, most notably the new
weighting technique, the inclusion of secondary decays, an appropriate
categorization of the sample, and the multidimensional likelihood fit.

The asymmetry for the 2017–2018 data is measured to be

${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=(-4.0\pm
1.5\pm 0.3\pm 0.2)\%.$

Treating the systematic uncertainties as uncorrelated between the data taking
periods, except for the shape of the fit functions that is considered to be
fully correlated, the asymmetry combining the results from the full Run 2 data
sample is obtained as

${{\mathcal{A}}^{{C\\!P}}}({{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}})=(-3.1\pm
1.2\pm 0.4\pm 0.2)\%.$

This measurement supersedes the previous LHCb result [10] and is in agreement
with all previous determinations [25, 26, 27]. It is the most precise
measurement of this quantity to date, and it is compatible with no $C\\!P$
asymmetry at the level of $2.4$ standard deviations.

## Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments
for the excellent performance of the LHC. We thank the technical and
administrative staff at the LHCb institutes. We acknowledge support from CERN
and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST
and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN
(Italy); NWO (Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MSHE
(Russia); MICINN (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC
(United Kingdom); DOE NP and NSF (USA). We acknowledge the computing resources
that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN
(Italy), SURF (Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and
Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-
GRID (Poland) and NERSC (USA). We are indebted to the communities behind the
multiple open-source software packages on which we depend. Individual groups
or members have received support from ARC and ARDC (Australia); AvH Foundation
(Germany); EPLANET, Marie Skłodowska-Curie Actions and ERC (European Union);
A*MIDEX, ANR, IPhU and Labex P2IO, and Région Auvergne-Rhône-Alpes (France);
Key Research Program of Frontier Sciences of CAS, CAS PIFI, CAS CCEPP,
Fundamental Research Funds for the Central Universities, and Sci. & Tech.
Program of Guangzhou (China); RFBR, RSF and Yandex LLC (Russia); GVA, XuntaGal
and GENCAT (Spain); the Leverhulme Trust, the Royal Society and UKRI (United
Kingdom).

## References

  * [1] LHCb collaboration, R. Aaij et al., _Observation of $C\\!P$ violation in charm decays_, Phys. Rev. Lett. 122 (2019) 211803, arXiv:1903.08726
  * [2] M. Chala et al., _$\Delta A_{CP}$ within the Standard Model and beyond_, JHEP 07 (2019) 161, arXiv:1903.10490
  * [3] Y. Grossman and S. Schacht, _The emergence of the $\Delta U=0$ rule in charm physics_, JHEP 07 (2019) 020, arXiv:1903.10952
  * [4] U. Nierste and S. Schacht, _CP violation in $D^{0}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$_, Phys. Rev. D92 (2015) 054036, arXiv:1508.00074
  * [5] H.-N. Li, C.-D. Lu, and F.-S. Yu, _Branching ratios and direct CP asymmetries in $D\rightarrow PP$ decays_, Phys. Rev. D86 (2012) 036012, arXiv:1203.3120
  * [6] H.-Y. Cheng and C.-W. Chiang, _Revisiting CP violation in $D\rightarrow P\\!P$ and $V\\!P$ decays_, Phys. Rev. D100 (2019) 093002, arXiv:1909.03063
  * [7] F. Buccella, A. Paul, and P. Santorelli, _$SU(3)_{F}$ breaking through final state interactions and $CP$ asymmetries in $D\rightarrow PP$ decays_, Phys. Rev. D99 (2019) 113001, arXiv:1902.05564
  * [8] J. Brod, A. L. Kagan, and J. Zupan, _Size of direct CP violation in singly Cabibbo-suppressed D decays_ , Phys. Rev. D86 (2012) 014023, arXiv:1111.5000
  * [9] Particle Data Group, P. A. Zyla et al., _Review of particle physics_ , Prog. Theor. Exp. Phys. 2020 (2020) 083C01
  * [10] LHCb collaboration, R. Aaij et al., _Measurement of the time-integrated $C\\!P$ asymmetry in ${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decays_, JHEP 11 (2018) 048, arXiv:1806.01642
  * [11] LHCb collaboration, R. Aaij et al., _Measurement of $C\\!P$ asymmetry in ${{D}^{0}}\\!\rightarrow{{K}^{+}}{{K}^{-}}$ decays_, Phys. Lett. B767 (2017) 177, arXiv:1610.09476
  * [12] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [13] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [14] R. Aaij et al., _Performance of the LHCb Vertex Locator_ , JINST 9 (2014) P09007, arXiv:1405.7808
  * [15] R. Aaij et al., _The LHCb trigger and its performance in 2011_, JINST 8 (2013) P04022, arXiv:1211.3055
  * [16] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [17] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [18] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [19] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [20] M. Clemencic et al., _The LHCb simulation application, Gauss: design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [21] G. Tuci, _Searching for confirmation of charm $C\\!P$ violation in ${K}^{0}_{\mathrm{S}}$ final states at LHCb_,  2021, PhD thesis, University of Pisa, Pisa, CERN-THESIS-2020-325
  * [22] W. D. Hulsbergen, _Decay chain fitting with a Kalman filter_ , Nucl. Instrum. Meth. A552 (2005) 566, arXiv:physics/0503191
  * [23] J. Friedman, T. Hastie, and R. Tibshirani, The Elements of Statistical Learning, Springer Series in Statistics, 2001
  * [24] N. L. Johnson, _Systems of frequency curves generated by methods of translation_ , Biometrika 36 (1949) 149
  * [25] CLEO collaboration, G. Bonvicini et al., _Search for CP violation in $D^{0}\rightarrow{{K}^{0}_{\mathrm{S}}}\pi^{0}$ and $D^{0}\rightarrow\pi^{0}\pi^{0}$ and $D^{0}\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decays_, Phys. Rev. D63 (2001) 071101, arXiv:hep-ex/0012054
  * [26] N. Dash et al., _Search for $C\\!P$ violation and measurement of the branching fraction in the decay ${D}^{0}$ $\rightarrow$ ${K}^{0}_{\mathrm{S}}$ ${K}^{0}_{\mathrm{S}}$_, Phys. Rev. Lett. 119 (2017) 171801, arXiv:1705.05966
  * [27] LHCb collaboration, R. Aaij et al., _Measurement of the time-integrated $C\\!P$ asymmetry in ${{D}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{S}}}{{K}^{0}_{\mathrm{S}}}$ decays_, JHEP 10 (2015) 055, arXiv:1508.06087

LHCb collaboration

R. Aaij32, C. Abellán Beteta50, T. Ackernley60, B. Adeva46, M. Adinolfi54, H.
Afsharnia9, C.A. Aidala85, S. Aiola25, Z. Ajaltouni9, S. Akar65, J.
Albrecht15, F. Alessio48, M. Alexander59, A. Alfonso Albero45, Z. Aliouche62,
G. Alkhazov38, P. Alvarez Cartelle55, S. Amato2, Y. Amhis11, L. An48, L.
Anderlini22, A. Andreianov38, M. Andreotti21, F. Archilli17, A. Artamonov44,
M. Artuso68, K. Arzymatov42, E. Aslanides10, M. Atzeni50, B. Audurier12, S.
Bachmann17, M. Bachmayer49, J.J. Back56, S. Baker61, P. Baladron Rodriguez46,
V. Balagura12, W. Baldini21,48, J. Baptista Leite1, R.J. Barlow62, S.
Barsuk11, W. Barter61, M. Bartolini24, F. Baryshnikov82, J.M. Basels14, G.
Bassi29, B. Batsukh68, A. Battig15, A. Bay49, M. Becker15, F. Bedeschi29, I.
Bediaga1, A. Beiter68, V. Belavin42, S. Belin27, V. Bellee49, K. Belous44, I.
Belov40, I. Belyaev41, G. Bencivenni23, E. Ben-Haim13, A. Berezhnoy40, R.
Bernet50, D. Berninghoff17, H.C. Bernstein68, C. Bertella48, A. Bertolin28, C.
Betancourt50, F. Betti20,d, Ia. Bezshyiko50, S. Bhasin54, J. Bhom35, L.
Bian73, M.S. Bieker15, S. Bifani53, P. Billoir13, M. Birch61, F.C.R. Bishop55,
A. Bitadze62, A. Bizzeti22,k, M. Bjørn63, M.P. Blago48, T. Blake56, F.
Blanc49, S. Blusk68, D. Bobulska59, J.A. Boelhauve15, O. Boente Garcia46, T.
Boettcher64, A. Boldyrev81, A. Bondar43, N. Bondar38,48, S. Borghi62, M.
Borisyak42, M. Borsato17, J.T. Borsuk35, S.A. Bouchiba49, T.J.V. Bowcock60, A.
Boyer48, C. Bozzi21, M.J. Bradley61, S. Braun66, A. Brea Rodriguez46, M.
Brodski48, J. Brodzicka35, A. Brossa Gonzalo56, D. Brundu27, A. Buonaura50, C.
Burr48, A. Bursche27, A. Butkevich39, J.S. Butter32, J. Buytaert48, W.
Byczynski48, S. Cadeddu27, H. Cai73, R. Calabrese21,f, L. Calefice15,13, L.
Calero Diaz23, S. Cali23, R. Calladine53, M. Calvi26,j, M. Calvo Gomez84, P.
Camargo Magalhaes54, A. Camboni45,84, P. Campana23, A.F. Campoverde Quezada6,
S. Capelli26,j, L. Capriotti20,d, A. Carbone20,d, G. Carboni31, R.
Cardinale24,h, A. Cardini27, I. Carli4, P. Carniti26,j, L. Carus14, K.
Carvalho Akiba32, A. Casais Vidal46, G. Casse60, M. Cattaneo48, G.
Cavallero48, S. Celani49, J. Cerasoli10, A.J. Chadwick60, M.G. Chapman54, M.
Charles13, Ph. Charpentier48, G. Chatzikonstantinidis53, C.A. Chavez
Barajas60, M. Chefdeville8, C. Chen3, S. Chen27, A. Chernov35, V. Chobanova46,
S. Cholak49, M. Chrzaszcz35, A. Chubykin38, V. Chulikov38, P. Ciambrone23,
M.F. Cicala56, X. Cid Vidal46, G. Ciezarek48, P.E.L. Clarke58, M. Clemencic48,
H.V. Cliff55, J. Closier48, J.L. Cobbledick62, V. Coco48, J.A.B. Coelho11, J.
Cogan10, E. Cogneras9, L. Cojocariu37, P. Collins48, T. Colombo48, L.
Congedo19,c, A. Contu27, N. Cooke53, G. Coombs59, G. Corti48, C.M. Costa
Sobral56, B. Couturier48, D.C. Craik64, J. Crkovská67, M. Cruz Torres1, R.
Currie58, C.L. Da Silva67, E. Dall’Occo15, J. Dalseno46, C. D’Ambrosio48, A.
Danilina41, P. d’Argent48, A. Davis62, O. De Aguiar Francisco62, K. De
Bruyn78, S. De Capua62, M. De Cian49, J.M. De Miranda1, L. De Paula2, M. De
Serio19,c, D. De Simone50, P. De Simone23, J.A. de Vries79, C.T. Dean67, D.
Decamp8, L. Del Buono13, B. Delaney55, H.-P. Dembinski15, A. Dendek34, V.
Denysenko50, D. Derkach81, O. Deschamps9, F. Desse11, F. Dettori27,e, B.
Dey73, P. Di Nezza23, S. Didenko82, L. Dieste Maronas46, H. Dijkstra48, V.
Dobishuk52, A.M. Donohoe18, F. Dordei27, A.C. dos Reis1, L. Douglas59, A.
Dovbnya51, A.G. Downes8, K. Dreimanis60, M.W. Dudek35, L. Dufour48, V. Duk77,
P. Durante48, J.M. Durham67, D. Dutta62, M. Dziewiecki17, A. Dziurda35, A.
Dzyuba38, S. Easo57, U. Egede69, V. Egorychev41, S. Eidelman43,v, S.
Eisenhardt58, S. Ek-In49, L. Eklund59,w, S. Ely68, A. Ene37, E. Epple67, S.
Escher14, J. Eschle50, S. Esen32, T. Evans48, A. Falabella20, J. Fan3, Y.
Fan6, B. Fang73, S. Farry60, D. Fazzini26,j, P. Fedin41, M. Féo48, P.
Fernandez Declara48, A. Fernandez Prieto46, J.M. Fernandez-tenllado Arribas45,
F. Ferrari20,d, L. Ferreira Lopes49, F. Ferreira Rodrigues2, S. Ferreres
Sole32, M. Ferrillo50, M. Ferro-Luzzi48, S. Filippov39, R.A. Fini19, M.
Fiorini21,f, M. Firlej34, K.M. Fischer63, C. Fitzpatrick62, T. Fiutowski34, F.
Fleuret12, M. Fontana13, F. Fontanelli24,h, R. Forty48, V. Franco Lima60, M.
Franco Sevilla66, M. Frank48, E. Franzoso21, G. Frau17, C. Frei48, D.A.
Friday59, J. Fu25, Q. Fuehring15, W. Funk48, E. Gabriel32, T. Gaintseva42, A.
Gallas Torreira46, D. Galli20,d, S. Gambetta58,48, Y. Gan3, M. Gandelman2, P.
Gandini25, Y. Gao5, M. Garau27, L.M. Garcia Martin56, P. Garcia Moreno45, J.
García Pardiñas26,j, B. Garcia Plana46, F.A. Garcia Rosales12, L. Garrido45,
C. Gaspar48, R.E. Geertsema32, D. Gerick17, L.L. Gerken15, E. Gersabeck62, M.
Gersabeck62, T. Gershon56, D. Gerstel10, Ph. Ghez8, V. Gibson55, H.K.
Giemza36, M. Giovannetti23,p, A. Gioventù46, P. Gironella Gironell45, L.
Giubega37, C. Giugliano21,f,48, K. Gizdov58, E.L. Gkougkousis48, V.V.
Gligorov13, C. Göbel70, E. Golobardes84, D. Golubkov41, A. Golutvin61,82, A.
Gomes1,a, S. Gomez Fernandez45, F. Goncalves Abrantes63, M. Goncerz35, G.
Gong3, P. Gorbounov41, I.V. Gorelov40, C. Gotti26, E. Govorkova48, J.P.
Grabowski17, R. Graciani Diaz45, T. Grammatico13, L.A. Granado Cardoso48, E.
Graugés45, E. Graverini49, G. Graziani22, A. Grecu37, L.M. Greeven32, P.
Griffith21,f, L. Grillo62, S. Gromov82, B.R. Gruberg Cazon63, C. Gu3, M.
Guarise21, P. A. Günther17, E. Gushchin39, A. Guth14, Y. Guz44,48, T. Gys48,
T. Hadavizadeh69, G. Haefeli49, C. Haen48, J. Haimberger48, T. Halewood-
leagas60, P.M. Hamilton66, Q. Han7, X. Han17, T.H. Hancock63, S. Hansmann-
Menzemer17, N. Harnew63, T. Harrison60, C. Hasse48, M. Hatch48, J. He6,b, M.
Hecker61, K. Heijhoff32, K. Heinicke15, A.M. Hennequin48, K. Hennessy60, L.
Henry25,47, J. Heuel14, A. Hicheur2, D. Hill49, M. Hilton62, S.E. Hollitt15,
J. Hu17, J. Hu72, W. Hu7, W. Huang6, X. Huang73, W. Hulsbergen32, R.J.
Hunter56, M. Hushchyn81, D. Hutchcroft60, D. Hynds32, P. Ibis15, M. Idzik34,
D. Ilin38, P. Ilten65, A. Inglessi38, A. Ishteev82, K. Ivshin38, R.
Jacobsson48, S. Jakobsen48, E. Jans32, B.K. Jashal47, A. Jawahery66, V.
Jevtic15, M. Jezabek35, F. Jiang3, M. John63, D. Johnson48, C.R. Jones55, T.P.
Jones56, B. Jost48, N. Jurik48, S. Kandybei51, Y. Kang3, M. Karacson48, M.
Karpov81, N. Kazeev81, F. Keizer55,48, M. Kenzie56, T. Ketel33, B. Khanji15,
A. Kharisova83, S. Kholodenko44, K.E. Kim68, T. Kirn14, V.S. Kirsebom49, O.
Kitouni64, S. Klaver32, K. Klimaszewski36, S. Koliiev52, A. Kondybayeva82, A.
Konoplyannikov41, P. Kopciewicz34, R. Kopecna17, P. Koppenburg32, M.
Korolev40, I. Kostiuk32,52, O. Kot52, S. Kotriakhova38,30, P. Kravchenko38, L.
Kravchuk39, R.D. Krawczyk48, M. Kreps56, F. Kress61, S. Kretzschmar14, P.
Krokovny43,v, W. Krupa34, W. Krzemien36, W. Kucewicz35,t, M. Kucharczyk35, V.
Kudryavtsev43,v, H.S. Kuindersma32, G.J. Kunde67, T. Kvaratskheliya41, D.
Lacarrere48, G. Lafferty62, A. Lai27, A. Lampis27, D. Lancierini50, J.J.
Lane62, R. Lane54, G. Lanfranchi23, C. Langenbruch14, J. Langer15, O.
Lantwin50,82, T. Latham56, F. Lazzari29,q, R. Le Gac10, S.H. Lee85, R.
Lefèvre9, A. Leflat40, S. Legotin82, O. Leroy10, T. Lesiak35, B.
Leverington17, H. Li72, L. Li63, P. Li17, Y. Li4, Y. Li4, Z. Li68, X. Liang68,
T. Lin61, R. Lindner48, V. Lisovskyi15, R. Litvinov27, G. Liu72, H. Liu6, S.
Liu4, X. Liu3, A. Loi27, J. Lomba Castro46, I. Longstaff59, J.H. Lopes2, G.H.
Lovell55, Y. Lu4, D. Lucchesi28,l, S. Luchuk39, M. Lucio Martinez32, V.
Lukashenko32, Y. Luo3, A. Lupato62, E. Luppi21,f, O. Lupton56, A. Lusiani29,m,
X. Lyu6, L. Ma4, R. Ma6, S. Maccolini20,d, F. Machefert11, F. Maciuc37, V.
Macko49, P. Mackowiak15, S. Maddrell-Mander54, O. Madejczyk34, L.R. Madhan
Mohan54, O. Maev38, A. Maevskiy81, D. Maisuzenko38, M.W. Majewski34, J.J.
Malczewski35, S. Malde63, B. Malecki48, A. Malinin80, T. Maltsev43,v, H.
Malygina17, G. Manca27,e, G. Mancinelli10, R. Manera Escalero45, D.
Manuzzi20,d, D. Marangotto25,i, J. Maratas9,s, J.F. Marchand8, U. Marconi20,
S. Mariani22,g,48, C. Marin Benito11, M. Marinangeli49, P. Marino49,m, J.
Marks17, P.J. Marshall60, G. Martellotti30, L. Martinazzoli48,j, M.
Martinelli26,j, D. Martinez Santos46, F. Martinez Vidal47, A. Massafferri1, M.
Materok14, R. Matev48, A. Mathad50, Z. Mathe48, V. Matiunin41, C. Matteuzzi26,
K.R. Mattioli85, A. Mauri32, E. Maurice12, J. Mauricio45, M. Mazurek48, M.
McCann61, L. Mcconnell18, T.H. Mcgrath62, A. McNab62, R. McNulty18, J.V.
Mead60, B. Meadows65, C. Meaux10, G. Meier15, N. Meinert76, D. Melnychuk36, S.
Meloni26,j, M. Merk32,79, A. Merli25, L. Meyer Garcia2, M. Mikhasenko48, D.A.
Milanes74, E. Millard56, M. Milovanovic48, M.-N. Minard8, L. Minzoni21,f, S.E.
Mitchell58, B. Mitreska62, D.S. Mitzel48, A. Mödden 15, R.A. Mohammed63, R.D.
Moise61, T. Mombächer15, I.A. Monroy74, S. Monteil9, M. Morandin28, G.
Morello23, M.J. Morello29,m, J. Moron34, A.B. Morris75, A.G. Morris56, R.
Mountain68, H. Mu3, F. Muheim58,48, M. Mukherjee7, M. Mulder48, D. Müller48,
K. Müller50, C.H. Murphy63, D. Murray62, P. Muzzetto27,48, P. Naik54, T.
Nakada49, R. Nandakumar57, T. Nanut49, I. Nasteva2, M. Needham58, I. Neri21,
N. Neri25,i, S. Neubert75, N. Neufeld48, R. Newcombe61, T.D. Nguyen49, C.
Nguyen-Mau49,x, E.M. Niel11, S. Nieswand14, N. Nikitin40, N.S. Nolte48, C.
Nunez85, A. Oblakowska-Mucha34, V. Obraztsov44, D.P. O’Hanlon54, R.
Oldeman27,e, M.E. Olivares68, C.J.G. Onderwater78, A. Ossowska35, J.M. Otalora
Goicochea2, T. Ovsiannikova41, P. Owen50, A. Oyanguren47, B. Pagare56, P.R.
Pais48, T. Pajero29,m,48, A. Palano19, M. Palutan23, Y. Pan62, G. Panshin83,
A. Papanestis57, M. Pappagallo19,c, L.L. Pappalardo21,f, C. Pappenheimer65, W.
Parker66, C. Parkes62, C.J. Parkinson46, B. Passalacqua21, G. Passaleva22, A.
Pastore19, M. Patel61, C. Patrignani20,d, C.J. Pawley79, A. Pearce48, A.
Pellegrino32, M. Pepe Altarelli48, S. Perazzini20, D. Pereima41, P. Perret9,
K. Petridis54, A. Petrolini24,h, A. Petrov80, S. Petrucci58, M. Petruzzo25,
T.T.H. Pham68, A. Philippov42, L. Pica29,n, M. Piccini77, B. Pietrzyk8, G.
Pietrzyk49, M. Pili63, D. Pinci30, F. Pisani48, A. Piucci17, Resmi P.K10, V.
Placinta37, J. Plews53, M. Plo Casasus46, F. Polci13, M. Poli Lener23, M.
Poliakova68, A. Poluektov10, N. Polukhina82,u, I. Polyakov68, E. Polycarpo2,
G.J. Pomery54, S. Ponce48, D. Popov6,48, S. Popov42, S. Poslavskii44, K.
Prasanth35, L. Promberger48, C. Prouve46, V. Pugatch52, H. Pullen63, G.
Punzi29,n, W. Qian6, J. Qin6, R. Quagliani13, B. Quintana8, N.V. Raab18, R.I.
Rabadan Trejo10, B. Rachwal34, J.H. Rademacker54, M. Rama29, M. Ramos
Pernas56, M.S. Rangel2, F. Ratnikov42,81, G. Raven33, M. Reboud8, F. Redi49,
F. Reiss13, C. Remon Alepuz47, Z. Ren3, V. Renaudin63, R. Ribatti29, S.
Ricciardi57, K. Rinnert60, P. Robbe11, A. Robert13, G. Robertson58, A.B.
Rodrigues49, E. Rodrigues60, J.A. Rodriguez Lopez74, A. Rollings63, P.
Roloff48, V. Romanovskiy44, M. Romero Lamas46, A. Romero Vidal46, J.D. Roth85,
M. Rotondo23, M.S. Rudolph68, T. Ruf48, J. Ruiz Vidal47, A. Ryzhikov81, J.
Ryzka34, J.J. Saborido Silva46, N. Sagidova38, N. Sahoo56, B. Saitta27,e, D.
Sanchez Gonzalo45, C. Sanchez Gras32, R. Santacesaria30, C. Santamarina
Rios46, M. Santimaria23, E. Santovetti31,p, D. Saranin82, G. Sarpis59, M.
Sarpis75, A. Sarti30, C. Satriano30,o, A. Satta31, M. Saur15, D. Savrina41,40,
H. Sazak9, L.G. Scantlebury Smead63, S. Schael14, M. Schellenberg15, M.
Schiller59, H. Schindler48, M. Schmelling16, B. Schmidt48, O. Schneider49, A.
Schopper48, M. Schubiger32, S. Schulte49, M.H. Schune11, R. Schwemmer48, B.
Sciascia23, A. Sciubba23, S. Sellam46, A. Semennikov41, M. Senghi Soares33, A.
Sergi24,48, N. Serra50, L. Sestini28, A. Seuthe15, P. Seyfert48, D.M.
Shangase85, M. Shapkin44, I. Shchemerov82, L. Shchutska49, T. Shears60, L.
Shekhtman43,v, Z. Shen5, V. Shevchenko80, E.B. Shields26,j, E. Shmanin82, J.D.
Shupperd68, B.G. Siddi21, R. Silva Coutinho50, G. Simi28, S. Simone19,c, N.
Skidmore62, T. Skwarnicki68, M.W. Slater53, I. Slazyk21,f, J.C. Smallwood63,
J.G. Smeaton55, A. Smetkina41, E. Smith14, M. Smith61, A. Snoch32, M.
Soares20, L. Soares Lavra9, M.D. Sokoloff65, F.J.P. Soler59, A. Solovev38, I.
Solovyev38, F.L. Souza De Almeida2, B. Souza De Paula2, B. Spaan15, E. Spadaro
Norella25,i, P. Spradlin59, F. Stagni48, M. Stahl65, S. Stahl48, P. Stefko49,
O. Steinkamp50,82, S. Stemmle17, O. Stenyakin44, H. Stevens15, S. Stone68,
M.E. Stramaglia49, M. Straticiuc37, D. Strekalina82, F. Suljik63, J. Sun27, L.
Sun73, Y. Sun66, P. Svihra62, P.N. Swallow53, K. Swientek34, A. Szabelski36,
T. Szumlak34, M. Szymanski48, S. Taneja62, F. Teubert48, E. Thomas48, K.A.
Thomson60, M.J. Tilley61, V. Tisserand9, S. T’Jampens8, M. Tobin4, S. Tolk48,
L. Tomassetti21,f, D. Torres Machado1, D.Y. Tou13, M. Traill59, M.T. Tran49,
E. Trifonova82, C. Trippl49, G. Tuci29,n, A. Tully49, N. Tuning32,48, A.
Ukleja36, D.J. Unverzagt17, E. Ursov82, A. Usachov32, A. Ustyuzhanin42,81, U.
Uwer17, A. Vagner83, V. Vagnoni20, A. Valassi48, G. Valenti20, N. Valls
Canudas45, M. van Beuzekom32, M. Van Dijk49, E. van Herwijnen82, C.B. Van
Hulse18, M. van Veghel78, R. Vazquez Gomez46, P. Vazquez Regueiro46, C.
Vázquez Sierra48, S. Vecchi21, J.J. Velthuis54, M. Veltri22,r, A.
Venkateswaran68, M. Veronesi32, M. Vesterinen56, D. Vieira65, M. Vieites
Diaz49, H. Viemann76, X. Vilasis-Cardona84, E. Vilella Figueras60, P.
Vincent13, G. Vitali29, A. Vollhardt50, D. Vom Bruch10, A. Vorobyev38, V.
Vorobyev43,v, N. Voropaev38, R. Waldi76, J. Walsh29, C. Wang17, J. Wang5, J.
Wang4, J. Wang3, J. Wang73, M. Wang3, R. Wang54, Y. Wang7, Z. Wang50, H.M.
Wark60, N.K. Watson53, S.G. Weber13, D. Websdale61, C. Weisser64, B.D.C.
Westhenry54, D.J. White62, M. Whitehead54, D. Wiedner15, G. Wilkinson63, M.
Wilkinson68, I. Williams55, M. Williams64,69, M.R.J. Williams58, F.F.
Wilson57, W. Wislicki36, M. Witek35, L. Witola17, G. Wormser11, S.A. Wotton55,
H. Wu68, K. Wyllie48, Z. Xiang6, D. Xiao7, Y. Xie7, A. Xu5, J. Xu6, L. Xu3, M.
Xu7, Q. Xu6, Z. Xu5, Z. Xu6, D. Yang3, S. Yang6, Y. Yang6, Z. Yang3, Z.
Yang66, Y. Yao68, L.E. Yeomans60, H. Yin7, J. Yu71, X. Yuan68, O.
Yushchenko44, E. Zaffaroni49, K.A. Zarebski53, M. Zavertyaev16,u, M. Zdybal35,
O. Zenaiev48, M. Zeng3, D. Zhang7, L. Zhang3, S. Zhang5, Y. Zhang5, Y.
Zhang63, A. Zhelezov17, Y. Zheng6, X. Zhou6, Y. Zhou6, X. Zhu3, V.
Zhukov14,40, J.B. Zonneveld58, S. Zucchelli20,d, D. Zuliani28, G. Zunica62.

1Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil
2Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil
3Center for High Energy Physics, Tsinghua University, Beijing, China
4Institute Of High Energy Physics (IHEP), Beijing, China
5School of Physics State Key Laboratory of Nuclear Physics and Technology,
Peking University, Beijing, China
6University of Chinese Academy of Sciences, Beijing, China
7Institute of Particle Physics, Central China Normal University, Wuhan, Hubei,
China
8Univ. Savoie Mont Blanc, CNRS, IN2P3-LAPP, Annecy, France
9Université Clermont Auvergne, CNRS/IN2P3, LPC, Clermont-Ferrand, France
10Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
11Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France
12Laboratoire Leprince-Ringuet, CNRS/IN2P3, Ecole Polytechnique, Institut
Polytechnique de Paris, Palaiseau, France
13LPNHE, Sorbonne Université, Paris Diderot Sorbonne Paris Cité, CNRS/IN2P3,
Paris, France
14I. Physikalisches Institut, RWTH Aachen University, Aachen, Germany
15Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
16Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
17Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
18School of Physics, University College Dublin, Dublin, Ireland
19INFN Sezione di Bari, Bari, Italy
20INFN Sezione di Bologna, Bologna, Italy
21INFN Sezione di Ferrara, Ferrara, Italy
22INFN Sezione di Firenze, Firenze, Italy
23INFN Laboratori Nazionali di Frascati, Frascati, Italy
24INFN Sezione di Genova, Genova, Italy
25INFN Sezione di Milano, Milano, Italy
26INFN Sezione di Milano-Bicocca, Milano, Italy
27INFN Sezione di Cagliari, Monserrato, Italy
28Universita degli Studi di Padova, Universita e INFN, Padova, Padova, Italy
29INFN Sezione di Pisa, Pisa, Italy
30INFN Sezione di Roma La Sapienza, Roma, Italy
31INFN Sezione di Roma Tor Vergata, Roma, Italy
32Nikhef National Institute for Subatomic Physics, Amsterdam, Netherlands
33Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, Netherlands
34AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
35Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
36National Center for Nuclear Research (NCBJ), Warsaw, Poland
37Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
38Petersburg Nuclear Physics Institute NRC Kurchatov Institute (PNPI NRC KI),
Gatchina, Russia
39Institute for Nuclear Research of the Russian Academy of Sciences (INR RAS),
Moscow, Russia
40Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
41Institute of Theoretical and Experimental Physics NRC Kurchatov Institute
(ITEP NRC KI), Moscow, Russia
42Yandex School of Data Analysis, Moscow, Russia
43Budker Institute of Nuclear Physics (SB RAS), Novosibirsk, Russia
44Institute for High Energy Physics NRC Kurchatov Institute (IHEP NRC KI),
Protvino, Russia, Protvino, Russia
45ICCUB, Universitat de Barcelona, Barcelona, Spain
46Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de
Santiago de Compostela, Santiago de Compostela, Spain
47Instituto de Fisica Corpuscular, Centro Mixto Universidad de Valencia -
CSIC, Valencia, Spain
48European Organization for Nuclear Research (CERN), Geneva, Switzerland
49Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
Lausanne, Switzerland
50Physik-Institut, Universität Zürich, Zürich, Switzerland
51NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
52Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
53University of Birmingham, Birmingham, United Kingdom
54H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
55Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
56Department of Physics, University of Warwick, Coventry, United Kingdom
57STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
58School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
59School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
60Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
61Imperial College London, London, United Kingdom
62Department of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
63Department of Physics, University of Oxford, Oxford, United Kingdom
64Massachusetts Institute of Technology, Cambridge, MA, United States
65University of Cincinnati, Cincinnati, OH, United States
66University of Maryland, College Park, MD, United States
67Los Alamos National Laboratory (LANL), Los Alamos, United States
68Syracuse University, Syracuse, NY, United States
69School of Physics and Astronomy, Monash University, Melbourne, Australia,
associated to 56
70Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
71Physics and Micro Electronic College, Hunan University, Changsha City,
China, associated to 7
72Guangdong Provencial Key Laboratory of Nuclear Science, Institute of Quantum
Matter, South China Normal University, Guangzhou, China, associated to 3
73School of Physics and Technology, Wuhan University, Wuhan, China, associated
to 3
74Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 13
75Universität Bonn - Helmholtz-Institut für Strahlen und Kernphysik, Bonn,
Germany, associated to 17
76Institut für Physik, Universität Rostock, Rostock, Germany, associated to 17
77INFN Sezione di Perugia, Perugia, Italy, associated to 21
78Van Swinderen Institute, University of Groningen, Groningen, Netherlands,
associated to 32
79Universiteit Maastricht, Maastricht, Netherlands, associated to 32
80National Research Centre Kurchatov Institute, Moscow, Russia, associated to
41
81National Research University Higher School of Economics, Moscow, Russia,
associated to 42
82National University of Science and Technology “MISIS”, Moscow, Russia,
associated to 41
83National Research Tomsk Polytechnic University, Tomsk, Russia, associated to
41
84DS4DS, La Salle, Universitat Ramon Llull, Barcelona, Spain, associated to 45
85University of Michigan, Ann Arbor, United States, associated to 68

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bHangzhou Institute for Advanced Study, UCAS, Hangzhou, China
cUniversità di Bari, Bari, Italy
dUniversità di Bologna, Bologna, Italy
eUniversità di Cagliari, Cagliari, Italy
fUniversità di Ferrara, Ferrara, Italy
gUniversità di Firenze, Firenze, Italy
hUniversità di Genova, Genova, Italy
iUniversità degli Studi di Milano, Milano, Italy
jUniversità di Milano Bicocca, Milano, Italy
kUniversità di Modena e Reggio Emilia, Modena, Italy
lUniversità di Padova, Padova, Italy
mScuola Normale Superiore, Pisa, Italy
nUniversità di Pisa, Pisa, Italy
oUniversità della Basilicata, Potenza, Italy
pUniversità di Roma Tor Vergata, Roma, Italy
qUniversità di Siena, Siena, Italy
rUniversità di Urbino, Urbino, Italy
sMSU - Iligan Institute of Technology (MSU-IIT), Iligan, Philippines
tAGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland
uP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
vNovosibirsk State University, Novosibirsk, Russia
wDepartment of Physics and Astronomy, Uppsala University, Uppsala, Sweden
xHanoi University of Science, Hanoi, Vietnam


# BIOPTIC - a target-agnostic potency-based small molecules search engine

Vlad Vinogradov
Optic Inc.
<EMAIL_ADDRESS>
&Ivan Izmailov
Optic Inc.
<EMAIL_ADDRESS>
&Simon Steshin
Optic Inc.
<EMAIL_ADDRESS>
&Kong T. Nguyen
Optic Inc.
<EMAIL_ADDRESS>

###### Abstract

Recent successes in virtual screening have been made possible by large models
and extensive chemical libraries. However, combining these elements is
challenging: the larger the model, the more expensive it is to run, making
ultra-large libraries unfeasible. To address this, we developed a target-
agnostic, potency-based molecule search model, which allows us to find
structurally dissimilar molecules with similar biological activities. We used
the best practices to design fast retrieval system, based on processor-
optimized SIMD instructions, enabling us to screen the ultra-large 40B Enamine
REAL library with 100% recall rate. We extensively benchmarked our model and
several state-of-the-art models for both speed performance and retrieval
quality of novel molecules.

_K_ eywords Deep Learning $\cdot$ Drug Discovery $\cdot$ Small Molecules
$\cdot$ Virtual Screening $\cdot$ Retrieval System

##  1 Introduction

Drug discovery and development process consists of hit identification, hit to
lead, lead optimization, and preclinical development stages. During this
process, scientists search for ligands that interact with the targets of
interest, and exhibit desired biological activity with favorable
pharmacokinetics and safety profiles. Despite technology improvements and huge
investments in the pharmaceutical industry, the current success rate of
clinical drug candidates remains quite low with only 10% of drug candidates
become marketed drugs[1]. One of the main reasons of drug candidate failures
is the lack of efficacy[2]. How can we significantly improve this process?

One of the modern approaches to drug discovery is virtual screening. This
procedure runs a statistical algorithm on a large number of molecules, acting
as a filter to narrow down the list of potential hits for experimental
testing[3]. There have been several recent successes in virtual screening[4,
5], thanks to the large amounts of data that have become available in the last
decades, which are required by data-greedy modern algorithms. Another key
advantage virtual screening is its ability to scan huge libraries, which is
not feasible in standard high-throughput screening (HTS) limited to hundreds
of thousands of molecules. Advancement of synthetic chemistry methodologies
and automation have led to the generation of huge virtual libraries, such as
Enamine REAL Space[6], consisting of 40 billion different potential molecules,
which is 5 orders of magnitude larger than the size of most HTS screening
decks[7].

Combining powerful data-greedy algorithms and ultra-large screening libraries
poses a significant challenge. It is realistic to apply fast filters (like
Lipinski’s Rule of 5[8], Rule of CNS drugs[9], etc) to huge libraries or to
use neural networks to screen small datasets, but screening huge libraries
using neural networks takes too much time and money. ML-based drug discovery
is a blooming field with a large number of publications. Most of the
publications focus on successful case studies rather than how to scale model
predictions to huge libraries.

Stepping out of drug discovery, the modern world is all about data, and many
industries have already transformed to handle large datasets and are reaping
the benefits. In practice, modern recommendation systems (YouTube, Instagram,
Twitter) and search engines (Google) are heavily bounded by computational
resources, so they have invented new methods to make searches and predictions
extremely fast and meaningful, leading to increased conversions and
revenues[10, 11, 12].

In this work, we apply the best practices of the most established software
engineering tools — recommendation systems and search engines — to scale
virtual screening of ultra-large libraries and demonstrate its performance in
efficacy prediction on BindingDB and its speed on the Enamine REAL Space 40B
library.

##  2 Related works

###  2.1 Deep docking

Docking is a traditional technique in computational chemistry that allows the
simulation of ligand binding to a receptor under certain assumptions. Even
with physically unrealistic assumptions, docking is not feasible for huge
libraries[13], so there are efforts to improve its performance[14].

In the recent papers[15], the authors suggest a hybrid architecture that uses
docking to create a synthetic dataset and then trains a fast neural network on
that dataset to make predictions for the entire library. This approach was
fast enough to screen a 40 billion library.

In our work, we use only ground truth binding data, without relying on docking
simulations. Furthermore, deep docking approach requires a known target
structure and re-training the model for each target, which consumes a lot of
resources. Instead, we employed search engine techniques to develop a global
target-agnostic model that captures ligands’ properties, which are independent
of targets. This allows us to have a single model that can find molecules with
similar activities even when they are structurally dissimilar.

###  2.2 DrugClip

The recent paper DrugClip[16] represents the first attempt to create a
retrieval system for binding affinity prediction. It employs a contrastive
CLIP-like architecture to train a model on protein-ligand pairs, uses the
UniMol 3D encoder, and suggests a new type of augmentation.

While it is a great approach to handle ultra-large libraries, DrugClip
requires structural knowledge of proteins, which is not always available.
Instead, we aimed to create a target-agnostic system capable of finding
molecules with similar activity for all possible targets simultaneously.

###  2.3 Chemprop

Chemprop[17, 18] is a graph neural network that uses the Directed Message
Passing Neural Network (D-MPNN) architecture. What makes this work outstanding
is its extensive practical validation: the model has been successfully used to
discover various antibiotics and senolytics[19, 5, 20]. Since real-life
success is the ultimate benchmark for any predictive model, Chemprop stands
out as one of the few models that have been tested and proven in practical
applications.

However, Chemprop is a supervised learning ‘model trained to predict specific
properties rather than to search for molecules with similar activities, and it
was not developed for integration with huge libraries. Nevertheless, it
remains a solid competitor in terms of accuracy. Therefore, we selected it as
a benchmark for comparison in our study.

##  3 Our contribution

  * •

We developed a single global model for target-agnostic potency-based molecule
search, which allows finding structurally dissimilar molecules with similar
activities.

  * •

We extensively benchmark our model against several state-of-the-art (SOTA)
models in hit identification settings for their ability to find novel,
structurally dissimilar molecules with similar activities.

  * •

We developed an efficient system based on processor-optimized SIMD
instructions and deployed it on a cluster of machines connected to a matrix of
SSD-disks that make it feasible to scan huge 40B libraries super fast.

##  4 Method

In this section, we describe the model’s architecture and the search engine
design. The benchmarks and demonstrations of the system are presented in the
next section.

###  4.1 Modeling

Our model is a SMILES-based transformer [21], which we train in two phases.
The Byte-Pair Encoding (BPE) tokenization [22] is applied to SMILES strings at
each phase. During the first phase, the model is pre-trained in an
unsupervised fashion on a large corpus of unlabeled molecular data from
PubChem and Enamine REAL Space. The training is performed via a masked
language modeling following the procedure from RoBERTa [23] (see fig. 1a).
This phase allows the model to learn the SMILES grammar.

In the second stage, the transformer model is augmented with a multi-pooling
module that aggregates hidden layer output in three different ways: taking
classification token output as is, and taking maximum and average along the
sequence axis of all the remaining tokens. These aggregated vectors are passed
into their corresponding linear layers for projection to a lower dimension
following up with concatenation to form the final embedding. This embedding is
further normalized and passed to the classification layer where the task is to
predict which targets a given molecule is active to. Here, the model is
trained on BindingDB. As the majority of labels are not available for the
BindingDB drug-target pairs, we mask them out to not introduce additional
noise into the training (fig. 1b).

Figure 1:  BIOPTIC’s training method includes two phases. In the initial phase
(a), a transformer neural network is pre-trained in a RoBERTa-style manner
using masked language modeling on SMILES strings of molecules. In the next
phase (b), the transformer model, augmented with a multi-pooling module and
linear projection layers, is fine-tuned to predict which targets a given
molecule is active to, using a masked binary cross-entropy loss function. A
label is 1 if the molecule is active against a certain target and 0 otherwise.
We mask out the unlabeled drug-target pairs.

We use binary cross-entropy loss with masking as an objective for learning
potency-based molecular representations:

$L=\sum_{i=1}^{B}{\sum_{j=1}^{K}{m_{ij}\left[y_{ij}\cdot\log\hat{y}_{ij}+(1-y_{ij})\cdot\log(1-\hat{y}_{ij})\right]}}$
(1)

for a batch of size $B$, where $\hat{y}_{ij}$ is the prediction of the
corresponding molecule $i$ being active to target $j$ (calculated from the
logits of the classification layer after applying sigmoid function), $y_{ij}$
is the corresponding ground truth value (0 or 1), and $K$ is the total number
of targets used in training. $m_{ij}$ represents a masking indicator, which
doesn’t penalize the model on unlabeled data pairs, allowing it to learn from
labeled drug-target interactions only.

###  4.2 Optimization

The chosen model is RoBERTa, with a vocabulary size of 500, 6 hidden layers
with a hidden size of 384, 8 attention heads, and an intermediate size of
1024, while the rest of the parameters are as described in [23]. The model has
a total of 8.7 million parameters. Initially, it is pre-trained in an
unsupervised fashion on a large corpus of data comprising 115 million unique
molecules from PubChem and a random set of 48 million molecules from the
Enamine REAL Space library, totaling over 160 million molecules. The training
is performed using masked language modeling following the procedure from
RoBERTa [23] (fig. 1).

After pre-training, the hidden states of the last layer are extracted from the
CLS token, and the rest of the tokens are max and average pooled. Three
independent linear layers with an output size of 20 are then applied. The
reduced representations are concatenated into a vector of size 60. The
embedding dimensionality is set to 60 as a trade-off between quality and
storage requirements. Embeddings are L2-normalized and passed to the
classification linear layer, where the number of classes depends on the
target-specific data split and averages around 6700 classes. For demonstrating
performance with Hi splits, BIOPTIC is optimized on each data split for a
binary classification objective for 300 epochs with the LARS optimizer [24],
using parameters: $lr=0.1$, weight decay of $1e-3$, $\epsilon=1e-8$ and a
trust coefficient of $1e-3$. To enlarge our dataset, we utilize a
randomization of the SMILES strings and random masking of 15% of tokens in 30%
of cases during BIOPTIC training.

###  4.3 Ultra-fast molecules retrieval system

Figure 2: BIOPTIC molecular retrieval system design. The infrastructure is
configured for processing the Enamine REAL Space 40 billion library. It
includes a) a one-time indexing job performed on a cluster of NVIDIA GPU
machines, and b) a searching job running on a cluster of CPU machines
connected to a matrix of storage disks. The search operation allows a user to
search multi-billion molecular libraries within seconds (see section 5.5).

The neural network produces molecular embeddings as $float16$ value vectors in
a 60-dimensional space. The number of dimensions is tuned to balance quality
and retrieval speed. With these compact representations, we construct an index
of a given molecular library and store it on SSD disks. To perform a search,
we calculate the cosine similarity of a query molecule’s embedding to each of
the indexed molecules’ embeddings in a brute-force fashion. Cosine similarity
is implemented in C++ as matrix multiplication using Single-Instruction-
Multiple-Data (SIMD) instructions. The main advantage of this method is the
100% recall, meaning that if there exists an item in the index similar to the
query, it must be found. This is crucial for large-scale libraries with
billions of molecules.

For the Enamine REAL Space of 40 billion molecules, the inference of the
neural network is performed on an AWS SageMaker cluster of 8 machines, each
with 8 x NVIDIA A100 GPUs in $float16$ precision. Processing also includes
standardized SMILES tokenization using the Tokenizers library [25]. The full
search job comprises the prediction steps and top-k selection steps running on
multiple machines in parallel, followed by a top-k aggregation step performed
using a PostgreSQL database on a single machine (see fig. 2 for the full
architecture overview). Multi-node parallel processing is an essential
component of the design because reading molecular embeddings from SSD disks is
a bottleneck. However, we opt out of using RAM for this task, as it would be
much more expensive when keeping almost 5 TB of embeddings for Enamine REAL
Space 40B in RAM. The search infrastructure is set up as follows: 27 nodes
with 2-core CPUs and 4 GB RAM connected to 270 SSD disks of 21 GB each, with 1
Gbit/s of throughput per disk. With this infrastructure, we are able to cost-
efficiently host clients willing to screen multi-billion molecular libraries
within seconds (see section 5.5 for speed evaluation).

##  5 Results

###  5.1 Hi-benchmark

The ultimate goal of model development is not to excel at synthetic standard
benchmarks retrospectively, but to provide meaningful chemical starting points
for prospective drug discovery projects. We found most benchmarks unsuitable
for this task, as they test models on a chemical space very similar to the
training set, which was shown to the models during the training phase[26].
This is contrary to hit identification stage, as we aim to apply the model in
novel chemical spaces, and we do not find trivial modifications of existing
drugs exciting.

For these reasons, we used Hi-splitter[26] (Hi as in Hit Identification) and
extended the Hi benchmark. This algorithm splits the dataset into training and
holdout sets in such a way that no molecule in the test set has an ECFP4
Tanimoto similarity > 0.40 to any training molecule, representing novel
chemical space (see Appendix B for test set nearest neighbors). To perform
well in these settings, models must generalize, as simple memorization would
be insufficient for good results.

We focus on binding affinity, which is not the focus of the Hi benchmark, so
we created our own benchmark. We selected seven very different targets
representing four protein families (see Appendix C). For each target, we
created a separate dataset with its own training and holdout set. This
approach helped us evaluate models extensively because we aimed to develop a
single model that finds structurally dissimilar molecules with the same
activities across various targets.

###  5.2 Data

We discuss data preprocessing in the Appendix A, but briefly: we took
BindingDB, preprocessed it, and binarized the binding affinity values such
that all Ki, Kd, IC50 less than 10 $\mu$M are considered positive. After that,
we created seven benchmarks from it. For each benchmark, we ran Hi-splitter to
split the dataset into training and holdout sets and further split the holdout
set into test and validation sets in a random fashion. To account for
differences in class balance for different targets (there are many more
actives than inactives for some targets, but not for others), we downsampled
the test and validation datasets to achieve 1:1 activity balance. We also
ensured that no molecule in the training set had an ECFP4 Tanimoto similarity
> 0.40 to the test and validation sets.

###  5.3 Hi-benchmark results

We benchmarked classic and SOTA models for their ability to identify
dissimilar molecules with similar activity. For baselines, we selected three
models: Chemprop, gradient boosting on ECFP4 fingerprints, and plain ECFP4
2048-bit fingerprint similarity search. The first two are popular choices for
molecular property prediction, but they are single-target models and require
retraining for each new target. The last is the standard method for searching
for similar molecules.

To ensure a fair comparison, we conducted a thorough hyperparameter search
using a validation set for each model (see Appendix E) and calculated metrics
on a test set that is dissimilar to the training set. The results are in the
table 1.

For BIOPTIC and fingerprint search, we used 10 most active compounds from the
training set as queries. To get QSAR predictions from the model, we calculated
cosine similarity between embeddings (or fingerprints) of queries and test set
molecules. We ranked test molecules by maximal similarity to the five queries,
with the most similar at the top of the ranking.

Target | Model | ROC AUC | AP | Precision@100 | R-Precision
---|---|---|---|---|---
ACHE | Chemprop | 73.2 | 76.8 | 91.0 | 67.5
GB | 67.1 | 65.2 | 75.0 | 61.8
FP | 61.4 | 57.2 | 59.0 | 60.2
BIOPTIC (last) | 66.2 | 66.0 | 77.0 | 63.3
BIOPTIC (best) | 67.2 | 68.1 | 80.0 | 65.2
AA2AR | Chemprop | 68.5 | 68.7 | 67.0 | 64.4
GB | 57.5 | 54.1 | 59.0 | 55.8
FP | 53.7 | 51.9 | 52.0 | 53.3
BIOPTIC (last) | 72.9 | 73.5 | 79.0 | 65.5
BIOPTIC (best) | 74.4 | 71.7 | 76.9 | 65.9
DRD2 | Chemprop | 59.6 | 55.4 | 53.0 | 54.9
GB | 64.2 | 64.8 | 65.0 | 63.1
FP | 49.5 | 53.5 | 53.5 | 47.2
BIOPTIC (last) | 65.9 | 60.3 | 63.0 | 62.6
BIOPTIC (best) | 69.1 | 63.3 | 68.9 | 65.9
EGFR | Chemprop | 82.0 | 85.2 | 94.0 | 77.3
GB | 73.3 | 77.7 | 92.0 | 67.8
FP | 40.9 | 45.0 | 39.0 | 41.7
BIOPTIC (last) | 82.5 | 84.0 | 95.0 | 75.4
BIOPTIC (best) | 82.4 | 84.1 | 95.9 | 75.4
5-HTT | Chemprop | 83.4 | 82.3 | 74.0 | 74.5
GB | 47.5 | 48.9 | 49.0 | 49.0
FP | 49.8 | 48.9 | 47.0 | 47.1
BIOPTIC (last) | 78.3 | 76.0 | 72.0 | 71.6
BIOPTIC (best) | 80.0 | 78.5 | 72.0 | 71.5
JAK2 | Chemprop | 81.1 | 83.0 | 70.0 | 77.3
GB | 56.5 | 62.1 | 52.0 | 54.7
FP | 51.4 | 49.1 | 51.0 | 47.6
BIOPTIC (last) | 77.8 | 78.9 | 68.0 | 69.0
BIOPTIC (best) | 76.7 | 78.3 | 66.0 | 73.8
KCNH2 | Chemprop | 73.3 | 71.8 | 86.0 | 65.6
GB | 64.0 | 65.7 | 83.0 | 60.7
FP | 54.2 | 53.6 | 53.0 | 54.0
BIOPTIC (last) | 61.4 | 59.2 | 61.0 | 58.2
BIOPTIC (best) | 63.6 | 61.5 | 68.0 | 59.7
Table 1: Benchmark results of the models. Here, checkpoints for Chemprop and
Gradient Boosting models are picked at the best Average Precision metric on
validation while BIOPTIC was tested in two scenarios: last checkpoint at the
end of the training and best checkpoint according to Average Precision on
validation. The best result is in bold. See the metrics description in the
Appendix D.

###  5.4 Effect of queries selection

Figure 3: Impact of query selection strategy on model performance.

In order to investigate the effect of query selection, we evaluated the
BIOPTIC model at the last training checkpoint for each target with different
query selection strategies that include:

  * •

1 most active (according to IC50 value)

  * •

5 most active (according to IC50 value)

  * •

10 most active (according to IC50 value)

  * •

all actives with IC50 $\leq$ 100 nM

  * •

10 diverse active molecules (based on agglomerative clustering into 10
clusters)

The results are in fig. 3. All queries are selected from the training set. The
"1 most active" strategy performs much worse on ACHE and EGFR, potentially due
to noise in the data. The "5 most active," "10 most active," and "10 diverse
active molecules" strategies are generally good.

###  5.5 Retrieval performance on ultra-large molecular libraries

We designed our system with the need to screen ultra-large libraries in mind.
To test the system, we evaluated speed performance under real-life conditions
for screening 40 billion Enamine REAL compounds. We measured three scenarios:

  1. 1.

Scan: a one-time embedding extraction for all the molecules, performed with
GPUs.

  2. 2.

Search: given query molecules and extracted embeddings, find the top-k
molecules with similar activity in the library, performed with CPUs only.

  3. 3.

RAM Search: the same as Search, but embeddings are stored in RAM instead of
slower SSDs. Instead of returning the best top-k, the search returns scores
for the whole library.

The scan operation runs at 30,000 molecules per second on a single NVIDIA A100
GPU and can be straightforwardly parallelized using more GPUs. For conducting
search operations on 40 billion molecules, we deployed a distributed system of
27 cheap nodes with 2-core CPUs and 4 GB RAM (AWS’s c7g.large instances)
connected to 270 SSD disks of 21 GB each, with 1 Gbit/s throughput per disk.

We found this step to be memory-bound by SSD throughput, so we measured
performance for the in-RAM scenario. If all the embeddings are already loaded
in RAM (which might be the case to handle a lot of search queries), we can
process up to 1 billion molecules per second on a single 32-core AWS Graviton3
processor, with the speed increasing with the number of cores according to the
demonstration in (see fig. 4).

See all time measurements in Table 2.

| 1M library | Enamine REAL 6B | Enamine REAL Space 40B
---|---|---|---
Scan |  $500$ msec |  $51$ min |  $5.6$ hours
Search |  $4$ msec |  $20$ sec |  $2$ min $15$ sec
RAM Search |  $1$ msec |  $6.6$ sec |  $43.9$ sec
Table 2: Speed performance of the BIOPTIC search system. Each search is
performed with one query. Figure 4:  RAM search performance grows linearly
with the number of cores.

###  5.6 Scaffold hopping from the box

To demonstrate diverse candidates selection performed with BIOPTIC, we chose
ACHE, AA2AR and EGFR as targets and ran search in Enamine REAL Space 40B
library. We used 10 diverse active compounds as queries and retrieved top-10k
candidates (selection lists). We also took 10k random compounds from Enamine
REAL Space and all known molecules tested for binding to the corresponding
targets from BindingDB. We then mapped Morgan fingerprints of all the
molecules into 2D coordinates with t-SNE. From each selection list, we then
removed all molecules that are closer than 0.4 Tanimoto similarity to any
BindingDB molecule and picked 10 diverse molecules with RDKit’s MaxMinPicker.
The results are shown in 5. It’s seen from the figures that BIOPTIC 1)
prioritizes molecules from a widely span corner of the Enamine chemical space,
2) picks molecules structurally different from what is already known, 3)
selects candidates of different scaffolds.

Figure 5:  Diversity of top candidates retrieved with BIOPTIC from Enamine
REAL Space 40 billion library for a) ACHE, b) AA2AR and c) EGFR.

##  6 Conclusion

With the exponential growth in size of small-molecule screening libraries,
it’s necessary to consider the efficiency of virtual screening models. A 40
billion compound library is already prohibitively large for many models, which
prompted us to develop a ultra-fast retrieval molecular search engine that can
find dissimilar molecules with similar activity. We designed and deployed our
BIOPTIC system using best practices of search engines. As a result of that,
our system can search a pre-indexed library of 40 billion compounds in 40
seconds using CPU-only machines.

We understand that in many cases, chemists are interested in finding
structurally novel molecules instead of minor modifications of known binders.
That’s why we benchmarked our system against several state-of-the-art models
on their ability to find active molecules with ECFP4 Tanimoto similarity < 0.4
to the training set for a wide variety of proteins representing different
families. We believe that our system can search for novel molecules with
desirable biological activity in ultra-large libraries in real-time.

## References

  * [1] Asher Mullard.  Parsing clinical success rates.  Nature Reviews Drug Discovery, 15(7):447–448, 2016.
  * [2] Duxin Sun, Wei Gao, Hongxiang Hu, and Simon Zhou.  Why 90% of clinical drug development fails and how to improve it?  Acta Pharmaceutica Sinica B, 12(7):3049–3062, 2022.
  * [3] Eugene N Muratov, Jürgen Bajorath, Robert P Sheridan, Igor V Tetko, Dmitry Filimonov, Vladimir Poroikov, Tudor I Oprea, Igor I Baskin, Alexandre Varnek, Adrian Roitberg, et al.  Qsar without borders.  Chemical Society Reviews, 49(11):3525–3564, 2020.
  * [4] Ai is a viable alternative to high throughput screening: a 318-target study.  Scientific Reports, 14(1):7526, 2024.
  * [5] Felix Wong, Erica J Zheng, Jacqueline A Valeri, Nina M Donghia, Melis N Anahtar, Satotaka Omori, Alicia Li, Andres Cubillos-Ruiz, Aarti Krishnan, Wengong Jin, et al.  Discovery of a structural class of antibiotics with explainable deep learning.  Nature, 626(7997):177–185, 2024.
  * [6] Oleksandr O Grygorenko, Dmytro S Radchenko, Igor Dziuba, Alexander Chuprina, Kateryna E Gubina, and Yurii S Moroz.  Generating multibillion chemical space of readily accessible screening compounds.  Iscience, 23(11), 2020.
  * [7] Hasan Aldewachi, Radhwan N Al-Zidan, Matthew T Conner, and Mootaz M Salman.  High-throughput screening platforms in the discovery of novel drugs for neurodegenerative diseases.  Bioengineering, 8(2):30, 2021.
  * [8] Christopher A Lipinski, Franco Lombardo, Beryl W Dominy, and Paul J Feeney.  Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings.  Advanced drug delivery reviews, 23(1-3):3–25, 1997.
  * [9] Zoran Rankovic.  Cns drug design: balancing physicochemical properties for optimal brain exposure.  Journal of medicinal chemistry, 58(6):2584–2608, 2015.
  * [10] Paul Covington, Jay Adams, and Emre Sargin.  Deep neural networks for youtube recommendations.  In Proceedings of the 10th ACM conference on recommender systems, pages 191–198, 2016.
  * [11] Scaling the instagram explore recommendations system — engineering.fb.com.  https://engineering.fb.com/2023/08/09/ml-applications/scaling-instagram-explore-recommendations-system/.  [Accessed 12-06-2024].
  * [12] Sergey Brin and Lawrence Page.  The anatomy of a large-scale hypertextual web search engine.  Computer Networks, 30:107–117, 1998.
  * [13] Thomas Gaillard.  Evaluation of autodock and autodock vina on the casf-2013 benchmark.  Journal of chemical information and modeling, 58(8):1697–1706, 2018.
  * [14] Francesco Gentile, Vibudh Agrawal, Michael Hsing, Anh-Tien Ton, Fuqiang Ban, Ulf Norinder, Martin E Gleave, and Artem Cherkasov.  Deep docking: a deep learning platform for augmentation of structure based drug discovery.  ACS central science, 6(6):939–949, 2020.
  * [15] Francesco Gentile, Michael Fernandez, Fuqiang Ban, Anh-Tien Ton, Hazem Mslati, Carl F Perez, Eric Leblanc, Jean Charle Yaacoub, James Gleave, Abraham Stern, et al.  Automated discovery of noncovalent inhibitors of sars-cov-2 main protease by consensus deep docking of 40 billion small molecules.  Chemical science, 12(48):15960–15974, 2021.
  * [16] Bowen Gao, Bo Qiang, Haichuan Tan, Yinjun Jia, Minsi Ren, Minsi Lu, Jingjing Liu, Wei-Ying Ma, and Yanyan Lan.  Drugclip: Contrasive protein-molecule representation learning for virtual screening.  Advances in Neural Information Processing Systems, 36, 2024.
  * [17] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  Journal of chemical information and modeling, 59(8):3370–3388, 2019.
  * [18] Esther Heid, Kevin P Greenman, Yunsie Chung, Shih-Cheng Li, David E Graff, Florence H Vermeire, Haoyang Wu, William H Green, and Charles J McGill.  Chemprop: A machine learning package for chemical property prediction.  Journal of Chemical Information and Modeling, 64(1):9–17, 2023.
  * [19] Jonathan M Stokes, Kevin Yang, Kyle Swanson, Wengong Jin, Andres Cubillos-Ruiz, Nina M Donghia, Craig R MacNair, Shawn French, Lindsey A Carfrae, Zohar Bloom-Ackermann, et al.  A deep learning approach to antibiotic discovery.  Cell, 180(4):688–702, 2020.
  * [20] Felix Wong, Satotaka Omori, Nina M Donghia, Erica J Zheng, and James J Collins.  Discovering small-molecule senolytics with deep neural networks.  Nature Aging, 3(6):734–750, 2023.
  * [21] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [22] Rico Sennrich, Barry Haddow, and Alexandra Birch.  Neural machine translation of rare words with subword units, 2016.
  * [23] Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized bert pretraining approach, 2019.
  * [24] Yang You, Igor Gitman, and Boris Ginsburg.  Large batch training of convolutional networks.  arXiv preprint arXiv:1708.03888, 2017.
  * [25] HuggingFace.  Tokenizers, 2024.  https://huggingface.co/docs/tokenizers.
  * [26] Simon Steshin.  Lo-hi: Practical ml drug discovery benchmark.  Advances in Neural Information Processing Systems, 36, 2024.
  * [27] Tuomo Kalliokoski, Christian Kramer, Anna Vulpetti, and Peter Gedeck.  Comparability of mixed ic50 data–a statistical analysis.  PloS one, 8(4):e61007, 2013.
  * [28] A Patrícia Bento, Anne Hersey, Eloy Félix, Greg Landrum, Anna Gaulton, Francis Atkinson, Louisa J Bellis, Marleen De Veij, and Andrew R Leach.  An open source chemical structure curation pipeline using rdkit.  Journal of Cheminformatics, 12:1–16, 2020.

##  Appendix A Data preprocessing

We constructed benchmarks using the database BindingDB_All_202310. This
section discusses the preprocessing steps.

We removed approximately 40 SMILES entries containing the asterisk (*). After
that, we binarized the activity values for Ki (nM), Kd (nM), and IC50 (nM). We
began by dividing all the IC50 values by a coefficient of 2.3 following [27].
Both Ki and IC50 values contain important information useful for training, but
they are not directly comparable. We used the 2.3 coefficient, as it is the
common difference between IC50 and Ki values in databases; we hypothesized
this could improve training. Some values contained greater (">") and lesser
("<") signs, so we considered a molecule active if its activity was less than
10,000 nM, and inactive if its activity was greater than 10,000 nM. All
ambiguous values (such as > 5,000 nM) were discarded. If a SMILES-Target pair
had several activity values, we preferred one of the binary values in this
order: Ki, Kd, IC50.

As targets, we selected fields UniProt (SwissProt) Primary ID of Target Chain
and UniProt (TrEMBL) Primary ID of Target Chain. We removed all data points
that had either both or none of these fields. After that, we standardized the
SMILES with ChEMBL’s molecular standardization pipeline [28] and removed data
points with the same SMILES and target, but conflicting binarized activity
values. Finally, we eliminated PAINS alerts using the medchem library.

##  Appendix B Hi-split nearest neighbors

Train molecule 1Train molecule 2Test neighbor of 1Test neighbor of 2 Figure 6:
The most similar test neighbors of the train molecules in the DRD2 dataset.

##  Appendix C Datasets analysis

To create a dataset for a specific target, we selected all molecules with
known activity to the target and split them into train and holdout sets using
Hi-splitter. We further randomly split the holdout set into test and
validation datasets. After that, we removed from the train set all molecules
that have an ECFP4 Tanimoto similarity > 0.4 to those in the test and
validation datasets. Next, we downsampled the test and validation datasets to
achieve a 1:1 class balance. The detailed statistics are in the Table 3.

Name |  Shorthand | UniProt ID |  Family | Type | Actives | Inactives
---|---|---|---|---|---|---
Acetylcholinesterase |  ACHE | P22303 |  Carboxylesterase | Globular | 3396 | 2535
Adenosine A2A receptor |  AA2AR | P29274 |  G protein-coupled receptor (GPCR) | Membrane | 6165 | 1002
D2 dopamine receptor |  DRD2 | P14416 |  G protein-coupled receptor (GPCR) | Membrane | 7912 | 1243
Epidermal growth factor receptor |  EGFR | P00533 |  Receptor tyrosine kinase | Membrane | 8130 | 2715
Serotonin transporter |  5-HTT | P31645 |  Solute carrier family | Membrane | 4452 | 783
Tyrosine-protein kinase JAK2 |  JAK2 | O60674 |  Non-receptor tyrosine kinase | Globular | 11986 | 1479
Potassium voltage-gated channel subfamily H member 2 |  KCNH2 | Q12809 |  Voltage-gated potassium channel | Membrane | 3644 | 5560
Table 3: Targets from BindingDB selected for evaluation. The selection
strategy is based on including different families and types of proteins while
maintaining enough active and inactive compounds. This table summarizes the
total number of active and inactive compounds present in the BindingDB after
all the filtering steps.

##  Appendix D Metrics description

Virtual screening is used as a filter, to select the top rated molecules. The
primary goal is to ensure that active molecules are ranked at the top,
optimizing the screening efficiency. Here, we discuss several key metrics used
for this purpose:

###  D.1 ROC AUC

The Receiver Operating Characteristic (ROC) Area Under the Curve (AUC) is a
standard ranking metric. It represents the area under the curve plotting the
True Positive Rate (TPR) against the False Positive Rate (FPR). An AUC of 0.5
indicates random ranking, while an AUC of 1.0 signifies perfect ranking. We
used the sklearn.metrics.roc_auc_score to compute this metric.

###  D.2 AP

Average Precision (AP) score summarizes a precision-recall curve by
calculating the weighted average of precision values at each threshold. The
weights are determined by the increase in recall from the previous threshold.
This metric provides a single-figure measure of quality across different
recall levels and prefer models with early recognition. We used the
sklearn.metrics.average_precision_score to compute this metric.

###  D.3 Precision@100

Precision@100 measures the fraction of active molecules within the top 100
molecules in the ranked list. This metric is particularly useful for assessing
the effectiveness of virtual screening when focusing on the top candidates.

###  D.4 R-Precision

R-Precision is equivalent to Precision@R, where R is the number of total
active molecules in the test set. This metric assesses the proportion of
active molecules retrieved when considering a number of top-ranked molecules
equal to the total number of actives. It helps in understanding the retrieval
performance relative to the available actives.

##  Appendix E Hyperparameter tuning

To ensure a fair comparison, we conducted hyperparameter tuning.

###  E.1 Chemprop

We used Chemprop 2.0.0 and performed 80 iterations of random search, sampling
these parameters:

    PARAMS_TO_SAMPLE = {
        "max_epochs": [60],
        "mp": ["bond", "atom"],
        "d_h": [150, 300, 600, 900, 1200, 2400],
        "depth": [3, 4, 5, 6],
        "dropout": [0.0, 0.0, 0.2, 0.3, 0.4, 0.5],
        "batch_norm": [False, True],
        "agg": ["mean", "sum", "norm"],
        "bias": [True, False],
        "init_lr": [1e-4],
        "max_lr": [1e-3, 3e-4, 1e-4],
        "final_lr": [1e-4],
        "batch_size": [32, 64, 96],
        "ffn_hidden_dim": [300, 600, 900],
        "ffn_n_layers": [1, 2, 3],
        "ffn_dropout": [0.0, 0.0, 0.1, 0.2, 0.3],
    }

The best checkpoint was selected based on the best validation average
precision.

###  E.2 GB

We used GradientBoostingClassifier from scikit-learn 1.2.2. We trained it on
ECFP4 2048 fingerprints. For each benchmark, we ran 30 iterations of random
search, sampling these parameters:

    params = {
        ’n_estimators’: [10, 50, 100, 150, 200, 250, 500],
        ’learning_rate’: [0.01, 0.1, 0.3, 0.5, 0.7, 1.0],
        ’subsample’: [0.4, 0.7, 0.9, 1.0],
        ’min_samples_split’: [2, 3, 5, 7],
        ’min_samples_leaf’: [1, 3, 5],
        ’max_depth’: [2, 3, 4],
        ’max_features’: [None, ’sqrt’]
    }

The best checkpoint was selected based on the best validation average
precision.