# Forward $J/\psi+J/\psi$ and $J/\psi+\psi^{\prime}$ production with High
Energy Factorization

S.P. Baranov1, A.V. Lipatov2,3, M.A. Malyshev2,4,
A.A. Prokhorov3, P.M. Zhang5

1P.N. Lebedev Institute of Physics, Moscow 119991, Russia
2Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow State University,
Moscow 119991, Russia
3Joint Institute for Nuclear Research, Dubna 141980, Moscow region, Russia
4Moscow Aviation Institute, Moscow 125993, Russia
5School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519082, China

Abstract

We calculate the cross sections of associated $J/\psi+\psi^{\prime}$ and
$J/\psi+J/\psi$ production in $pp$ collisions at $\sqrt{s}=13$ TeV in the
forward kinematic region. The High Energy Factorization
($k_{T}$-factorization) framework supplemented with the Catani-Ciafaloni-
Fiorani-Marchesini evolution of gluon densities in a proton is applied. We
demonstrate that latest data on $J/\psi+J/\psi$ production and first
experimental data on $J/\psi+\psi^{\prime}$ events taken very recently by the
LHCb Collaboration can be described well by the color singlet terms and
contributions from the double parton scattering (DPS) with the standard choice
for $\sigma_{\rm eff}$ parameter. The relative production rate
$\sigma(J/\psi+\psi^{\prime})/\sigma(J/\psi+J/\psi)$ is found to be sensitive
to the DPS terms as well as to feeddown contributions.

Keywords: heavy quarkonia, high-energy factorization, CCFM evolution, TMD
gluon densities, double parton scattering

Inclusive and associated production of quarkonium states in hadron-hadron
collisions at high energies attracts much attention from both theoretical and
experimental sides. It probes Quantum Chromodynamics (QCD) in perturbative and
non-perturbative regimes and provides important information about the
interaction dynamics. New data, as soon as they appear, trig theoretical
activity aimed at their understanding and description. Our present note is
motivated by the first measurement[1] of the associated $J/\psi+\psi^{\prime}$
production and very recent data[2] on $J/\psi+J/\psi$ production in $pp$
collisions taken by the LHCb Collaboration in the forward kinematic region at
$\sqrt{s}=13$ TeV. These data include contributions from both single and
double parton scattering production mechanisms (SPS and DPS, respectively).
The data on $J/\psi+\psi^{\prime}$ events present the differential cross
sections as functions of several kinematic variables, namely, the rapidity and
azimuthal angle differences between the $J/\psi$ and $\psi^{\prime}$ mesons
$\Delta y(J/\psi,\psi^{\prime})$ and $\Delta\phi(J/\psi,\psi^{\prime})$, the
transverse momentum $p_{T}(J/\psi,\psi^{\prime})$, rapidity
$y(J/\psi,\psi^{\prime})$ and invariant mass $M(J/\psi,\psi^{\prime})$ of the
$J/\psi+\psi^{\prime}$ system. Similar observables as well as transverse
momentum, rapidity of either $J/\psi$ meson and transverse momentum asymmetry
${\cal A}_{T}$ of the two $J/\psi$ mesons have been investigated in the latest
$J/\psi+J/\psi$ analysis[2]. Moreover, the relative production rate
$\sigma(J/\psi+\psi^{\prime})/\sigma(J/\psi+J/\psi)$ has been reported[1] for
the first time.

A commonly accepted framework for the description of heavy quarkonia
production and decays is provided by non-relativistic QCD (NRQCD)
approximation[3, 4]. Explicit calculations[5] show that the main role in the
production of quarkonium pairs at forward rapidities is played by the color
singlet (CS) mechanism, whereas the color octet (CO) contributions are much
smaller. The tree-level next-to-leading order (NLO∗) CS calculations for SPS
mechanism performed with the helac-onia tool[6, 7] tend to overestimate the
LHCb data[1], especially at low $\Delta y(J/\psi,\psi^{\prime})$, low
$M(J/\psi,\psi^{\prime})$ and moderate $p_{T}(J/\psi,\psi^{\prime})$, though
still remain consistent with them within the large theoretical uncertainties.
The latest data[2] on forward $J/\psi+J/\psi$ events can also be described by
the NLO∗ CS calculations within the large uncertainties, although the DPS
contributions were not subtracted from the measurements111The leading-order
(LO) NRQCD predictions $J/\psi+J/\psi$ production are known[8] and the NLO∗
contributions to the CS and CO mechanisms have been calculated[9]..

In the present note we analyse the LHCb data[1, 2] on double charmonia
production in the framework of the $k_{T}$-factorization[10] or, equivalently,
the High Energy Factorization[11] formalism. This formalism has certain
technical advantages in the ease of including higher-order radiative
corrections that can be taken into account in the form of transverse momentum
dependent (TMD) gluon distributions in a proton222See, for example, review[12]
for more information.. Recently, it has been successfully applied to the
double $J/\psi$ production at the LHC[13, 14, 15]. In particular, it was
demonstrated[14, 15] that the early data on $J/\psi+J/\psi$ events collected
at $\sqrt{s}=7$ and $13$ TeV at forward rapidities can be described well by
the sum of color-singlet SPS and DPS mechanisms. It was also confirmed that
the CO contributions at the LHC conditions can be safely omitted in the region
of small invariant masses $M(J/\psi,J/\psi)$ in the forward rapidity region.
Now, in addition to $J/\psi$ pairs, we consider $J/\psi+\psi^{\prime}$
production. Our study is even more stimulated by the fact that previous
theoretical attempts[9, 6, 7, 16] were not fairly successful.

Figure 1: Examples of Feynman diagrams, contributing to the CS production of
charmonium pairs.

We preserve full consistency with our previous investigations[13, 14, 15] and
employ the $k_{T}$-factorization QCD approach[10, 11], as it was mentioned
above. This approach is based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL)[17]
or Catani-Ciafaloni-Fiorani-Marchesini (CCFM)[18] gluon evolution equations
and can be used as a convenient alternative to explicit higher-order pQCD
calculations.

In the NRQCD framework, the colliding gluons produce heavy quark-antiquark
pairs which further evolve into real mesons. At forward rapidities, the
leading CS contribution to associated $J/\psi+J/\psi$ or
$J/\psi+\psi^{\prime}$ production is represented by an
$\mathcal{O}(\alpha^{4}_{s})$ partonic subprocess

$\displaystyle g^{*}(k_{1})+g^{*}(k_{2})\to
c\bar{c}\left[{}^{3}S_{1}^{(1)}\right](p_{1})+c\bar{c}\left[{}^{3}S_{1}^{(1)}\right](p_{2}),$
(1)

where the four-momenta of all particles are indicated in the parentheses. This
subprocess includes also feeddown contributions (namely,
$J/\psi+\psi^{\prime}$ and/or $\psi^{\prime}+\psi^{\prime}$ production), where
excited charmonium state $\psi^{\prime}$ decays into $J/\psi$ meson. Other
subleading subprocesses which involve the production of intermediate $P$-wave
states and/or additional gluons (see[15] for more details) were found to be
small and not taken into account here. In accordance with the
$k_{T}$-factorization, the initial gluons have nonzero transverse momenta
$k_{1}^{2}\equiv-{\mathbf{k}}_{1T}^{2}\neq 0$,
$k_{2}^{2}\equiv-{\mathbf{k}}_{2T}^{2}\neq 0$ and are off-shell. Some examples
of the relevant Feynman diagrams are shown in Fig. 1.

The gauge invariant expression for off-shell (depending on the virtualities of
incoming gluons) production amplitude (1) has been calculated earlier[13]. It
contains spin and color projection operators[19, 20, 21, 22, 23] which
guarantee the proper quantum numbers of the final state charmonia. In the
NRQCD approximation which we are using the meson mass must be strictly equal
to the sum of the constituent quark masses. However, this cannot be fulfilled
simultaneously for both $J/\psi$ and $\psi^{\prime}$. So that, for the
associated $J/\psi+\psi^{\prime}$ production we used a compromise value
$m_{c}=(m_{J/\psi}+m_{\psi^{\prime}})/4$ (see discussion below). In all other
respects our present theoretical scheme is identical to that used in previous
studies[13, 14, 15]. We only note that the initial gluon spin density matrix
is taken in the form of so called ”non-sense gauge”, namely,
$\sum\epsilon^{\mu}\epsilon^{*\nu}=\mathbf{k}^{\mu}_{T}\mathbf{k}^{\nu}_{T}/\mathbf{k}_{T}^{2}$.
This expression converges to the ordinary $-g^{\mu\nu}/2$ in the collinear
limit $\mathbf{k}_{T}\to 0$ after averaging over the azimuthal angle.

According to the $k_{T}$-factorization prescription[10, 11], the contribution
from CS production mechanism (1) to $J/\psi+\psi^{\prime}$ cross section is
calculated as a convolution of the corresponding off-mass shell production
amplitudes and TMD gluon densities in a proton:

$\displaystyle\sigma(pp\to J/\psi+\psi^{\prime}+X)=\int{1\over
16\pi(x_{1}x_{2}s)^{2}}|{\cal\bar{A}}(g^{*}+g^{*}\to
J/\psi+\psi^{\prime})|^{2}\times$ $\displaystyle\times
f_{g}(x_{1},\mathbf{k}_{1T}^{2},\mu^{2})f_{g}(x_{2},\mathbf{k}_{2T}^{2},\mu^{2})d{\mathbf{k}}_{1T}^{2}d{\mathbf{k}}_{2T}^{2}d{\mathbf{p}}_{1T}^{2}dy_{1}dy_{2}{d\phi_{1}\over
2\pi}{d\phi_{2}\over 2\pi}{d\psi_{1}\over 2\pi},$ (2)

where $\psi_{1}$ is the azimuthal angle of the outgoing $J/\psi$ meson,
$\phi_{1}$ and $\phi_{2}$ are the azimuthal angles of the initial off-shell
gluons having the longitudinal momentum fractions $x_{1}$ and $x_{2}$, $y_{1}$
and $y_{2}$ are the center of mass rapidities of the produced particles. Note
that the cross section of associated $J/\psi+J/\psi$ production can be
calculated in a similar way. Here $f_{g}(x,{\mathbf{k}}_{T}^{2},\mu^{2})$ is
the TMD gluon density in a proton taken at the scale $\mu^{2}$. For the
latter, we have tried two recent sets, referred to as JH’2013 set 2[24] and
LLM’2022[25]. These gluon densities have been obtained from a numerical
solution of the CCFM equation and are now widely used in phenomenological
applications (see, for example,[26, 27, 28, 29, 30] and references therein).
The parameters of (rather empirical) input distributions employed in the
JH’2013 gluon were derived from a fit to the HERA data on the proton structure
functions $F_{2}(x,Q^{2})$ and $F_{2}^{c}(x,Q^{2})$ at small $x$. An
analytical expression for the input gluon density in the very recent LLM’2022
set was suited to the best description of the LHC data on the charged hadron
production at low transverse momenta in the framework of modified soft quark-
gluon string model[31, 32] with taking into account the gluon saturation
effects important at low scales. The necessary phenomenological parameters
were deduced from the LHC and HERA data on several hard QCD processes (see[25]
for more information). Both the JH’2013 set 2 and LLM’2022 gluon distributions
are available from the popular tmdlib package[33] and Monte-Carlo generator
pegasus[34].

Figure 2: Differential cross sections of forward $J/\psi+J/\psi$ production at
$\sqrt{s}=13$ TeV as functions of their rapidity difference $\Delta
y(J/\psi,J/\psi)$, azimuthal angle difference $\Delta\phi(J/\psi,J/\psi)$,
transverse momentum assymetry ${\cal A}_{T}$, transverse momentum
$p_{T}(J/\psi,J/\psi)$, rapidity $y(J/\psi,J/\psi)$, invariant mass
$M(J/\psi,J/\psi)$ and transverse momentum $p_{T}(J/\psi)$ and rapidity
$y(J/\psi)$ of either of $J/\psi$ meson. Separately shown contrbutions from
the DPS production mechanism. The data are of LHCb[2].

In addition to the above, we also take into account the contributions from DPS
mechanism, which is now widely discussed in the literature. In contrast with
SPS subprocess (1), where every quarkonium pair is produced in a single gluon-
gluon collision, the DPS events originate from two independent parton
interactions and are expected to be important at forward rapidities (see, for
example,[35, 36, 37]). We apply a commonly used factorization formula:

$\displaystyle\sigma_{\rm DPS}(pp\to{\cal Q}_{1}+{\cal Q}_{2}+X)={m\over
2}{\sigma(pp\to{\cal Q}_{1}+X)\sigma(pp\to{\cal Q}_{2}+X)\over\sigma_{\rm
eff}},$ (3)

where ${\cal Q}_{i}=J/\psi$, $\psi^{\prime}$ or $\chi_{cJ}$ and $\sigma_{\rm
eff}$ is the effective cross section which represents the degree of overlap in
the transverse space between two partonic interactions that constitute the DPS
process. Note that $m=1$ if ${\cal Q}_{1}={\cal Q}_{2}$ and $m=2$ if ${\cal
Q}_{1}$ and ${\cal Q}_{2}$ are different, thus preventing double counting for
identical particles. The effective DPS cross section can be considered as a
normalization constant which incorporates all ”DPS unknowns” into a single
phenomenological parameter. The derivation of formula (3) relies on two
approximations: first, the double parton distribution function can be
decomposed into longitudinal and transverse components and, second, the
longitudinal component reduces to the diagonal product of two independent
single parton densities. The latter is generally acceptable for the conditions
of LHCb experiment. The inclusive cross sections $\sigma(pp\to J/\psi+X)$
and/or $\sigma(pp\to\psi^{\prime}+X)$ involved in (3) are calculated within
the NRQCD framework supplemented with $k_{T}$-factorization (see, for
example,[38] and references therein). The effective cross section is chosen as
$\sigma_{\rm eff}=15$ mb, which is a commonly accepted value333A very close
value $\sigma_{\rm eff}=13.8$ mb has been extracted[14] from a fit to the LHCb
data on the double $J/\psi$ production. However, some other works on the
double quarkonia production may seem to favor much smaller values:
$\sigma_{\rm eff}\sim 5$ mb[39, 40, 41, 42]. We argue that this could be due
to an incomplete treatment of the SPS contributions, see[14] for more
details.. Our calculation incorporates all possible feeddown contributions to
the DPS cross section except the $\chi_{c0}$ production due to its low
branching fraction to the $J/\psi$.

The meson masses were taken in our calculations as $m(J/\psi)=3.096$ GeV,
$m(\chi_{c1})=3.511$ GeV, $m(\chi_{c2})=3.556$ GeV and
$m(\psi^{\prime})=3.686$ GeV (see[43]). A complete list of used CS and CO
LDMEs can be found[15]. The branching fractions are ${Br}(\chi_{c1}\to
J/\psi+\gamma)=0.339$, ${Br}(\chi_{c2}\to J/\psi+\gamma)=0.192$ and
${Br}(\psi^{\prime}\to J/\psi+X)=0.529$. The renormalization and factorization
scales were calculated as
$\mu_{R}^{2}=\mu_{F}^{2}=\xi^{2}(\hat{s}+\mathbf{Q}_{T}^{2})$, where
$\mathbf{Q}_{T}$ is the total transverse momentum of the initial off-shell
gluon pair. This choice observes full consistency with our previous work[14]
and is dictated by the CCFM evolution. It comes from the kinematics of gluon
radiation, namely from the gluon angular ordering condition. The quantity
$\hat{s}+\mathbf{Q}_{T}^{2}$ describes the phase space available for the
(angular ordered) gluon radiation and sets the upper bound both for the
emitted and exchanged transverse momentum (see[18] for more information). We
use a two-loop formula for QCD strong coupling with $N_{f}=4$ active quark
flavours and $\Lambda_{\rm QCD}=200$ MeV. To estimate the theoretical
uncertainties of our calculations, we varied the parameter $\xi$ between $1/2$
and $2$.

Our simulations closely follow the experimental setup[1, 2]. In particular,
the produced $J/\psi$ and $\psi^{\prime}$ mesons are required to be in the
region of $2<y<4.5$ and $p_{T}<14$ GeV, where $y$ is the rapidity of a final
state particle. All these conditions have been implemented in our numerical
program.

Figure 3: Differential cross sections of associated $J/\psi$ and
$\psi^{\prime}$ production as functions of their rapidity difference $\Delta
y(J/\psi,\psi^{\prime})$, azimuthal angle difference
$\Delta\phi(J/\psi,\psi^{\prime})$, transverse momentum
$p_{T}(J/\psi,\psi^{\prime})$, rapidity $y(J/\psi,\psi^{\prime})$ and
invariant mass $M(J/\psi,\psi^{\prime})$. Separately shown are the
contrbutions from DPS production mechanism. The data are of LHCb[1].

The results of our calculations of the differential cross sections for forward
$J/\psi+J/\psi$ and $J/\psi+\psi^{\prime}$ production are shown in Figs. 2 and
3, respectively. The theoretical uncertainty bands (represented by the shaded
regions) are shown for JH’2013 set 2 gluon density and related with scale
uncertainties, which have been estimated in a usual way, by varying the $\xi$
parameter around its default value $\xi=1$ by a factor of $2$. This was
accompanied with using the JH’2013 set $2+$ and JH’2013 set $2-$ gluon
densities in place of default distribution, in accordance with[24]. We also
include the mass uncertainties for CS $J/\psi+\psi^{\prime}$ subprocess by
varying the c-quark mass between $m(J/\psi)/2<m_{c}<m(\psi^{\prime})/2$. These
effects are added in quadratures with scale uncertainties and represented by
the light shaded bands in Fig. 3. To highlight the role of DPS mechanism, we
separately show the corresponding contributions. One can see that our
predictions agree well with the latest LHCb data. The only exception is seen
at larger invariant masses $M(J/\psi,J/\psi)$ and $M(J/\psi,\psi^{\prime})$,
where the calculated cross sections systematically underestimate the
experimental results. Note that still additional contributions related to
multiple gluon radiation in the initial gluon evolution cascade444Such terms
are extremely important at central rapidities[14, 15]. could play a role at
large $M(J/\psi,J/\psi)$ and $M(J/\psi,\psi^{\prime})$. However, an accurate
account of all these terms needs rather lengthy numerical calculations and is
therefore left out of our present scope.

We find special interest in the correlation observables
$\Delta\phi(J/\psi,\psi^{\prime})$, $\Delta\phi(J/\psi,J/\psi)$ and in the
$p_{T}$ balance observables $p_{T}(J/\psi,J/\psi)$ and
$p_{T}(J/\psi,\psi^{\prime})$ as they are known to be particularly sensitive
to the non-collinear gluon evolution dynamics in a proton (see, for
example,[44] and references therein). We point out a good description (within
the uncertainties) of all these distributions with the JH’2013 set 2 gluon
density. The LLM-based predictions, in general, tend to slightly underestimate
the LHCb data[1, 2] and JH’2013 set 2 predictions. In fact, the measured
fiducial cross sections of forward $J/\psi+J/\psi$ and $J/\psi+\psi^{\prime}$
production are $\sigma(J/\psi+J/\psi)=16.36\pm 0.92$ nb and
$\sigma(J/\psi+\psi^{\prime})=4.5\pm 0.8$ nb. They can be compared with the
JH’2013 set 2 and LLM’2022 results, which are $\sigma^{\rm
JH}(J/\psi+J/\psi)=16.2^{+7.2}_{-3.5}$ (scale unc.) nb, $\sigma^{\rm
JH}(J/\psi+\psi^{\prime})=6.1^{+3.7}_{-1.5}$ (scale unc.)${}^{+6.2}_{-2.6}$
(mass unc.) nb and $\sigma^{\rm LLM}(J/\psi+J/\psi)=12.1^{+4.9}_{-2.3}$ (scale
unc.) nb, $\sigma^{\rm LLM}(J/\psi+\psi^{\prime})=4.5^{+2.5}_{-0.9}$ (scale
unc.)${}^{+4.3}_{-1.8}$ (mass unc.) nb, respectively. Moreover, one can see a
clear difference between the JH’2013 set 2 and LLM’2022 predictions in the
shape of azimuthal angle correlations and transverse momenta correlations of
the produced mesons. It demonstrates again that such observables could be
useful to discriminate the different TMD evolution scenarios.

The estimated contributions from the DPS production mechanism are typically
small and only become important in some specific regions like low
$p_{T}(J/\psi,J/\psi)$, $p_{T}(J/\psi,\psi^{\prime})$ and large
$M(J/\psi,J/\psi)\geq 15$ GeV, $M(J/\psi,\psi^{\prime})\geq 15$ GeV, $|\Delta
y(J/\psi,J/\psi)|\geq 1.5$, $|\Delta y(J/\psi,\psi^{\prime})|\geq 1.5$. Our
calculations clearly show that taking these terms into account improves the
quality of the data description and becomes necessary in the forward kinematic
region. The role of DPS and feeddown contributions is exhibited in Fig. 4,
where the ratio $\sigma(J/\psi+\psi^{\prime})/\sigma(J/\psi+J/\psi)$ is shown
as a function of the different kinematic variables. We find that the relative
production rate is very sensitive to the DPS and feeddown yields555Note that
only the SPS terms are affected by the choice of TMD gluon density in a
proton. The DPS predictions are stable, because when we switch to a different
set of TMD gluons we have to accordingly change the LDME values, thus making
the net result the same.. Finally, we can conclude that the first LHCb
measurement[1] of the associated $J/\psi+\psi^{\prime}$ production and the
very recent experimental data[2] on $J/\psi+J/\psi$ production at
$\sqrt{s}=13$ TeV can be well described by the color singlet SPS mechanism and
DPS contributions calculated with the standard choice of $\sigma_{\rm eff}$
parameter.

Figure 4: Relative production rate
$\sigma(J/\psi+\psi^{\prime})/\sigma(J/\psi+J/\psi)$ calculated at
$\sqrt{s}=13$ TeV as a function of several kinematic variables. The data are
of LHCb[1].

Acknowledgements. A.V.L. and M.A.M. would like to thank School of Physics and
Astronomy, Sun Yat-sen University (Zhuhai, China) for warm hospitality. P.M.Z.
was partially supported by the National Natural Science Foundation of China
(Grant No. $12375084$).

## References

  * [1] LHCb Collaboration, arXiv:2311.15921 [hep-ex].
  * [2] LHCb Collaboration, JHEP 03, 088 (2024).
  * [3] G. Bodwin, E. Braaten, G. Lepage, Phys. Rev. D 51, 1125 (1995).
  * [4] P. Cho, A.K. Leibovich, Phys. Rev. D 53, 150 (1996);
P. Cho, A.K. Leibovich, Phys. Rev. D 53, 6203 (1996).

  * [5] J.-P. Lansberg, H.-S. Shao, N. Yamanaka, Y.-J. Zhang, Eur. Phys. J. C 79, 1006 (2019).
  * [6] H.-S. Shao, Comput. Phys. Commun. 184, 2562 (2013).
  * [7] H.-S. Shao, Comput. Phys. Commun. 198, 238 (2016).
  * [8] Z.-G. He, B.A. Kniehl, Phys. Rev. Lett. 115, 022002 (2015).
  * [9] J.-P. Lansberg, H.-S. Shao, Phys. Rev. Lett. 111, 122001 (2013).
  * [10] L.V. Gribov, E.M. Levin, M.G. Ryskin, Phys. Rep. 100, 1 (1983);
E.M. Levin, M.G. Ryskin, Yu.M. Shabelsky, A.G. Shuvaev, Sov. J. Nucl. Phys.
53, 657 (1991).

  * [11] S. Catani, M. Ciafaloni, F. Hautmann, Nucl. Phys. B 242, 97 (1990);
S. Catani, M. Ciafaloni, F. Hautmann, Nucl. Phys. B 366, 135 (1991);
J.C. Collins, R.K. Ellis, Nucl. Phys. B 360, 3 (1991).

  * [12] R. Angeles-Martinez, A. Bacchetta, I.I. Balitsky, D. Boer, M. Boglione, R. Boussarie, F.A. Ceccopieri, I.O. Cherednikov, P. Connor, M.G. Echevarria, G. Ferrera, J. Grados Luyando, F. Hautmann, H. Jung, T. Kasemets, K. Kutak, J.P. Lansberg, A. Lelek, G.I. Lykasov, J.D. Madrigal Martinez, P.J. Mulders, E.R. Nocera, E. Petreska, C. Pisano, R. Placakyte, V. Radescu, M. Radici, G. Schnell, I. Scimemi, A. Signori, L. Szymanowski, S. Taheri Monfared, F.F. Van der Veken, H.J. van Haevermaet, P. Van Mechelen, A.A. Vladimirov, S. Wallon, Acta Phys. Polon. B 46, 2501 (2015).
  * [13] S.P. Baranov, Phys. Rev. D 84, 054012 (2011).
  * [14] A.A. Prokhorov, A.V. Lipatov, M.A. Malyshev, S.P. Baranov, Eur. Phys. J. C 80, 1046 (2020).
  * [15] S.P. Baranov, A.V. Lipatov, A.A. Prokhorov, Phys. Rev. D 106, 034020 (2022).
  * [16] Z.-G. He, B.A. Kniehl, M.A. Nefedov, V.A. Saleev, Phys. Rev. Lett. 123, 162002 (2019).
  * [17] E.A. Kuraev, L.N. Lipatov, V.S. Fadin, Sov. Phys. JETP 44, 443 (1976);
E.A. Kuraev, L.N. Lipatov, V.S. Fadin, Sov. Phys. JETP 45, 199 (1977);
I.I. Balitsky, L.N. Lipatov, Sov. J. Nucl. Phys. 28, 822 (1978).

  * [18] M. Ciafaloni, Nucl. Phys. B 296, 49 (1988);
S. Catani, F. Fiorani, G. Marchesini, Phys. Lett. B 234, 339 (1990);
S. Catani, F. Fiorani, G. Marchesini, Nucl. Phys. B 336, 18 (1990);
G. Marchesini, Nucl. Phys. B 445, 49 (1995).

  * [19] C.-H. Chang, Nucl. Phys. B 172, 425 (1980).
  * [20] E.L. Berger, D. Jones, Phys. Rev. D 23, 1521 (1981).
  * [21] R. Baier, R. Rückl, Phys. Lett. B 102, 364 (1981).
  * [22] H. Krasemann, Z. Phys. C 1, 189 (1979).
  * [23] G. Guberina, J. Kühn, R. Peccei, R. Rückl, Nucl. Phys. B 174, 317 (1980).
  * [24] F. Hautmann, H. Jung, Nucl. Phys. B 883, 1 (2014).
  * [25] A.V. Lipatov, G.I. Lykasov, M.A. Malyshev, Phys. Rev. D 107, 014022 (2023).
  * [26] K. Golec-Biernat, L. Motyka, T. Stebel, Phys. Rev. D 103, 034013 (2021).
  * [27] A.V. Lipatov, M.A. Malyshev, H. Jung, Phys. Rev. D 100, 034028 (2019).
  * [28] A.V. Lipatov, M.A. Malyshev, Phys. Rev. D 103, 094021 (2020).
  * [29] A.V. Lipatov, G.I. Lykasov, M.A. Malyshev, Phys. Lett. B 839, 137780 (2023).
  * [30] A.V. Lipatov, M.A. Malyshev, Phys. Rev. D 108, 014022 (2023).
  * [31] V.A. Bednyakov, G.I. Lykasov, V.V. Lyubushkin, Europhys. Lett. 92, 31001 (2010).
  * [32] V.A. Bednyakov, A.A. Grinyuk, G.I. Lykasov, M. Poghosyan, Int. J. Mod. Phys. A 27, 1250042 (2012).
  * [33] N.A. Abdulov, A. Bacchetta, S.P. Baranov, A. Bermudez Martinez, V. Bertone, C. Bissolotti, V. Candelise, L.I. Estevez Banos, M. Bury, P.L.S. Connor, L. Favart, F. Guzman, F. Hautmann, M. Hentschinski, H. Jung, L. Keersmaekers, A.V. Kotikov, A. Kusina, K. Kutak, A. Lelek, J. Lidrych, A.V. Lipatov, G.I. Lykasov, M.A. Malyshev, M. Mendizabal, S. Prestel, S. Sadeghi Barzani, S. Sapeta, M. Schmitz, A. Signori, G. Sorrentino, S. Taheri Monfared, A. van Hameren, A.M. van Kampen, M. Vanden Bemden, A. Vladimirov, Q. Wang, H. Yang, Eur. Phys. J. C 81, 752 (2021).
  * [34] A.V. Lipatov, M.A. Malyshev, S.P. Baranov, Eur. Phys. J. C 80, 330 (2020).
  * [35] S.P. Baranov, A.M. Snigirev, N.P. Zotov, Phys. Lett. B 705, 116 (2011).
  * [36] C.H. Kom, A. Kulesza, W.J. Stirling, Phys. Rev. Lett. B 107, 082002 (2011).
  * [37] S.P. Baranov, A.M. Snigirev, N.P. Zotov, A. Szczurek, W. Schäfer, Phys. Rev. D 87, 034035 (2013).
  * [38] S.P. Baranov, A.V. Lipatov, Phys. Rev. D 100, 114021 (2019).
  * [39] D0 Collaboration, Phys. Rev. D 90, 111101(R) (2014).
  * [40] D0 Collaboration, Phys. Rev. Lett 116, 082002 (2016).
  * [41] ATLAS Collaboration, Eur. Phys. J. C 77, 76 (2017).
  * [42] CMS Collaboration, JHEP 1705, 013 (2017).
  * [43] PDG Collaboration, Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * [44] S.P. Baranov, A.V. Lipatov, M.A. Malyshev, Eur. Phys. J. C 78, 820 (2018).


11institutetext: Jožef Stefan International Postgraduate School, Ljubljana,
Slovenia 22institutetext: Jožef Stefan Institute, Ljubljana, Slovenia
33institutetext: Institute of Pathophysiology, Faculty of Medicine, University
of Ljubljana, Ljubljana, Slovenia 44institutetext: University of Nova Gorica,
Vipava, Slovenia 55institutetext: Temida d.o.o, Ljubljana, Slovenia

#  COVID-19 therapy target discovery with context-aware literature mining

Matej Martinc 1122 Blaž Škrlj 1122 Sergej Pirkmajer 33 Nada Lavrač 112244
Bojan Cestnik 5511 Martin Marzidovšek 1122 Senja Pollak 22

###### Abstract

The abundance of literature related to the widespread COVID-19 pandemic is
beyond manual inspection of a single expert. Development of systems, capable
of automatically processing tens of thousands of scientific publications with
the aim to enrich existing empirical evidence with literature-based
associations is challenging and relevant. We propose a system for
contextualization of empirical expression data by approximating relations
between entities, for which representations were learned from one of the
largest COVID-19-related literature corpora. In order to exploit a larger
scientific context by transfer learning, we propose a novel embedding
generation technique that leverages SciBERT language model pretrained on a
large multi-domain corpus of scientific publications and fine-tuned for domain
adaptation on the CORD-19 dataset. The conducted manual evaluation by the
medical expert and the quantitative evaluation based on therapy targets
identified in the related work suggest that the proposed method can be
successfully employed for COVID-19 therapy target discovery and that it
outperforms the baseline FastText method by a large margin.

###### Keywords:

Knowledge discovery Literature mining Representation learning Contextual
embeddings COVID-19.

The final reviewed publication was published in Proceedings of the 23rd
International Conference on Discovery Science (DS 2020), Thessaloniki, Greece,
October 19–21, 2020 and is available online at
https://doi.org/10.1007/978-3-030-61527-7˙8.

##  1 Introduction

Scientific knowledge for a specific domain is in most cases given in an
unstructured form, as a set of scientific papers covering a variety of
findings, experiments and methodologies related to a specific scientific field
or problem. The current speed and quantity of scientific research production
makes manual inspection of the literature from a specific field virtually
impossible. The recent trend of inter-disciplinary research complicates things
even more, as it would require from a researcher to understand all the
aspects, from which a specific research problem can be covered in order to
“connect all the dots” and advance the field by the discovery of the so-called
latent scientific knowledge.

To solve this problem, several automated strategies for uncovering this
knowledge have been proposed. Somewhat older studies proposed literature-based
discovery (LBD) [8] focusing especially on cross-domain literature mining,
which aims at finding interesting bridging terms (b-terms) or bridging links
revealing the potentially new connections between separate domain corpora of
interest. On the other hand, more recent approaches to latent knowledge
discovery from the scientific literature employ word embeddings [26]. For
example, a study by [34] showed that latent knowledge regarding future
discoveries is to a large extent embedded in past publications by retrieving
information from the scientific literature with the usage of Word2Vec
embeddings [26].

The latest development in the natural language processing (NLP) is a new type
of embeddings called contextual embeddings. ELMo (Embeddings from Language
Models) [29] and BERT (Bidirectional Encoder Representations from
Transformers) [12] are the most prominent representatives of this type of
contextual embeddings, and have been also adapted to scientific literature
[3]. The main difference between these novel contextual embeddings and older
“static” embeddings is that in these embeddings a different vector is
generated for each context a word appears in, i.e., for each specific word
usage in the corpus. These new contextual embeddings solve the problems with
word polysemy and other changes in word meaning given different context. On
the other hand, it is not entirely clear how to generate a meaningful general
word representation from the word usage embeddings. This means that the usage
of contextual embeddings for LBD is not entirely straight forward, since they
can not be used in the same way as the traditional static embeddings, and have
at least to our knowledge not been used for the task at hand.

In this work, we explore how contextual embeddings can be leveraged for the
task of discovering latent scientific knowledge in the very topical scientific
literature about the COVID-19 disease. More specifically, we are interested in
the discovery of new COVID-19 therapy targets from the targets discovered in
the past research. The novelty of this work is two-fold:

  * •

The paper contributes a new methodology of generating general word
representations from contextual embeddings, proposes an entire workflow for
acquisition of novel COVID-19 therapy targets and shows that our method of
using contextual embeddings for LBD outperforms the baseline method of using
static embeddings by a large margin.

  * •

Medically, the paper contributes to identifying new potential COVID-19 therapy
targets, motivated by a recent proof-of-concept study that used a state-of-
the-art omics approach to identify new possible targets for existing drugs,
such as ribavirin [5].

##  2 COVID-19 medical background and recent therapy targets

In late 2019 a novel coronavirus disease (COVID-19), caused by severe acute
respiratory syndrome coronavirus 2 (SARS-CoV-2), emerged in China [38, 39].
COVID-19 quickly spread and was declared a pandemic by the World Health
Organization.

While new targeted therapies and vaccines against SARS-CoV-2 virus are being
actively developed, their potential use in the clinics is not imminent.
Therefore, until effective pharmacological therapies and/or vaccines are
available, medicine needs to resort to other approaches to treat patients with
COVID-19 or prevent transmission of SARS-CoV-2. One approach is to identify
which among the antiviral drugs that were developed to treat other viral
diseases might be effective against SARS-CoV-2. A preliminary report suggests
that remdesivir seems to be the most promising candidate among these drugs
[2]. Another approach is to identify drugs that are used for other purposes
but also exert antiviral effects. The most prominent example among these is
hydroxychloroquine, which is used for chronic treatment of rheumatic diseases
but also suppresses SARS-CoV-2 in vitro [22]. Identifying a known drug with
well-characterized adverse effects would certainly save time and lives before
more specific treatments are developed. However, repurposing of existing drugs
is also a challenge as highlighted by a recent controversy with
hydroxychloroquine [7, 25] and new candidate drugs and/or therapeutic targets
are needed.

##  3 Related work

The related work is divided into three Sections, namely related work on
Literature-based discovery in Section 3.1, related work on text representation
learning in Section 3.2 and selected overview of recent NLP research on
COVID-19 in Section 3.3.

###  3.1 Literature-based discovery

Literature-based discovery (LBD) aims to generate new knowledge by combining
what is already known in the literature. It has been used to (semi-
automatically) identify new connections between genes, drugs and diseases,
etc. [18]. Traditionally, LBD has been addressed as finding interesting
bridging terms revealing the potentially new connections between separate
domain corpora of interest [8]. Swanson [33] developed one of the early LBD
approaches, the so-called ABC model, to detecting interesting b-terms to
uncover the possible cross-domain relations among previously unrelated
concepts.

On the other hand, a more recent state-of-the art tool LION LBD [31] enables
researchers to navigate published information and supports hypothesis
generation and testing. The system is built with a particular focus on the
molecular biology of cancer. LBD has led to discovery of potential treatments
in other domains, including multiple sclerosis [19], and has been applied
successfully in drug development and repurpusing [11]. Recent LBD approaches
benefit from word embeddings. One is the study by [34] already mentioned in
Section 1 and the other is the work by [9], who proposed graph-based, neural
network methods to perform open and closed LBD and demonstrated improved
performance on existing tasks.

###  3.2 Text representation and embeddings

Recently, the embedding approach became a prevalent way to build
representations for many different types of entities, e.g., texts, graphs,
electronic health records, images, relations, recommendations, etc. Text
embeddings use large corpora of documents to extract vector representations
for words, sentences, and documents. The first neural word embeddings like
Word2vec [26] produced one vector for each word, irrespective of its polysemy.
These so-called static embeddings have been further developed and the most
popular static embeddings currently in use besides Word2Vec are GloVe (Global
vectors for word representation) [28] and FastText [4]. Recent developments
like ELMo [29] and BERT [12] take a context of a sentence into account and
produce different word vectors for different contexts of each word. Another
novelty of these approaches is the employment of the transfer learning
technique, which has recently become a well established procedure in the field
of NLP. This procedure relies on a language model pretraining on very large
unlabeled textual resources and after that transfer of the knowledge obtained
by the language model onto a specific downstream task by further fine-tuning
the model.

###  3.3 Text mining and NLP research related to COVID-19

With regard to biomedical research on COVID-19, time is a central factor as
scientists try to design treatments and vaccines amid the pandemic caused by
the SARS-CoV-2 virus, therefore leveraging LBD and its potential to reduce
scientific discovery time could prove crucial.

Many search platforms emerged for retrieving COVID-19 related papers. For
example, Neural Covidex111https://covidex.ai/ is based on neural ranking
architecture and provides information access capabilities to the COVID-19 Open
Research Dataset (CORD-19) (see Section 4.1). SciSight [17] in contrast to
standard targeted search facilitates finding connections between biomedical
concepts that are not obvious from reading individual papers. It displays a
network of top related terms mined from the corpus, based on the co-appearance
in the same sentence.

Studies that can generate new knowledge about COVID-19 by applying embeddings
are still scarce but do exist. For example, a recent study has projected
Covid-related medical texts in a 3D human atlas space that helps to navigate
the literature [14]. The objective was to learn semantically aware groundings
of sentences with five different BERT models [12].

##  4 Background knowledge and resources

We describe the CORD-19 corpus (Section 4.1) and embeddings technology
(Section 4.2) used in this study.

###  4.1 CORD19 database

The scientific literature considered in this work has been recently introduced
as the CORD-19 corpus222https://www.kaggle.com/allen-institute-for-
ai/CORD-19-research-challenge. CORD-19 is a resource of over 135,000 scholarly
articles, including over 68,000 with full text, about COVID-19, SARS-CoV-2,
and related coronaviruses. This freely available data set is provided to the
global research community to apply recent advances in NLP and other AI
techniques to generate new insights in support of the ongoing fight against
this infectious disease.

We use the corpus version 12, published on May 1st 2020, from which we extract
only full text scholarly articles converted into xml from a pdf format. This
results in altogether 48,410 papers, which are summarized in Table 1.

Table 1: CORD-19 dataset statistics.

Origin | Number of papers | Number of tokens
---|---|---
Commercial use subset | 9,918 | 46,206,453
Non-commercial use subset | 2,584 | 10,732,608
PMC custom license subset | 32,450 | 156,247,363
bioRxiv (not peers reviewed) | 2,670 | 8,968,183
medRxiv subset (not peer reviewed) | 788 | 3,285,558
All | 48,410 | 225,440,165

###  4.2 Considered embeddings

We use FastText [4] embeddings as a baseline in this study. The main advantage
of FastText embeddings is its word representation as a sum of n-grams, which
allows the model to, in addition to leveraging semantic relations, also
leverage morphological information.

One of the most oftenly used models for the generation of contextual
embeddings is the BERT model [12] that was originally pretrained on the Google
Books Corpus (800 million tokens) and Wikipedia (2,500 million tokens). This
pretraining is however not entirely appropriate for the text mining tasks on
the scientific literature due to specificities of the scientific language and
vocabulary. For this reason, in this research we opted for SciBERT [3], a
version of BERT pretrained on a large multi-domain corpus of scientific
publications, a random sample of 1.14M papers from Semantic Scholar. SciBERT
model has 12 encoder layers with the attention mechanism and a hidden layer
size of 768.

##  5 Methodology

In this section, we present the methodology of the proposed approach by
explaining how we obtain word representations, how we acquire therapy target
candidates and how we evaluate the approach.

###  5.1 Word representations

First, we fine-tune SciBERT as a masked language model for domain adaptation
on the lowercased CORD-19 dataset. Next, we generate word representations for
each word in the vocabulary. Figure 1 visualizes the process described below.
The documents from the corpus are split into sequences of byte-pair encoded
tokens [20] of a maximum length of 256 tokens and fed into the fine-tuned
SciBERT model. For each of these sequences of length $n$, we create a sequence
embedding by summing the last four encoder output layers. The resulting
sequence embedding of size $n$ times embeddings size represents a
concatenation of contextual embeddings for the $n$ tokens in the input
sequence. By chopping it into $n$ pieces, we acquire a representation, i.e. a
contextual token embedding, for each word used in the corpus. Note that these
representations vary according to the context in which the token appears,
meaning that the same word has a different representation in each specific
context (sequence).

Finally, the resulting embeddings are aggregated on the token level (i.e. for
every token in the corpus vocabulary, we create a list of all their contextual
embeddings) and are averaged, in order to get one representation for each
token in the vocabulary. We enforce a constraint that a list of contextual
embeddings for a specific token should contain at least five elements,
otherwise the specific token is discarded. This is done in order to remove
tokens that do no appear in the corpus enough times for the model to learn a
meaningful representation (e.g., mostly tokens that contain typos or very rare
technical terms). Since the byte-pair input encoding scheme [20] employed by
the SciBERT model does not necessarily generate tokens that correspond to
words but rather generate tokens that correspond to parts of words, we also
propose the following on the fly reconstruction mechanism that allows us to
get word representations from byte pair tokens. If a word is split into more
than one byte pair token, we take an embedding for each byte pair token
constituting a word and build a word embedding by averaging these byte pair
tokens. The resulting average is used as a context specific word
representation.

The final result are static embeddings for each word in the vocabulary,
capable of leveraging a broader semantic knowledge due to the SciBERT being
pretrained on a large corpus of scientific articles. As a baseline, we also
train a FastText skip-gram model with an embedding dimension of 100 (which is
the default) on the lowercased CORD-19 dataset. Once again we enforce the
constraint that a word should appear in the corpus at least five times.

Figure 1: Extraction of word usage embeddings from BERT. Note that only the
last 4 out of 12 BERT encoder layers are used for the embedding generation.
This was done in accordance with the previous studies that suggested that the
last four layers carry the bulk of the semantic information obtained by the
model [24].

###  5.2 Synonym resolution

Once embeddings are generated, we conduct synonym resolution with the help of
a list of 19,302 gene names and their most common synonyms [30]. The
embedddings belonging to the synonyms of the same gene are averaged in order
to combine contextual information of different identifiers referring to the
same gene and in order to avoid possible mismatches due to different naming.

###  5.3 Candidate acquisition

The main idea of our approach is to leverage semantic similarity in order to
derive new scientific knowledge from an already existing one. For this to
work, some initial seed concepts need to be acquired and used as a starting
point. We explore two possibilities for this:

  * •

Seed concepts recommended by the expert: The experts with a medical background
were asked to recommend genes and/or proteins with a known and confirmed link
to COVID-19. The final consensus was to focus on angiotensin-converting enzyme
2 (ACE2) and transmembrane protease serine 2 (TMPRSS2). ACE2, a receptor for
the spike S protein, is important because SARS-CoV-2 uses it to enter the host
cell [16]. TMPRSS2 promotes SARS-CoV-2 entry into the cell by priming the
spike S protein [16]. Blockage of binding of SARS-CoV-2 to ACE2 or inhibition
of TMPRSS2 are therefore two possible approaches to treat COVID-19.

  * •

Seed concepts found in the literature: Due to the abundance of recent research
on COVID-19 it is also possible to find seed concepts in the related research.
We opted for a study by [5] in which a set of COVID-19 therapy targets were
identified. The considered list of altogether 2802 potential targets333Note
that the original list contains 2715 targets (see Supplementary Table 1 in
[5]). Some of them are however represented as a set of similar genes/proteins
belonging to the same family. On the other hand, we treat each individual
gene/protein as a separate target, which results in a set of 2802 targets. is
the result of a large-scale screening for active proteins, and offers a
starting set of candidates obtained empirically. The list is ranked according
to the increase or decrease of production of a specific protein at a specific
time point. We explore what is the optimal number of seed candidates by
exponentially enlarging the size of the seed candidate set. Sampling from the
list is conducted according to the ranking of the protein candidates, i.e., we
sample 2, 4, 8, 16, 32 and 64 best ranked seed candidates according to the
increase in their production 24 hours after the infection (column Ratio 24h in
Supplementary Table 1 in the study by [5]).

Once seed concepts are acquired, we calculate their embeddings and look for
semantically similar concepts by finding the concepts that are the closest to
seed concepts according to the cosine distance between the embeddings444Note
that these concepts obtained according to semantic similarity are not
necessarily proteins/genes but rather any word in the embedding vocabulary..
More specifically, we find a set of 2802 closest candidate concepts for each
gene/protein in each seed candidate set, and the acquired candidates are
ranked according to the cosine similarity. Finally, we calculate the average
ranking for each candidate (i.e. by averaging ranks for each seed concept in
the set) and therefore obtain $\textrm{NumOfCandidatesInSet}*2802$ closest
candidates for each of the seed concept sets with possible duplicates
originating from different seed concepts.

Since the initial experiments showed that many of the most similar concepts
are in fact variations of the same base concept (e.g., the closest neighbours
to ACE2 being ACE, ACE2M, ACE2S…) and since we are interested in maximizing
the variety of the acquired candidates, we conduct an additional filtering
according to the normalized Levenshtein distance defined as:

$\textsc{normLD}=1-\frac{LD}{\max(\textrm{len}(w_{1}),\textrm{len}(w_{2}))},$

where normLD stands for normalized Levenshtein distance, $LD$ for Levenshtein
distance, $w_{1}$ is either a seed concept or a concept already in the list of
acquired neighbours and $w_{2}$ is the new candidate neighbour. Concepts for
which normalized Levenshtein difference is bigger than 0.7 are discarded555The
normalized Levenshtein difference threshold of 0.7 was chosen empirically..
The filtering is conducted in order from the top of the list (neighbours with
the best average rank) to the bottom.

At the end of the candidate acquisition process, we cut the ranked list of
neighbours at 2802 target candidates for each of the distinct seed concept
sets used in the evaluation.

###  5.4 Evaluation

The methods for discovering new therapy targets are evaluated in two
evaluation settings, quantitative and qualitative.

####  5.4.1 Quantitative Evaluation

We evaluate if therapy target candidates acquired in the previous step have
been confirmed as targets in the study by [5], i.e. how many of them appear in
the list of 2802 candidates they identified666Note that the study by [5] is
not included in the CORD-19 corpus used for training the embeddings, since it
was published on May 14th 2020 and we use the CORD-19 version published on May
1st 2020.. Note that in this setting we only evaluate the proposed method on
the previously existing knowledge, therefore in the quantitative evaluation we
can not asses, if the method has managed to discover some potentially useful
and previously undiscovered knowledge.

We are interested in precision at rank $k$. This means that only the
candidates ranked equal to or higher than $k$ are considered and the rest are
disregarded. Precision is the ratio of the number of relevant candidates
divided by the number of candidates returned by the system, or more formally:

$\textrm{precision}=\frac{|\textrm{relevant}~{}~{}\textrm{candidates}@k|}{|\textrm{returned}~{}~{}\textrm{candidates}|}$

Recall@$k$ is the ratio of the number of relevant candidates ranked equal to
or higher than $k$ by the system divided by the number of correct ground truth
candidates:

$\textrm{recall}=\frac{|\textrm{relevant}~{}~{}\textrm{candidates}@k|}{|\textrm{correct}~{}~{}\textrm{candidates}|}$

We measure precision and recall at k=100 and k=2802 in order to investigate
how different number of retrieved candidates for each seed concept set affects
the precision and recall of the methods. More specifically, we are trying to
confirm or deny a hypothesis that larger k values degrade the overall
precision of the method.

The relevance of the candidate is determined according to two matching
criteria. First one is the exact match, where the candidate is deemed relevant
if it appears in the list of identified targets in the study by [5]. The
second is the fuzzy match, where we check if the targets belong to the same
“family” as a specific confirmed target. This strategy was proposed by the
medical experts and checks whether the prefix of the specific gene (characters
in the gene name that appear before the first digit in the name) matches a
prefix of a specific gene name in the list. We enforce an additional
constraint that the matching prefixes need to be at least three characters
long for a successful match in order to minimize the false positive rate.

####  5.4.2 Qualitative Evaluation

We generated two distinct therapy target candidate lists using the proposed
SciBERT based embedding method. First one contained 100 closest neighbours to
the protein ACE2 according to the cosine distance between embeddings, and the
second one contained 100 closest neighbours to the protein TMPRSS2. Both lists
were given to the medical expert who inspected the list for possible
previously undiscovered candidates.

##  6 Results

Here we present the results of the quantitative and qualitative evaluation.

###  6.1 Results of the quantitative evaluation

The results of the quantitative evaluation are presented in Table 2. In column
ACE2 + TMPRSS2 we present results when these two proteins are used as seed
concepts, and in column UBA2 + NCKAP1 we present results when these two
proteins, which were chosen according to the largest value of the Ratio 24h
criterion (see Section 5.3) are used as seed concepts. Left part of the Table
presents results for the proposed approach based on SciBERT and the right part
of the Table presents results for the baseline FastText approach in terms of
precision and recall at two distinct $k$ values (100 and 2802). EXACT
indicates that exact matching is used and FUZZY indicates fuzzy matching (see
Section 5.4.1).

Table 2: Results (precision@k and recall@k) of the quantitative evaluation for
two seeds by the expert and two seeds from the literature. Best result in each
row is bolded.

| SciBERT | FastText
---|---|---
| ACE2 + TMPRSS2 | UBA2 + NCKAP1 | ACE2 + TMPRSS2 | UBA2 + NCKAP1
EXACT P@100 | 0.110 | 0.220 | 0.040 | 0.170
EXACT R@100 | 0.004 | 0.008 | 0.001 | 0.006
EXACT P@2802 | 0.097 | 0.118 | 0.025 | 0.076
EXACT R@2802 | 0.097 | 0.118 | 0.025 | 0.076
FUZZY P@100 | 0.290 | 0.490 | 0.070 | 0.380
FUZZY R@100 | 0.010 | 0.017 | 0.002 | 0.014
FUZZY P@2802 | 0.222 | 0.252 | 0.092 | 0.183
FUZZY R@2802 | 0.222 | 0.252 | 0.092 | 0.183

SciBERT based method outperforms the FastText baseline by a large margin in
both seed therapy target acquisition scenarios and according to all the
criteria. Using UBA2 + NCKAP1 works better than using ACE2 + TMPRSS2,
achieving the best fuzzy precision@100 of 0.490 and the best exact
precision@100 of 0.220. FastText baseline also works fairly well in this
scenario, achieving fuzzy precision@100 of 0.380 and the best exact
precision@100 of 0.170. When more (2802) candidates are obtained, the recall
increases for both methods but at an expense of a significant drop in
precision for both methods and for almost all configurations. The only
exception is the increase in fuzzy precision by about 2 percentage points when
FastText method and ACE2 + TMPRSS2 seed concepts are used. The most likely
reason for the drop is that at larger k values some of the target candidates
acquired by the method might be semantically too dissimilar to the seed
targets, since more candidates per each seed therapy target need to be
acquired in order to get the required amount of semantic neighbours (e.g., for
k=2802, we get about 1401 semantic neighbours for each of the seed genes).

This raises the question of how many seed terms should be supplied to the
system for the best performance when a large number of target candidates is
required as output. Figure 2 shows the relation between the achieved
recall@2802 (exact and fuzzy) of both methods when we increase the number of
seed candidates (see Section 5.3 for details about our sampling procedure).
For SciBERT based method, the best fuzzy and exact recalls are achieved when
32 seed candidates are used (28.2% and 14.1% respectively). On the other hand,
the FastText based method shows a spike in performance when 4 seed concepts
are used. This indicates that for some reason the two seed candidates ranked
third and fourth (ENO1 and ATP5O, respectively) according to the Ratio 24
criterion have a very positive effect on the FastText model but not SciBERT.
While we do not have a clear explanation for this phenomenon, it is
hypothesized that it might be connected with morphological similarity between
these two genes and other genes in the list of candidates proposed by [5],
since FastText can also leverage morphological similarity. Spikes asides, the
general trend for both methods and both recalls is quite similar. There is a
gradual increase in performance for up to 32 seed candidates and after that
the performance decreases.

Figure 2: Relation between recall and the number of seed candidates.

###  6.2 Results of the qualitative evaluation

Nine genes/proteins were the same in the ACE2 and TMPRSS2 lists, indicating
they might be important for pathogenesis of COVID-19. The role of these
genes/proteins in pathogenesis of COVID-19 has not been established, but
indirect evidence supports this notion at least for some of them. Indeed, most
of these genes/proteins have been previously linked to viral diseases,
including those caused by SARS-CoV (a virus, which causes SARS, and is related
to SARS-CoV-2), and other coronaviruses (Table 3). Furthermore, METAP2 and
DPP7, which we identified as potentially relevant for COVID-19, were altered
in cells infected with SARS-CoV-2, although the difference for DPP7 did not
reach the level of statistical significance [5].

Interestingly, three proteins in Table 3 (PTGS2, CRTH2, and PLA2R1) are linked
to infection with coronaviruses as well as metabolism of phospholipids and/or
prostaglandin synthesis and action. Furthermore, both the ACE2 and TMPRSS2
lists contain genes/proteins, such as PLA2 (phospholipase A2, PLA2G2D (Group
IID secretory phospholipase A2), and SPLA2 (secretory PLA2), which do not
match directly, but are involved in the same or related cellular processes.
Notably, increased expression of Pla2g2d in older mice was shown to be linked
with increased mortality due to SARS-CoV infection [36]. In addition, a recent
proteomic analysis has demonstrated that protein abundance of PLAA
(phospholipase A2-activating protein), PLA2G4A (cytosolic phospholipase A2),
and PLA2G2 (Group IIA phospholipase A2) is altered in cultured cells infected
with SARS-CoV-2 [5], which gives further credence to the idea that
phospholipid metabolism is important under these conditions. In summary, taken
together with published experimental data, our analysis suggests that
phospholipases and/or prostaglandins might represent a target for treatment of
COVID-19.

Table 3: Genes/proteins (in alphabetical order) which are common to the
TMPRSS2 and ACE2 list and their (putative) relevance to COVID-19.

Gene | Protein | Relevance to COVID-19
---|---|---
ATP2B2 (PMCA2) | Plasma membrane $Ca^{2+}$-transporting ATPase | ?
CRTH2 (PTGDR2) | Prostaglandin D2 (PGD2) receptor |  | PGD2 is important for survival of mice infected
---
with neurotropic coronavirus. Increased production of PGD2
is linked to increased mortality in aged mice. PGD2blockade
improves survival in mice infected with SARS-CoV [35, 37].
DPP7 (DPP2) | Dipeptidyl peptidase 2 |  | DPP7 is associated with the magnitude of
---
the antibody response to influenza vaccination [15].
MECP2 | Methyl-CpG-binding protein 2 |  | MECP2 duplication in humans is associated with IgA/IgG2
---
antibody deficiency and severe infections. Mice overexpressing
MECP2 are hypersensitive to influenza A virus [1, 10].
METAP2 (P67EIF2) |  | Methionine aminopeptidase 2 (Initiation
---
factor 2-associated 67 kDa glycoprotein)
| Plays a role in regulation of protein
---
synthesis during vaccinia virus infection [6].
PLA2R1 | Secretory phospholipase A2 (PLA2) receptor |  | Restricted activity of PLA2 is associated with improved
---
survival in mice infected with HCoV-OC43. Inhibition
of cytosolic PLA2 suppresses replication of HCoV-229E [13, 27].
PTGS2 (COX2) | Prostaglandin G/H synthase 2 (cyclooxygenase-2) | SARS-CoV induces cyclooxygenase 2 [23].
SOX2 | Transcription factor SOX-2 |  | SOX2+ cells are important for regeneration of airway
---
epithelium after severe influenza infection in mice [32].
SSTR2 (SST2) | Somatostatin receptor type 2 | ?

##  7 Conclusions and further work

In this paper we presented a method for discovering new COVID-19 therapy
targets by leveraging contextual embeddings, which outperforms the method
based on FastText embeddings. We explored the best tactics for acquiring seed
targets from the related work if expert knowledge is not available. The
results of the manual qualitative evaluation by the expert indicate that at
least two groups of novel therapy target candidates have been discovered.

The proposed method outperforms the baseline FastText method by a large
margin, which can be explained by the fact that SciBERT is also leveraging
knowledge gained during the pretraining on the large corpus of scientific
literature, which enables the model to generate vector representations that
reflect this wider semantic context. The drawback is however the difference in
the amount of computational resources required by the two methods. We also
acknowledge that the proposed method, which constructs static embeddings from
the SciBERT contextual embeddings is not the only possibility for construction
of meaningful semantic representations. Other possibilities and models (e.g.,
BioBERT [21]) will be explored in the future work. The quantitative evaluation
indicates that the precision and recall of the method are still relatively low
in most cases. This can on one side indicate that COVID-19 topic is not
researched enough to confirm relations between COVID-19 and some candidates
found by the proposed method. Another indication of this is the qualitative
study, which confirmed that some of the proposed candidates found by the
system have research potential but have not yet been explicitly confirmed as
being related to COVID-19 in the existing literature.

On the other hand, low precision most likely also indicates that there is
still a large amount of proposed candidates, which play no role in the
advancement and prevention of the COVID-19 disease. Some of these false
positives can be attributed to inadequate synonym resolution since the list
used for that task (see Section 5.2) most likely covers only a small
percentage of genes and their synonyms found in the CORD-19 corpus. Other
mistakes can be contributed to the byte pair encoding scheme SciBERT employs.
Since the model generates embeddings for subword tokens instead for an entire
words (see how we deal with this problem in Section 5.1), some words with
similar roots or affixes can perhaps appear closer in the semantic space as
they should according to their semantic relatedness because of the
morphological resemblance. We will address this issues in the future work.

#### Acknowledgements

This work was funded by the Slovenian Research Agency (ARRS) through core
research programme _Knowledge Technologies_ (P2-0103), research project
_Semantic Data Mining for Linked Open Data_ (financed under the ERC
Complementary Scheme, N2-0078), and a young researcher grant (BŠ). The work
was also supported by EU Horizon 2020 research and innovation programme under
grant agreement No 825153, project EMBEDDIA (Cross-Lingual Embeddings for
Less-Represented Languages in European News Media). The publication reflects
only the authors’ views and the EC is not responsible for any use that may be
made of the information it contains.

## References

  * [1] Bauer, M., Kölsch, U., Krüger, R., Unterwalder, N., Hameister, K., Kaiser, F.M., Vignoli, A., Rossi, R., Botella, M.P., Budisteanu, M., et al.: Infectious and immunologic phenotype of mecp2 duplication syndrome. Journal of Clinical Immunology 35(2), 168–181 (2015)
  * [2] Beigel, J.H., Tomashek, K.M., Dodd, L.E., Mehta, A.K., Zingman, B.S., Kalil, A.C., Hohmann, E., Chu, H.Y., Luetkemeyer, A., Kline, S., et al.: Remdesivir for the treatment of covid-19—preliminary report. New England Journal of Medicine (2020)
  * [3] Beltagy, I., Cohan, A., Lo, K.: Scibert: Pretrained contextualized embeddings for scientific text. arXiv preprint arXiv:1903.10676 (2019)
  * [4] Bojanowski, P., Grave, E., Joulin, A., Mikolov, T.: Enriching word vectors with subword information. Transactions of the Association for Computational Linguistics 5, 135–146 (Dec 2017)
  * [5] Bojkova, D., Klann, K., Koch, B., Widera, M., Krause, D., Ciesek, S., Cinatl, J., Münch, C.: Proteomics of sars-cov-2-infected host cells reveals therapy targets. Nature pp. 1–8 (May 2020). https://doi.org/10.1038/s41586-020-2332-7
  * [6] Bose, A., Saha, D., Gupta, N.K.: Viral infection: I. Regulation of protein synthesis during vaccinia viral infection of animal cells. Archives of Biochemistry and Biophysics 342(2), 362–372 (1997)
  * [7] Boulware, D.R., Pullen, M.F., Bangdiwala, A.S., Pastick, K.A., Lofgren, S.M., Okafor, E.C., Skipper, C.P., Nascene, A.A., Nicol, M.R., Abassi, M., et al.: A randomized trial of hydroxychloroquine as postexposure prophylaxis for covid-19. New England Journal of Medicine (2020)
  * [8] Bruza, P., Weeber, M.: Literature-based Discovery. Springer Science & Business Media (2008)
  * [9] Crichton, G., Baker, S., Guo, Y., Korhonen, A.: Neural networks for open and closed literature-based discovery. PLOS ONE 15(5), 1–16 (05 2020)
  * [10] Cronk, J.C., Herz, J., Kim, T.S., Louveau, A., Moser, E.K., Sharma, A.K., Smirnov, I., Tung, K.S., Braciale, T.J., Kipnis, J.: Influenza a induces dysfunctional immunity and death in mecp2-overexpressing mice. JCI Insight 2(2) (2017)
  * [11] Deftereos, S.N., Andronis, C., Friedla, E.J., Persidis, A., Persidis, A.: Drug repurposing and adverse event prediction using high-throughput literature analysis. Wiley Interdisc. Reviews: Systems Biology and Medicine 3(3), 323–334 (2011)
  * [12] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint: 1810.04805 (2018)
  * [13] Do Carmo, S., Jacomy, H., Talbot, P.J., Rassart, E.: Neuroprotective effect of apolipoprotein d against human coronavirus oc43-induced encephalitis in mice. Journal of Neuroscience 28(41), 10330–10338 (2008)
  * [14] Grujicic, D., Radevski, G., Tuytelaars, T., Blaschko, M.B.: Self-supervised context-aware covid-19 document exploration through atlas grounding (2020)
  * [15] HIPC-I Consortium, et al.: Multicohort analysis reveals baseline transcriptional predictors of influenza vaccination responses. Science immunology 2(14) (2017)
  * [16] Hoffmann, M., Kleine-Weber, H., Schroeder, S., Krüger, N., Herrler, T., Erichsen, S., Schiergens, T.S., Herrler, G., Wu, N.H., Nitsche, A., et al.: Sars-cov-2 cell entry depends on ace2 and tmprss2 and is blocked by a clinically proven protease inhibitor. Cell (2020)
  * [17] Hope, T., Portenoy, J., Vasan, K., Borchardt, J., Horvitz, E., Weld, D.S., Hearst, M.A., West, J.: Scisight: Combining faceted navigation and research group detection for covid-19 exploratory scientific search. arXiv preprint: 2005.12668 (2020)
  * [18] Korhonen, A., Guo, Y., Baker, S., Yetisgen-Yildiz, M., Stenius, U., Narita, M., Liò, P.: Improving literature-based discovery with advanced text mining. In: DI Serio, C., Liò, P., Nonis, A., Tagliaferri, R. (eds.) Computational Intelligence Methods for Bioinformatics and Biostatistics. pp. 89–98. Springer, Cham (2015)
  * [19] Kostoff, R.N., Briggs, M.B., Lyons, T.J.: Literature-related discovery (LRD): Potential treatments for multiple sclerosis. Technological Forecasting and Social Change 75(2), 239–255 (2008)
  * [20] Kudo, T., Richardson, J.: Sentencepiece: A simple and language independent subword tokenizer and detokenizer for neural text processing. arXiv preprint:1808.06226 (2018)
  * [21] Lee, J., Yoon, W., Kim, S., Kim, D., Kim, S., So, C.H., Kang, J.: BioBERT: a pre-trained biomedical language representation model for biomedical text mining. Bioinformatics 36(4), 1234–1240 (2020)
  * [22] Liu, J., Cao, R., Xu, M., Wang, X., Zhang, H., Hu, H., Li, Y., Hu, Z., Zhong, W., Wang, M.: Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting sars-cov-2 infection in vitro. Cell Discovery 6(1), 1–4 (2020)
  * [23] Liu, M., Gu, C., Wu, J., Zhu, Y.: Amino acids 1 to 422 of the spike protein of sars associated coronavirus are required for induction of cyclooxygenase-2. Virus genes 33(3), 309–317 (2006)
  * [24] Martinc, M., Novak, P.K., Pollak, S.: Leveraging contextual embeddings for detecting diachronic semantic shift. arXiv preprint arXiv:1912.01072 (2019)
  * [25] Mehra, M.R., Desai, S.S., Kuy, S., Henry, T.D., Patel, A.N.: Retraction: Cardiovascular disease, drug therapy, and mortality in covid-19. New England Journal of Medicine (2020)
  * [26] Mikolov, T., Chen, K., Corrado, G., Dean, J.: Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781 (2013)
  * [27] Müller, C., Hardt, M., Schwudke, D., Neuman, B.W., Pleschka, S., Ziebuhr, J.: Inhibition of cytosolic phospholipase a2$\alpha$ impairs an early step of coronavirus replication in cell culture. Journal of virology 92(4), e01463–17 (2018)
  * [28] Pennington, J., Socher, R., Manning, C.D.: Glove: Global vectors for word representation. In: Proc. of the 2014 conference on empirical methods in natural language processing (EMNLP). pp. 1532–1543 (2014)
  * [29] Peters, M.E., Neumann, M., Iyyer, M., Gardner, M., Clark, C., Lee, K., Zettlemoyer, L.: Deep contextualized word representations. arXiv preprint:1802.05365 (2018)
  * [30] Povey, S., Lovering, R., Bruford, E., Wright, M., Lush, M., Wain, H.: The hugo gene nomenclature committee (hgnc). Human genetics 109(6), 678–680 (2001)
  * [31] Pyysalo, S., Baker, S., Ali, I., Haselwimmer, S., Shah, T., Young, A., Guo, Y., Högberg, J., Stenius, U., Narita, M., Korhonen, A.: LION LBD: a literature-based discovery system for cancer biology. Bioinformatics 35(9), 1553–1561 (10 2018)
  * [32] Ray, S., Chiba, N., Yao, C., Guan, X., McConnell, A.M., Brockway, B., Que, L., McQualter, J.L., Stripp, B.R.: Rare sox2+ airway progenitor cells generate krt5+ cells that repopulate damaged alveolar parenchyma following influenza virus infection. Stem cell reports 7(5), 817–825 (2016)
  * [33] Swanson, D.R.: Medical literature as a potential source of new knowledge. Bulletin of the Medical Library Association 78(1), 29 (1990)
  * [34] Tshitoyan, V., Dagdelen, J., Weston, L., Dunn, A., Rong, Z., Kononova, O., Persson, K.A., Ceder, G., Jain, A.: Unsupervised word embeddings capture latent knowledge from materials science literature. Nature 571(7763), 95–98 (2019)
  * [35] Vijay, R., Fehr, A.R., Janowski, A.M., Athmer, J., Wheeler, D.L., Grunewald, M., Sompallae, R., Kurup, S.P., Meyerholz, D.K., Sutterwala, F.S., et al.: Virus-induced inflammasome activation is suppressed by prostaglandin d2/dp1 signaling. Proceedings of the National Academy of Sciences 114(27), E5444–E5453 (2017)
  * [36] Vijay, R., Hua, X., Meyerholz, D.K., Miki, Y., Yamamoto, K., Gelb, M., Murakami, M., Perlman, S.: Critical role of phospholipase a2 group iid in age-related susceptibility to severe acute respiratory syndrome–cov infection. Journal of Experimental Medicine 212(11), 1851–1868 (2015)
  * [37] Zhao, J., Zhao, J., Legge, K., Perlman, S.: Age-related increases in pgd 2 expression impair respiratory dc migration, resulting in diminished t cell responses upon respiratory virus infection in mice. The Journal of clinical investigation 121(12), 4921–4930 (2011)
  * [38] Zhou, P., Yang, X.L., Wang, X.G., Hu, B., Zhang, L., Zhang, W., Si, H.R., Zhu, Y., Li, B., Huang, C.L., et al.: A pneumonia outbreak associated with a new coronavirus of probable bat origin. nature 579(7798), 270–273 (2020)
  * [39] Zhu, N., Zhang, D., Wang, W., Li, X., Yang, B., Song, J., Zhao, X., Huang, B., Shi, W., Lu, R., et al.: A novel coronavirus from patients with pneumonia in china, 2019. New England Journal of Medicine (2020)