††thanks: Corresponding author

# Quantum lattice dynamics and their importance in ternary superhydride
clathrates

Roman Lucrezi  Institute of Theoretical and Computational Physics, Graz
University of Technology, NAWI Graz, 8010 Graz, Austria  Eva Kogler  Institute
of Theoretical and Computational Physics, Graz University of Technology, NAWI
Graz, 8010 Graz, Austria  Simone Di Cataldo  Institute of Theoretical and
Computational Physics, Graz University of Technology, NAWI Graz, 8010 Graz,
Austria  Dipartimento di Fisica, Sapienza Università di Roma, 00185 Rome,
Italy  Markus Aichhorn  Institute of Theoretical and Computational Physics,
Graz University of Technology, NAWI Graz, 8010 Graz, Austria  Lilia Boeri
Dipartimento di Fisica, Sapienza Università di Roma, 00185 Rome, Italy  Enrico
Fermi Research Center, Via Panisperna 89 A, 00184, Rome, Italy  Christoph Heil
<EMAIL_ADDRESS>Institute of Theoretical and Computational Physics,
Graz University of Technology, NAWI Graz, 8010 Graz, Austria

(December 19, 2022)

###### Abstract

Taking a detailed look at the dynamic stability of the recently predicted
BaSiH8 phase [Lucrezi et al., npj Comp. Mat. 8, 119 (2022)], we find that the
inclusion of quantum ionic effects, anharmonicity, and phonon-phonon
interaction within the framework of the stochastic self-consistent harmonic
approximation (SSCHA) leads to an increase in the critical dynamical pressure
to $20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ as compared to
$5\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ in a purely harmonic treatment of
classical nuclei. We identify the change in the crystal structure due to
quantum ionic effects to be the main driving force for this increase and
demonstrate that this can already be understood at the harmonic level by
considering zero-point energy corrections to the total electronic energy of
the system. Importantly, the previously determined critical pressure of
kinetic stability $p_{\text{kin}}=$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$
still poses a stricter bound for the synthesizability of BaSiH8 and similar
hydride materials than the dynamical stability and therefore constitutes a
more rigorous and accurate estimate for the experimental realizability of
these structures. With respect to our methodology, we demonstrate the use of
highly accurate machine-learned interatomic potentials to access considerably
larger population sizes and supercell dimensions, enabling carefully
converging and speeding up the SSCHA calculations.

##  I Introduction

The discovery of high-temperature superconductivity in H3S at extreme
pressures [1] stimulated an intense hunt for novel hydride compounds with even
higher $T_{\text{c}}$’s, spearheaded by computational material discovery [2,
3, 4, 5, 6]. One prominent example is LaH10, which has been shown to
superconduct up to temperatures of $265\text{\,}\mathrm{K}$ at pressures of
$\sim$$190\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ [7, 8]. While it is very
tempting to continue searching for materials with record-breaking
$T_{\text{c}}$’s [9, 10, 11], lowering the required stabilization pressures is
even more important in view of technological applications [12, 13, 14, 15, 16,
17, 18, 19, 11, 20].

In a recent paper [15] some of us proposed a strategy to bring the
stabilization pressures of high-$T_{\text{c}}$ hydrides closer to ambient
pressure, based on the concept of an optimized chemical precompression. In
fact, we identified LaBH8, a hydride superconductor with a $T_{\text{c}}$ $>$
$100\text{\,}\mathrm{K}$, dynamically stable at an unprecedentedly low
pressure. In a follow-up work, we showed that other hydride superconductors
with the same $Fm\bar{3}m$ $XY$H8 structural template can be identified with
even lower critical pressures of stability, such as SrSiH8 and BaSiH8 [19].
The latter is particularly interesting, as it remains dynamically stable down
to $3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$. We note in passing that these
estimates of dynamical stability were based on anharmonic frozen-phonon
calculations.

All of the mentioned $Fm\bar{3}m$ $XY$H8 compounds are thermodynamically
stable only at pressures above
$\sim$$100\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, and metastable below. A
conceivable route to synthesize these materials would hence be to obtain them
at high pressures where they are thermodynamically stable, and quench them to
lower pressures.

The standard criterion employed in literature to estimate how far a metastable
phase can be quenched down in pressure is dynamical (phonon) stability.
However, dynamical stability indicates only that a structure is in a local
minimum of the potential energy surface. To estimate its actual lifetime
(kinetic stability) one needs also to estimate the height of the barriers that
separate the current minimum from other minima. In Ref. [19], some of us
introduced a rigorous method to assess the kinetic stability pressure
$p_{\text{kin}}$ by explicitly calculating the energy barrier protecting the
metastable $Fm\bar{3}m$ structure from decomposition as a function of
pressure, using the variable-cell nudged elastic band method [21]. For BaSiH8,
for example, we found a $p_{\text{kin}}$ of
$\sim$$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, significantly higher than
the dynamical value $p_{\text{dyn}}$ =
$3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$.

It was argued in a recent work [22] that quantum lattice effects treated
within the stochastic self-consistent harmonic approximation (SSCHA)
drastically increase the dynamical stabilization pressure $p_{\text{dyn}}$ for
LaBH8 and it was further suggested that a similar increase in $p_{\text{dyn}}$
should be expected for other $Fm\bar{3}m$ $XY$H8 hydrides.

To investigate this, we apply the SSCHA formalism to BaSiH8, which, so far,
has the lowest $p_{\text{dyn}}$ among all $Fm\bar{3}m$ $XY$H8 hydrides. In
SSCHA, a major bottleneck is represented by the need to use large supercells
and large numbers of individuals if one wants to fully converge the
calculation. We overcome this problem by employing machine-learned moment
tensor potentials (MTP) [23, 24] that allow us to obtain total energies,
forces, and stresses with DFT accuracy but at a fraction of the computational
cost [25, 26, 27]. To our knowledge, this work represents the first
combination of MTPs with the SSCHA method. In addition, we introduce a method
to discern the contribution of quantum ionic (QI) effects from those of
anharmonic and phonon-phonon (ph-ph) effects.

We find that $p_{\text{dyn}}$ increases from
$3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ to about
$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ within the SSCHA and that this
rise can almost entirely be attributed to QI effects, with actual anharmonic
and ph-ph effects playing only a subordinate role. In fact, the same crystal
structure that minimizes the free energy within SSCHA can already be obtained
in density-functional theory (DFT) by including zero-point energies (ZPE).

Importantly, we demonstrate that even after including QI, anharmonic, and ph-
ph effects within the framework of SSCHA, the actual limit of stability is
still set by $p_{\text{kin}}$
($\sim$$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$), as stated in our previous
work [19].

##  II Results

###  II.1 Ab-initio machine-learned interatomic potentials

In the self-consistent harmonic approximation, the full anharmonic and
interacting lattice is mapped onto an auxiliary harmonic system and the free
energy $\mathcal{F}$ of the full system is approximated by the minimum of the
free energy of the auxiliary harmonic system [28, 29, 30, 31]. In the SSCHA,
this minimization is performed stochastically via Monte-Carlo summation and
importance sampling over several consecutive ensembles (populations) of a
large number of individuals. Each individual here corresponds to a supercell
structure with displaced atomic positions, where the supercell size determines
the density of phonon wave vectors in the Brillouin zone [32]. More details
are provided in the Method section, the Supplementary Material (SM) 111The
Supplemental Material is available on the journal website., and in Refs. [34,
35, 36, 37, 38, 39, 40].

In practice, to calculate accurate phonon frequencies within SSCHA, in
particular for slow-converging soft modes, one needs to consider population
sizes of several ten or hundred thousands individuals. In addition, one also
needs to converge the supercell size.

Doing this fully at a DFT level is computationally prohibitive, which is why
we made use of MTPs in this work. For every pressure, MTPs were trained on DFT
results of 50 structures randomly chosen out of the SSCHA random-displacement
individuals in 2$\times$2$\times$2 supercells. We then validated the trained
MTPs for all other individuals by comparing the total energies, forces and
stress components. This validation is shown in Fig. 1, demonstrating the
exceptional accuracy of the used MTPs (see Supplementary Fig. 2 in the SM for
other pressures, as well as forces and stresses). As can be appreciated in
this figure, the root-mean-squared error (RMSE) is below
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$,
i.e., at the same level as the error in DFT. The inset also shows that the
potential energy surface of the slow-converging, $T_{2g}$ mode at $\Gamma$ is
reproduced very nicely with the MTPs.

As a final validation, we compare the SSCHA phonon dispersions obtained using
only DFT with those employing only MTPs (see Supplementary Fig. 3 in the SM)
and find very good agreement, with only minor differences in the $T_{2g}$ mode
at $\Gamma$ and the $E_{g}$ at $X$. To fully converge these modes within
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, we increased the populations sizes
within MTP-SSCHA up to $100\,000\text{\,}$ individuals compared to
$10\,000\text{\,}$ for the DFT-SSCHA calculations.

Figure 1: MTP validation: MTP total energy ($E_{\mathrm{MTP}}$) versus DFT
total energy ($E_{\mathrm{DFT}}$) for all $19\,250\text{\,}$ individuals of a
SSCHA calculation for a lattice constant
$a=$6.242\text{\,}\mathrm{\text{\AA}}$$ (blue scatter plot).
$E_{\mathrm{ref}}$ is the DFT total energy of the high-symmetry structure with
undisplaced H positions. The diagonal black line serves as guide to the eye.
The inset shows the full frozen-phonon potential of the lowest $T_{2g}$ mode
at $\Gamma$ obtained with DFT (solid black line) and MTP (dotted blue line),
as well as the harmonic potential (dashed grey line).

The use of MTPs does not only substantially speed up the calculations 222We
found MTPs to be a factor of about $20\,000\text{\,}$ faster than DFT in the
case of 2$\times$2$\times$2 supercells of BaSiH8., but also gives access to
larger supercells. In this work, we performed additional SSCHA calculations
using MTPs on $n\times n\times n$ supercells with $n=1,2,3,4$ at all studied
pressures. The convergence of the free energy, the structural parameters, and
the phonon dispersions with respect to the supercell size is provided in
Supplementary Figs. 4-6 in the SM. An overview of all performed SSCHA runs is
given in Supplementary Tab. 1 in the SM. Unless stated otherwise, all SSCHA
results presented in the following have been obtained with MTPs for
$100\,000\text{\,}$ individuals in 4$\times$4$\times$4 supercells.

###  II.2 Structural parameters and electronic dispersion

The $Fm\bar{3}m$ phase of BaSiH8 has a face-centered cubic unit cell with 10
atoms in the primitive cell, where Ba and Si occupy Wyckoff $4a/b$ sites and
the H occupy $32f$ sites. The eight H atoms form rhombicuboctahedral cages
around the Ba atoms and cubic cages around the Si atoms. The structure has
only two free parameters, namely the lattice constant $a$ and the Wyckoff
coordinate of the $32f$ sites $x$, defining the H-H distance
$d_{\text{H-H}}=2a\cdot x$ (side length of the cubic cage) and the H-Si
distance $d_{\text{H-Si}}=a\sqrt{3}\cdot x$ (half the space diagonal of the
cubic cage).

Relaxing the structure within DFT to target pressures of 10, 20, 25, and
$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, we obtained lattice constants
between $6.5\text{\,}\mathrm{\text{\AA}}$ and
$6.2\text{\,}\mathrm{\text{\AA}}$, and H-Si distances of about
$1.6\text{\,}\mathrm{\text{\AA}}$, as shown in Table 1. An extensive list of
the structural parameters from ambient pressure up to
$100\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, as well as the fit to the Birch-
Murnaghan equation of state can be found in Supplementary Note 1 in the SM.

Starting from the atomic positions obtained in DFT and the harmonic dynamical
matrices obtained in density-functional perturbation theory (DFPT)
calculations at each pressure, we performed constant-volume SSCHA relaxation
calculations 333Within SSCHA, the stress tensor is obtained as
$P_{\alpha\beta}=-V_{\text{sc}}^{-1}\left(\partial\mathcal{F}/\partial\varepsilon_{\alpha\beta}\right)|_{\varepsilon=0}$,
where $V_{\text{sc}}$ is the supercell volume, $\mathcal{F}$ the free energy
functional, and $\varepsilon_{\alpha\beta}$ the strain tensor [37]. The
pressure is then obtained as $\tilde{p}=\sum_{\alpha}P_{\alpha\alpha}/3$.. The
corresponding parameters, indicated by $\tilde{x},\tilde{d}$, and $\tilde{p}$,
are reported in Table 1.

Table 1: Structural parameters: Lattice constant $a$, Wyckoff parameter $x$, H-Si distance $d_{\text{H-Si}}$, and pressure $p$ after the relaxation with respect to the DFT total energy and after the constant-volume relaxation within SSCHA ($\tilde{x},\tilde{d}$, and $\tilde{p}$). $a\,/\,\text{\AA}$ | $x$ | $d_{\text{H-Si}}\,/\,\text{\AA}$ | $p\,/\,\text{GPa}$ | $\tilde{x}$ | $\tilde{d}_{\text{H-Si}}\,/\,\text{\AA}$ | $\tilde{p}\,/\,\text{GPa}$
---|---|---|---|---|---|---
6.541 | 0.1434 | 1.625 | 10 | 0.1459 | 1.653 | 12.2
6.323 | 0.1471 | 1.611 | 20 | 0.1498 | 1.640 | 22.6
6.242 | 0.1483 | 1.603 | 25 | 0.1510 | 1.633 | 27.9
6.171 | 0.1494 | 1.597 | 30 | 0.1521 | 1.626 | 33.1

We observe an elongation of $d_{\text{H-H/Si}}$ of about
$30\text{\,}\mathrm{m\text{\AA}}$ (2%) for all pressures and an increase in
pressure of about 2 to $3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, i.e. $\sim
20\%$ at $10\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ and $\sim 10\%$ at
$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$. The change in atomic positions
introduces only small changes in the electronic structure, as demonstrated in
Fig. 2, where we compare the electronic bands and densities of states (DOS)
for $x$ and $\tilde{x}$. The largest differences are found above and below the
Fermi energy, whereas electronic bands and DOS at the Fermi energy, and hence
the Fermi surface, remain essentially unchanged.

Figure 2: Difference in electronic properties: a electronic bands, and b
density of states for the structure with H positions defined by $x$ (DFT
minimum, blue line) and $\tilde{x}$ (SSCHA minimum, red line) for
$a=$6.242\text{\,}\mathrm{\text{\AA}}$$.

###  II.3 Phonon dispersions and lattice instability

Moving on, we evaluate and compare the DFPT and SSCHA phonon dispersions at
all studied pressures, as shown in Fig. 3. Similar to the results for LaBH8 in
Ref. [22], we find that the high-energy optical modes are strongly
renormalized to lower frequencies. In particular, a significant softening
occurs for the threefold degenerate $T_{2g}$ mode at $\Gamma$ (harmonic values
around $50\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ in Fig. 3), which becomes
imaginary and indicates a (dynamic) lattice instability for lattice constants
$a>$6.323\text{\,}\mathrm{\text{\AA}}$$, corresponding to
$p=$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$ and
$\tilde{p}=$22.6\text{\,}\mathrm{GPa}$$. Thus, the inclusion of quantum
lattice effects within the SSCHA shifts the dynamical stability pressure from
the anharmonic frozen-phonon value $p_{\text{dyn}}$ =
$3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ to
$\tilde{p}_{\text{dyn}}=$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$ (see
Supplementary Fig. 7 in the SM). This
$\sim$$17\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ difference is substantial,
but considerably smaller than the
$\sim$$40\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ shift reported for LaBH8
[22].

Figure 3: Harmonic and SSCHA phonon dispersions: Phonon dispersions for
various lattice constants $a$ along a high-symmetry path of the BZ. The dashed
black lines correspond to harmonic calculations and the solid red lines to the
SSCHA results. The indicated pressures correspond to the DFT (SSCHA) pressures
for relaxed atomic positions defined by $x$ ($\tilde{x}$) [43].

We want to note in passing that we also calculated the fourth order
corrections to the phonon frequencies in SSCHA and find, in contrast to Ref.
[22], but in accordance with other works employing SSCHA [40, 39, 35, 38, 36],
only minor differences to the results obtained only up to third order. The
maximum phonon energy differences are on the order of about
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ for all pressures. The phonon
dispersions in the 2$\times$2$\times$2 supercells 444Calculating fourth-order
corrections on larger supercells is computationally unfeasible with our
currently available resources. This calculation in an $n\times n\times n$
supercell requires 4D arrays with $(3\cdot n^{3}\cdot<\text{number of atoms in
uc}>)^{4}$ elements. For 10 atoms in the unit cell and $n=3$, this yields
$810^{4}$ double-precision entries, for which
$3.13\text{\,}\mathrm{T}\mathrm{B}$ of RAM are needed on a single shared-
memory node. with second (auxiliary), third, and fourth order terms for all
studied pressures are shown in Supplementary Fig. 8 in the SM.

###  II.4 Different effects contributing to frequency shifts

The observed changes in the phonon dispersions when employing SSCHA and the
resulting different dynamical stabilization pressures result from a
combination of several effects that are not included at the level of standard
DFT and DFPT. These are most importantly the vibrational contributions of the
ions to the free energy, phonon anharmonicity, and ph-ph interactions. In the
following, we will present an attempt to disentangle and determine the
importance of each of these effects for BaSiH8.

QI effects:  First, we want to look at the contributions to the total energy
originating from the so-called zero-point vibrations, i.e, vibrations of the
ions around their equilibrium positions due to the quantum mechanical
treatment of the nuclei, absent in the classical, clamped-nuclei picture [45].
In the Born-Oppenheimer approximation, the total energy
$E_{\text{tot}}[\mathbf{R}]$ (at $T=$0\text{\,}\mathrm{K}$$) for ionic
positions $\mathbf{R}$ is given by the sum of the internal electronic energy
$E_{\text{el}}[\mathbf{R}]$ and the ZPE contributions of the nuclei
$E_{\text{ZP}}[\mathbf{R}]$. In most solids, $E_{\text{ZP}}$ is much smaller
than $E_{\text{el}}$ and can be safely neglected. However, due to the small
mass of H and the resulting high phonon frequencies in hydrides, the ZPE can
become substantial and thus cause a modification of the equilibrium crystal
structure.

At the harmonic level, the true ZPE can be approximated via
$E_{\text{ZP}}[\mathbf{R}]\approx\int_{0}^{\infty}\mathrm{d}\omega\rho_{\mathbf{R}}(\omega)\hbar\omega/2$,
where $\rho_{\mathbf{R}}(\omega)$ is the DFPT phonon density of states and
$\hbar\omega/2$ the ZPE of a quantum harmonic oscillator. At constant volume,
the only free parameter in the $Fm\bar{3}m$ structure is the Wyckoff parameter
$x$, for which we have plotted $E_{\text{tot}}$, $E_{\text{el}}$, and
$E^{\text{harm}}_{\text{ZP}}$ in Fig. 4 for a lattice constant of
$a=$6.242\text{\,}\mathrm{\text{\AA}}$$. The results for other lattice
constants, i.e., pressures, are provided in Supplementary Fig. 9 in the SM.

As can be appreciated in this figure, the inclusion of the ZPE, even at the
harmonic level, shifts the position of the minimum of the total energy
considerably and puts it almost exactly at the minimum position $\tilde{x}$
predicted by the SSCHA. The differences in $d_{\text{H-Si}}$ between the SSCHA
calculations and the ZPE analysis are, in fact, of the order of
$1\text{\,}\mathrm{m}\mathrm{\text{\AA}}$, i.e., well within the observed
stochastic noise in SSCHA. We want to note that the same is true for LaBH8
(see Supplementary Fig. 10 of the SM). Furthermore, inclusion of ZPE reduces
the total energy at its minimum by
$\sim$$30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{u}\mathrm{c}$,
agreeing very nicely with the result from SSCHA
($\sim$$27\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{u}\mathrm{c}$).

This is a quite remarkable result as it demonstrates the importance of QI
effects of the light H ions on the dynamic stability of the hydride materials,
and shows that the minimum structure from SSCHA can already be obtained at the
level of harmonic ZPE corrections, at least for this class of materials.

Figure 4: Electronic total energy, ZPE and resulting total energy as a
function of H-Si distance, where the DFT minimum $x$ and the SSCHA minimum
$\tilde{x}$ are marked explicitly, the latter coinciding with the minimum
position of $E_{\text{tot}}$. The three energy curves are plotted relative to
their respective values at $x$. The solid lines represent a cubic spline for
$E_{\text{el}}$ and second order polynomial fits for the other energies.
Figure 5: Different effects contributing to shifts in the phonon dispersions:
a Harmonic dispersion for atomic positions defined by $x$ (dashed black) and
$\tilde{x}$ (dotted blue). b SSCHA-obtained dispersions for fixed atomic
positions defined by $x$ (dotted green lines) and harmonic reference (dashed
black). c Dispersions for SSCHA-relaxed atomic positions defined by
$\tilde{x}$ (solid red), obtained from $C_{\text{sum}}$ (dotted grey, see text
for more details), and for the harmonic case (dashed black).

Having established that the ZPE has a crucial effect on the structure, we
investigate the effect of the changed structure on the phonon dispersions. In
Fig. 5a, we present the harmonic dispersions for atomic positions defined by
$x$ and $\tilde{x}$. We observe large differences for the high-energy optical
modes above $150\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$, but also for the low
$T_{2g}$ mode at $\Gamma$. The energy shifts for these modes are between 15
and $25\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$.

Anharmonicity and ph-ph interaction effects: Having established the QI effects
on the structure and the phonon dispersions, we now want to assess the
contributions of phonon anharmonicity (anh) and ph-ph interactions. To do
that, we perform a SSCHA calculation while keeping the ions fixed at the DFT
equilibrium positions, thus qualitatively 555In the case of fixed ions, we
obtain non-vanishing forces within SSCHA. As the individual force components
on the H atoms are still small
($\sim$$150\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{\AA}}$
in the cell with $a=$6.242\text{\,}\mathrm{\text{\AA}}$$, for example), a
fixed-ion calculation seems applicable for a qualitative insight. The forces
on atom $i$ are obtained as
$\mathbf{F}_{i}=-\partial\mathcal{F}/\partial\mathcal{\mathbf{R}}_{i}$.
removing structural effects on the phonon frequencies. The phonon dispersions
obtained from this calculation are presented in Fig. 5b, where we find that
all H-dominated optical modes in the whole BZ experience a sizable frequency
renormalization. It is worth noting that fixed-ions SSCHA calculations
indicate the onset of dynamical instability just between 5 and
$10\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ (8.4 and
$12.2\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, respectively, for the SSCHA
pressure $\tilde{p}$, see Supplementary Fig. 11 in the SM), which is
reasonably close to the frozen-phonon anharmonic (harmonic) $p_{\text{dyn}}$ =
$3(5)\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ [19].

Combining QI, anh, and ph-ph effects: As we made the attempt of separating the
different contributions from QI effects, anharmonicity, and ph-ph
interactions, it is tempting to combine the individual contributions to the
phonon dispersion as a simple sum and compare it to the full SSCHA
calculation. We approach this task in real space via adding the force
constants $C_{\text{sum}}=C_{\text{harm}}+\Delta C_{\text{QI}}+\Delta
C_{\text{anh/phph}}$, where $C_{\text{harm}}$ is the force constant matrix for
harmonic phonons of structure $x$, $\Delta C_{\text{QI}}$ are the force
constant contributions due to the structural changes based on including the
ZPE (QI effects), and $\Delta C_{\text{anh/phph}}$ are the force constant
contributions to the anharmonic and ph-ph interaction effects (see
Supplementary Note 3 in the SM).

The dispersions obtained from $C_{\text{sum}}$ are shown in Fig. 5c, where we
also present as reference the phonon dispersions from a full SSCHA
calculation. Again, we find very good agreement between these results,
providing further, a posteriori justification and support for the
qualitatively introduced separation Ansatz. In Table 2, we summarize the
considered effects, methods, and the corresponding physical description of the
ions.

Table 2: Overview of ionic treatment: The separate cases are classified according to the structural and phonon treatment. The phonons are obtained either via DFPT or SSCHA. The ground-state (GS) structure is determined by minimizing either the electronic energy within DFT or the total energy including the ZPE (using DFPT or SSCHA.) ${}^{\text{\hskip 21.68231ptphonon}\rightarrow}_{\downarrow\text{structure}}$ | DFPT | SSCHA
---|---|---
$x$ $\hat{=}$ min($E_{\text{el}}$) | Classical ions | Interact. quantum ions
DFT | in el. GS | in el. GS
| standard |  anh+ph-ph
$\tilde{x}$ $\hat{=}$ min($E_{\text{tot}}$) | Quantum ions | Interact. quantum ions
ZPE/SSCHA | in lattice GS | in lattice GS
| QI |  QI+anh+ph-ph

###  II.5 Superconductivity

As BaSiH8 is potentially a very promising high-$T_{\text{c}}$ superconductor,
we also want to assess the implications of the above mentioned effects on its
superconducting (SC) properties. To do that, we solved the anisotropic Migdal-
Eliashberg (ME) equations as implemented in EPW [47] for the four cases in
Tab. 2. Details about the calculation within EPW are provided in the Method
section and in the SM, at this point we only want to highlight that for each
case we used the corresponding force constants to compute the dynamical
matrices, and computed the electron-phonon ($ep$) coupling matrix elements as
the self-consistent first-order variation of the potential using the
equilibrium positions as defined in Tab. 2 666A more rigorous treatment would
require the replacement of the polarization vectors
$\mathbf{e}_{\mathbf{q}\nu}$, as well as either the inclusion of the force
term arising in the DF(P)T description of the structure defined by
$\tilde{x}$, or an $\mathcal{F}$-based description of the $ep$ matrix elements
within the SSCHA.. In Table 3, we summarize the obtained values for quantities
characterizing the SC state, i.e., the $ep$ coupling strength $\lambda$, the
logarithmic average of the phonon frequencies $\omega_{\text{log}}$, and the
SC critical temperature $T_{\text{c}}$ 777The $ep$ coupling strength
$\lambda(\omega)=2\int_{0}^{\omega}\mathrm{d}\omega^{\prime}\frac{\alpha^{2}F(\omega^{\prime})}{\omega^{\prime}}$
and the logarithmic average phonon frequency
$\omega_{\text{log}}=\exp\left(\frac{2}{\lambda}\int_{0}^{\infty}\mathrm{d}\omega\frac{\alpha^{2}F(\omega)\ln\omega}{\omega}\right)$
are obtained from the Eliashberg spectral function $\alpha^{2}F(\omega)$. The
corresponding Eliashberg spectral functions $\alpha^{2}F(\omega)$ and the
cumulative coupling strengths $\lambda(\omega)$ are shown in Supplementary
Fig. 12 in the SM. We want to stress that the provided values for
$T_{\text{c}}$ are obtained by the solution of the full ME equations. The
distribution of the SC gap function $\Delta_{\mathbf{k}}$ indicates no change
in the distinct two-gap shape calculated in Ref. [19], and is therefore not
shown here.

Table 3: SC properties from ME theory: Critical temperature $T_{\text{c}}$, $ep$ coupling strength $\lambda$, and logarithmic phonon frequency average $\omega_{\text{log}}$ for the cases discussed in the text. effect | (struct., phon.) | $\omega_{\mathrm{log}}\,/\,\mathrm{meV}$ | $\lambda$ |  $T_{\text{c}}$ / K
---|---|---|---|---
standard | ($x$, harm.) | 54 | 1.25 | 84
QI | ($\tilde{x}$, harm.) | 47 | 1.43 | 82
anh+ph-ph | ($x$, SSCHA) | 54 | 1.38 | 90
QI+anh+ph-ph | ($\tilde{x}$, SSCHA) | 28 | 2.44 | 96

The differences in $\omega_{\text{log}}$, $\lambda$, and $T_{\text{c}}$ are in
the order of 10-15% except for the full SSCHA calculation, where we see a
considerable increase in $\lambda$ to almost double the harmonic value, but
also a decrease in $\omega_{\text{log}}$, compensating the enhancement of
$\lambda$. The resulting $T_{\text{c}}$ is increased from 84 to
$96\text{\,}\mathrm{K}$, showing that the full inclusion of all discussed
effects results only in $\sim$15% change in $T_{\text{c}}$ for BaSiH8.

##  III Discussion

In this work, we study the effects of quantum lattice dynamics within the
SSCHA framework on the structure and the dynamical stability of the
$Fm\bar{3}m$ phase of BaSiH8. The SSCHA structure relaxation suggests a 2%
elongation of the H-H and H-Si bonds for the studied pressure range of 10 to
$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$
($\sim$30\text{\,}\mathrm{m\text{\AA}}$$).

In the phonon dispersions, we find an overall softening of the high optical
modes, as well as a dynamic lattice instability characterized by imaginary
SSCHA phonon frequencies in the $T_{2g}$ mode at $\Gamma$ below
$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, setting the estimate for the
critical dynamical pressure to
$\tilde{p}_{\text{dyn}}\approx$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$. We
have further demonstrated the importance of QI effects over anharmonicity and
ph-ph interactions, and found that the change in structure, and consequently
in pressure, can already be understood by considering harmonic ZPE corrections
to the total electronic energy of the system alone (which can be obtained much
faster and easier than performing a full SSCHA calculation).

We are now left with the question: what is the stability boundary of
$Fm\bar{3}m$-BaSiH8? In our previous work on BaSiH8 [19], we challenged the
common practice of assuming the range of metastability of high-pressure
hydride phases to coincide with the range of (an)harmonic dynamical stability,
which systematically underestimates the stabilization pressures needed to
synthesize these materials in reality [50, 51, 52, 53]. Dynamical stability is
only a prerequisite for thermodynamic metastability, which is characterized by
the existence of a distinctive enthalpy barrier that protects a metastable
phase from decomposition into other phases (kinetic stability). We calculated
the enthalpy transition path to the thermodynamic groundstate at different
pressures (corresponding to a decomposition of the $Fm\bar{3}m$ BaSiH8 phase
into BaSiH6 \+ H2 in molecular form), and could estimate the barrier height
from the intermediate structures. In combination with the calculated convex
hulls for the B-S-H system, we can argue with confidence that the $Fm\bar{3}m$
BaSiH8 phase could be synthesized above
$100\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$, and retained down to
$\sim$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$, where a distinctive
enthalpy barrier still exists. At lower pressures, metastable $Fm\bar{3}m$
BaSiH8 will decompose, even though (anharmonic) lattice dynamics calculations
predict it to be stable. Hence, kinetic stability poses a stricter bound for
synthesizability than dynamical stability.

In conclusion, employing ab-initio machine-learned MTPs, we were able to
perform SSCHA calculations for BaSiH8 at various pressures for supercells up
to $4\times 4\times 4$ and more than $100\,000\text{\,}$ individuals. The
inclusion of QI effects, anharmonicity, and ph-ph interactions within the
SSCHA-framework increases the pressure of dynamical stability from
$p_{\text{dyn}}\approx$3\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$ to
$\tilde{p}_{\text{dyn}}\approx$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$. We
identified the change in structure due to QI effects to be the main driving
force here, something that can already be captured to good approximation at
the level of harmonic zero-point energy corrections.

Most importantly, the determined
$\tilde{p}_{\text{dyn}}\approx$20\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$ is
still below $p_{\text{kin}}\approx$30\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$$
posed by the concept of kinetic stability, thus the latter represents a much
stricter bound for the stability and realizability in these materials.

##  IV Methods

###  IV.1 DF(P)T calculations

All DFT and DFPT calculations of electronic and vibrational properties were
carried out using the plane-wave pseudopotential code Quantum Espresso [54],
scalar-relativistic optimized norm-conserving Vanderbilt pseudopotentials
[55], and the Pbe-Gga exchange and correlation functional [56]. The unit cell
calculations are done in the fcc primitive unit cell with 10 atoms, a
12$\times$12$\times$12 $\mathbf{k}$-grid, and a plane-wave cutoff energy of
$80\text{\,}\mathrm{R}\mathrm{y}$. The 2$\times$2$\times$2 supercell
calculations were done on a 6$\times$6$\times$6 $\mathbf{k}$-grid. Further
details are provided in Supplementary Method 1 in the SM.

###  IV.2 SSCHA calculations

The calculations in the SSCHA are done in the constant-volume relaxation mode,
i.e. minimizing the free energy with respect to the average atomic positions
$\mathcal{R}$ and the force constants $\Phi$, as implemented in the SSCHA
python package [40]. We use the DFT equilibrium atomic positions and the DFPT
dynamical matrices on a 2$\times$2$\times$2 $\mathbf{q}$-grid as initial
guesses for $\mathcal{R}$ and $\Phi$, respectively. The starting point for the
larger supercells is obtained by interpolating the previously converged
auxiliary dynamical matrices.

The total energies, forces, and stress tensors for the individuals are
obtained from DFT calculations or from machine-learned interatomic potentials
in the framework of MTPs, see below. At the end of a minimization run, a new
population with higher number of individuals $N$ is generated from the
minimized trial density matrix until convergence. We set two stopping criteria
for the minimization loops: a Kong-Liu ratio for the effective sample size
888The effective sample size is calculated as
$N_{\text{eff}}=(\sum_{i}\rho_{i}^{2})/(\sum_{i}\rho_{i})^{2}$, where the
$\rho_{i}$ are the importance sampling weighting factors. of 0.2, and ratio of
$<10^{-7}$ between the free energy gradient with respect to the auxiliary
dynamical matrix and its stochastic error. In calculations based on DFT we
increased $N$ up to $10^{4}$ individuals, for the MTP cases up to $10^{5}$.

The anharmonic phonon dispersions are obtained from the positional free-energy
Hessians without the fourth-order term, if not specified otherwise explicitly.
The final atomic positions are obtained from the converged average atomic
positions $\mathcal{R}$ and the pressure as derivative of the converged free
energy with respect to a strain tensor. The free energy difference between the
last two populations in the 2$\times$2$\times$2 cells is well below
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{u}\,\mathrm{c}\,$,
and well below
$0.1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{u}\,\mathrm{c}\,$
for higher cells. The total forces in the last population are well below
$10^{-6}\,\mathrm{meV}\text{\AA}^{-1}$. The physical phonon frequency
differences between the last two populations are below
$5\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ for the DFT cases and well below
$1\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ for the MTP case ($T_{2g}$ and
$A_{2u}$ converge slower, see SM). All calculations were carried out at zero
temperature. We note in passing that with these settings, we could reproduce
all LaBH8 results from Ref. [22].

Further details are provided in Supplementary Method 2 in the SM.

###  IV.3 Moment tensor potentials

The MTPs were trained and evaluated using the MLIP package [23, 58]. We choose
a functional form of level 26, eight radial basis functions,
$R_{\text{cut}}=$5.0\text{\,}\mathrm{\text{\AA}}$$ and
$R_{\text{min}}=$1.2\text{\,}\mathrm{\text{\AA}}$$, and trained on 50
structures in a 2$\times$2$\times$2 supercell randomly chosen out of all
individuals of the DFT SSCHA calculations for a given pressure. We validated
the potentials on all individuals in the DFT SSCHA calculations and find a
RMSE on the total energy of
0.5-$0.6\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{m}$,
45-$50\text{\,}\mathrm{meV}\text{\,}{\mathrm{\text{\AA}}}^{-1}$ for the force
components, and 0.3-$0.4\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ for the
diagonal stress tensor components. We further validated the MTPs on 30
randomly chosen individuals in a 3$\times$3$\times$3 supercell and achieve
similar RMSEs. The validations and RMSEs for each pressure are shown in
Supplementary Fig. 2 and Supplementary Note 2 in the SM.

###  IV.4 ZPE and total energy

The internal electronic energy $E_{\mathrm{el}}[\mathbf{R}(x)]$ is obtained
from DFT calculations at fixed volume by varying the Wyckoff parameter $x$ of
the H positions. The explicit positions $\mathbf{R}(x)$ are given in the SM.
The phonon density of states $\rho_{\mathbf{R}(x)}(\omega)$ is obtained using
DFPT on a 2$\times$2$\times$2 $\mathbf{q}$-grid, interpolated on a
16$\times$16$\times$16 $\mathbf{q}$-grid. Smooth ZPE and total energy curves
are obtained by second-order polynomial fits in $x$. Due to the shift out of
the DFT equilibrium structure, forces on the individual H atoms arise at the
DFT level. Around the total energy minimum, the force components are in the
order of $150\text{\,}\mathrm{meV}\text{\,}{\mathrm{\text{\AA}}}^{-1}$, i.e.,
small enough to warrant the use of linear-response theory to gain qualitative
and systematic insights. DFT diagonal stress tensor components (pressures) are
decreased by about $2\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$.

###  IV.5 Migdal-Eliashberg theory

The Wannier interpolation of the $ep$ matrix elements onto dense
$\mathbf{k}$\- and $\mathbf{q}$-grids, and the subsequent self-consistent
solution of the fully anisotropic Migdal-Eliashberg equations were done in Epw
[59, 47], for all the cases in Tab. 3. We used coarse 6$\times$6$\times$6 and
fine 30$\times$30$\times$30 $\mathbf{k}$\- and $\mathbf{q}$ grids, a Matsubara
frequency cutoff of $1\text{\,}\mathrm{e}\mathrm{V}$, and a Morel-Anderson
pseudopotential $\mu^{*}=0.10$. Further details are provided in Supplementary
Method 3 in the SM.

## Acknowledgments

We thank Pedro Pires Ferreira, Carla Verdi, and Luigi Ranalli for fruitful
discussions. This work was supported by the Austrian Science Fund (FWF),
projects P30269 and P32144. Calculations have been performed at the dCluster
of TU Graz, as well as on the Vienna Scientific Cluster (VSC4 and VSC5). L.B.
acknowledges support from Fondo Ateneo-Sapienza 2018-2022. S.D.C. acknowledges
computational resources from CINECA, proj. IsC90-HTS-TECH_C and IsC99-ACME-C.

## Author contributions

R.L., E.K., and S.D.C. performed the calculations, M.A. introduced the idea of
MTP, and C.H. and L.B. conceived and supervised the project. All authors
contributed to the discussion of the results and participated in preparing the
manuscript.

## Competing interests

The authors declare no competing interests.

## Data availability

The authors confirm that the data supporting the findings of this study are
available within the article and its supplementary materials. Further
information is available upon request.

## References

  * Drozdov  _et al._ [2015] A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system, Nature 525, 73 (2015).
  * Boeri and Bachelet [2019] L. Boeri and G. B. Bachelet, Viewpoint: the road to room-temperature conventional superconductivity, Journal of Physics: Condensed Matter 31, 234002 (2019).
  * Lilia  _et al._ [2022] B. Lilia, R. Hennig, P. Hirschfeld, G. Profeta, A. Sanna, E. Zurek, W. E. Pickett, M. Amsler, R. Dias, M. I. Eremets, C. Heil, R. J. Hemley, H. Liu, Y. Ma, C. Pierleoni, A. N. Kolmogorov, N. Rybin, D. Novoselov, V. Anisimov, A. R. Oganov, C. J. Pickard, T. Bi, R. Arita, I. Errea, C. Pellegrini, R. Requist, E. K. U. Gross, E. R. Margine, S. R. Xie, Y. Quan, A. Hire, L. Fanfarillo, G. R. Stewart, J. J. Hamlin, V. Stanev, R. S. Gonnelli, E. Piatti, D. Romanin, D. Daghero, and R. Valenti, The 2021 room-temperature superconductivity roadmap, Journal of Physics: Condensed Matter 34, 183002 (2022).
  * Flores-Livas  _et al._ [2020] J. A. Flores-Livas, L. Boeri, A. Sanna, G. Profeta, R. Arita, and M. Eremets, A perspective on conventional high-temperature superconductors at high pressure: Methods and materials, Physics Reports 856, 1 (2020).
  * Hilleke  _et al._ [2022] K. P. Hilleke, T. Bi, and E. Zurek, Materials under high pressure: a chemical perspective, Applied Physics A 128, 441 (2022).
  * Hilleke and Zurek [2022] K. P. Hilleke and E. Zurek, Tuning chemical precompression: Theoretical design and crystal chemistry of novel hydrides in the quest for warm and light superconductivity at ambient pressures, Journal of Applied Physics 131, 070901 (2022).
  * Drozdov  _et al._ [2019] A. P. Drozdov, P. P. Kong, V. S. Minkov, S. P. Besedin, M. A. Kuzovnikov, S. Mozaffari, L. Balicas, F. F. Balakirev, D. E. Graf, V. B. Prakapenka, E. Greenberg, D. A. Knyazev, M. Tkacz, and M. I. Eremets, Superconductivity at 250 K in lanthanum hydride under high pressures, Nature 569, 528 (2019).
  * Somayazulu  _et al._ [2019] M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe, M. Baldini, Y. Meng, V. V. Struzhkin, and R. J. Hemley, Evidence for Superconductivity above 260 K in Lanthanum Superhydride at Megabar Pressures, Phys. Rev. Lett. 122, 027001 (2019).
  * Peng  _et al._ [2017] F. Peng, Y. Sun, C. J. Pickard, R. J. Needs, Q. Wu, and Y. Ma, Hydrogen Clathrate Structures in Rare Earth Hydrides at High Pressures: Possible Route to Room-Temperature Superconductivity, Phys. Rev. Lett. 119, 107001 (2017).
  * Grockowiak  _et al._ [2022] A. D. Grockowiak, M. Ahart, T. Helm, W. A. Coniglio, R. Kumar, K. Glazyrin, G. Garbarino, Y. Meng, M. Oliff, V. Williams, N. W. Ashcroft, R. J. Hemley, M. Somayazulu, and S. W. Tozer, Hot Hydride Superconductivity Above 550 K, Frontiers in Electronic Materials 2, 837651 (2022).
  * Di Cataldo  _et al._ [2022] S. Di Cataldo, W. von der Linden, and L. Boeri, First-principles search of hot superconductivity in La-$X$-H ternary hydrides, npj Computational Materials 8, 2 (2022).
  * Pickard  _et al._ [2020] C. J. Pickard, I. Errea, and M. I. Eremets, Superconducting Hydrides Under Pressure, Annual Review of Condensed Matter Physics 11, 57 (2020).
  * Lv  _et al._ [2020] J. Lv, Y. Sun, H. Liu, and Y. Ma, Theory-orientated discovery of high-temperature superconductors in superhydrides stabilized under high pressure, Matter and Radiation at Extremes 5, 068101 (2020).
  * Di Cataldo  _et al._ [2020] S. Di Cataldo, W. Von Der Linden, and L. Boeri, Phase diagram and superconductivity of calcium borohyrides at extreme pressures, Phys. Rev. B 102, 014516 (2020).
  * Di Cataldo  _et al._ [2021] S. Di Cataldo, C. Heil, W. von der Linden, and L. Boeri, LaBH8: Towards high-${T}_{c}$ low-pressure superconductivity in ternary superhydrides, Phys. Rev. B 104, L020511 (2021).
  * Shipley  _et al._ [2021] A. M. Shipley, M. J. Hutcheon, R. J. Needs, and C. J. Pickard, High-throughput discovery of high-temperature conventional superconductors, Phys. Rev. B 104, 054501 (2021).
  * Zhang  _et al._ [2022] Z. Zhang, T. Cui, M. J. Hutcheon, A. M. Shipley, H. Song, M. Du, V. Z. Kresin, D. Duan, C. J. Pickard, and Y. Yao, Design Principles for High-Temperature Superconductors with a Hydrogen-Based Alloy Backbone at Moderate Pressure, Physical Review Letters 128, 047001 (2022).
  * Liang  _et al._ [2021] X. Liang, A. Bergara, X. Wei, X. Song, L. Wang, R. Sun, H. Liu, R. J. Hemley, L. Wang, G. Gao, and Y. Tian, Prediction of high-${T}_{c}$ superconductivity in ternary lanthanum borohydrides, Phys. Rev. B 104, 134501 (2021).
  * Lucrezi  _et al._ [2022] R. Lucrezi, S. Di Cataldo, W. von der Linden, L. Boeri, and C. Heil, In-silico synthesis of lowest-pressure high-$T_{\text{c}}$ ternary superhydrides, npj Computational Materials 8, 119 (2022).
  * Sun  _et al._ [2022] Y. Sun, S. Sun, X. Zhong, and H. Liu, Prediction for high superconducting ternary hydrides below megabar pressure, Journal of Physics: Condensed Matter 34, 505404 (2022).
  * Qian  _et al._ [2013] G.-R. Qian, X. Dong, X.-F. Zhou, Y. Tian, A. R. Oganov, and H.-T. Wang, Variable cell nudged elastic band method for studying solid–solid structural phase transitions, Computer Physics Communications 184, 2111 (2013).
  * Belli and Errea [2022] F. Belli and I. Errea, Impact of ionic quantum fluctuations on the thermodynamic stability and superconductivity of LaBH8, Phys. Rev. B 106, 134509 (2022).
  * Shapeev [2016] A. V. Shapeev, Moment tensor potentials: A class of systematically improvable interatomic potentials, Multiscale Modeling & Simulation 14, 1153 (2016).
  * Gubaev  _et al._ [2019] K. Gubaev, E. V. Podryabinkin, G. L. Hart, and A. V. Shapeev, Accelerating high-throughput searches for new alloys with active learning of interatomic potentials, Computational Materials Science 156, 148 (2019).
  * Ranalli  _et al._ [2022] L. Ranalli, C. Verdi, L. Monacelli, M. Calandra, G. Kresse, and C. Franchini, Temperature-dependent anharmonic phonons in quantum paraelectric KTaO3 by first principles and machine-learned force fields (2022), arXiv preprint.
  * Pedrielli  _et al._ [2022] A. Pedrielli, P. E. Trevisanutto, L. Monacelli, G. Garberoglio, N. M. Pugno, and S. Taioli, Understanding anharmonic effects on hydrogen desorption characteristics of MgnH2n nanoclusters by ab initio trained deep neural network, Nanoscale 14, 5589 (2022).
  * Zuo  _et al._ [2020] Y. Zuo, C. Chen, X. Li, Z. Deng, Y. Chen, J. Behler, G. Csányi, A. V. Shapeev, A. P. Thompson, M. A. Wood, and S. P. Ong, Performance and cost assessment of machine learning interatomic potentials, The Journal of Physical Chemistry A 124, 731 (2020).
  * Born [1951] M. Born, Die Gültigkeitsgrenze der Theorie der idealen Kristalle und ihre Überwindung, in  _Festschrift zur Feier des Zweihundertjährigen Bestehens der Akademie der Wissenschaften in Göttingen: I. Mathematisch-Physikalische Klasse_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 1951) pp. 1–16.
  * Hooton [1955] D. J. Hooton, LI. A new treatment of anharmonicity in lattice thermodynamics: I, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 46, 422 (1955).
  * Koehler [1966] T. R. Koehler, Theory of the Self-Consistent Harmonic Approximation with Application to Solid Neon, Phys. Rev. Lett. 17, 89 (1966).
  * Werthamer [1970] N. R. Werthamer, Self-consistent phonon formulation of anharmonic lattice dynamics, Phys. Rev. B 1, 572 (1970).
  * Errea  _et al._ [2014] I. Errea, M. Calandra, and F. Mauri, Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides, Phys. Rev. B 89, 064302 (2014).
  * Note [1] The Supplemental Material is available on the journal website.
  * Errea  _et al._ [2013] I. Errea, M. Calandra, and F. Mauri, First-principles theory of anharmonicity and the inverse isotope effect in superconducting palladium-hydride compounds, Phys. Rev. Lett. 111, 177002 (2013).
  * Bianco  _et al._ [2017] R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Phys. Rev. B 96, 014111 (2017).
  * Ribeiro  _et al._ [2018] G. A. S. Ribeiro, L. Paulatto, R. Bianco, I. Errea, F. Mauri, and M. Calandra, Strong anharmonicity in the phonon spectra of PbTe and SnTe from first principles, Phys. Rev. B 97, 014306 (2018).
  * Monacelli  _et al._ [2018] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Pressure and stress tensor of complex anharmonic crystals within the stochastic self-consistent harmonic approximation, Phys. Rev. B 98, 024106 (2018).
  * Aseginolaza  _et al._ [2019] U. Aseginolaza, R. Bianco, L. Monacelli, L. Paulatto, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Phonon collapse and second-order phase transition in thermoelectric SnSe, Phys. Rev. Lett. 122, 075901 (2019).
  * Errea  _et al._ [2020] I. Errea, F. Belli, L. Monacelli, A. Sanna, T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri, and J. A. Flores-Livas, Quantum crystal structure in the 250-kelvin superconducting lanthanum hydride, Nature 578, 66 (2020).
  * Monacelli  _et al._ [2021] L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, The stochastic self-consistent harmonic approximation: calculating vibrational properties of materials with full quantum and anharmonic effects, Journal of Physics: Condensed Matter 33, 363001 (2021).
  * Note [2] We found MTPs to be a factor of about $20\,000\text{\,}$ faster than DFT in the case of 2$\times$2$\times$2 supercells of BaSiH8.
  * Note [3] Within SSCHA, the stress tensor is obtained as $P_{\alpha\beta}=-V_{\text{sc}}^{-1}\left(\partial\mathcal{F}/\partial\varepsilon_{\alpha\beta}\right)|_{\varepsilon=0}$, where $V_{\text{sc}}$ is the supercell volume, $\mathcal{F}$ the free energy functional, and $\varepsilon_{\alpha\beta}$ the strain tensor [37]. The pressure is then obtained as $\tilde{p}=\sum@\slimits@_{\alpha}P_{\alpha\alpha}/3$.
  * [43] The small imaginary dip near $\Gamma$ in b is due to interpolation, all phonon frequencies at wave vectors commensurate with the supercell are positive.
  * Note [4] Calculating fourth-order corrections on larger supercells is computationally unfeasible with our currently available resources. This calculation in an $n\times n\times n$ supercell requires 4D arrays with $(3\cdot n^{3}\cdot<\text{number of atoms in uc}>)^{4}$ elements. For 10 atoms in the unit cell and $n=3$, this yields $810^{4}$ double-precision entries, for which $3.13\text{\,}\mathrm{T}\mathrm{B}$ of RAM are needed on a single shared-memory node.
  * Born and Oppenheimer [1927] M. Born and R. Oppenheimer, Zur Quantentheorie der Molekeln, Annalen der Physik 389, 457 (1927).
  * Note [5] In the case of fixed ions, we obtain non-vanishing forces within SSCHA. As the individual force components on the H atoms are still small ($\sim$$150\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}\mathrm{/}\mathrm{\text{\AA}}$ in the cell with $a=$6.242\text{\,}\mathrm{\text{\AA}}$$, for example), a fixed-ion calculation seems applicable for a qualitative insight. The forces on atom $i$ are obtained as $\mathbf{F}_{i}=-\partial\mathcal{F}/\partial\mathcal{\mathbf{R}}_{i}$.
  * Poncé  _et al._ [2016] S. Poncé, E. Margine, C. Verdi, and F. Giustino, EPW: Electron–phonon coupling, transport and superconducting properties using maximally localized Wannier functions, Comput. Phys. Commun. 209, 116 (2016).
  * Note [6] A more rigorous treatment would require the replacement of the polarization vectors $\mathbf{e}_{\mathbf{q}\nu}$, as well as either the inclusion of the force term arising in the DF(P)T description of the structure defined by $\tilde{x}$, or an $\mathcal{F}$-based description of the $ep$ matrix elements within the SSCHA.
  * Note [7] The $ep$ coupling strength $\lambda(\omega)=2\intop\ilimits@_{0}^{\omega}\mathrm{d}\omega^{\prime}\frac{\alpha^{2}F(\omega^{\prime})}{\omega^{\prime}}$ and the logarithmic average phonon frequency $\omega_{\text{log}}=\mathop{exp}\nolimits\left(\frac{2}{\lambda}\intop\ilimits@_{0}^{\infty}\mathrm{d}\omega\frac{\alpha^{2}F(\omega)\mathop{ln}\nolimits\omega}{\omega}\right)$ are obtained from the Eliashberg spectral function $\alpha^{2}F(\omega)$.
  * Heil  _et al._ [2019] C. Heil, S. di Cataldo, G. B. Bachelet, and L. Boeri, Superconductivity in sodalite-like yttrium hydride clathrates, Phys. Rev. B 99, 220502 (2019).
  * Kong  _et al._ [2021] P. Kong, V. S. Minkov, M. A. Kuzovnikov, A. P. Drozdov, S. P. Besedin, S. Mozaffari, L. Balicas, F. F. Balakirev, V. B. Prakapenka, S. Chariton, D. A. Knyazev, E. Greenberg, and M. I. Eremets, Superconductivity up to 243 K in the yttrium-hydrogen system under high pressure, Nature Communications 12, 5075 (2021).
  * Errea  _et al._ [2016] I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system, Nature 532, 81 (2016).
  * Einaga  _et al._ [2016] M. Einaga, M. Sakata, T. Ishikawa, K. Shimizu, M. I. Eremets, A. P. Drozdov, I. A. Troyan, N. Hirao, and Y. Ohishi, Crystal structure of the superconducting phase of sulfur hydride, Nature Physics 12, 835 (2016).
  * Giannozzi  _et al._ [2017] P. Giannozzi _et al._ , Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys. Condens. Matter. 29, 465901 (2017).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Note [8] The effective sample size is calculated as $N_{\text{eff}}=(\sum@\slimits@_{i}\rho_{i}^{2})/(\sum@\slimits@_{i}\rho_{i})^{2}$, where the $\rho_{i}$ are the importance sampling weighting factors.
  * Novikov  _et al._ [2020] I. S. Novikov, K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, The MLIP package: moment tensor potentials with MPI and active learning, Machine Learning: Science and Technology 2, 025002 (2020).
  * Margine and Giustino [2013] E. R. Margine and F. Giustino, Anisotropic Migdal-Eliashberg theory using Wannier functions, Phys. Rev. B 87, 024505 (2013).


# Evaluating COVID-19 Sequence Data Using Nearest-Neighbors Based Network
Model

Sarwan Ali 2022 IEEE International Conference on Big Data (Big Data) —
978-1-6654-8045-1/22/$31.00 ©2022 European Union  Department of Computer
Science
Georgia State University, Atlanta, USA
<EMAIL_ADDRESS>
###### Abstract

The SARS-CoV-2 coronavirus is the cause of the COVID-19 disease in humans.
Like many coronaviruses, it can adapt to different hosts and evolve into
different lineages. It is well-known that the major SARS-CoV-2 lineages are
characterized by mutations that happen predominantly in the spike protein.
Understanding the spike protein structure and how it can be perturbed is vital
for understanding and determining if a lineage is of concern. These are
crucial to identifying and controlling current outbreaks and preventing future
pandemics. Machine learning (ML) methods are a viable solution to this effort,
given the volume of available sequencing data, much of which is unaligned or
even unassembled. However, such ML methods require fixed-length numerical
feature vectors in Euclidean space to be applicable. Similarly, euclidean
space is not considered the best choice when working with the classification
and clustering tasks for biological sequences. For this purpose, we design a
method that converts the protein (spike) sequences into the sequence
similarity network (SSN). We can then use SSN as an input for the classical
algorithms from the graph mining domain for the typical tasks such as
classification and clustering to understand the data. We show that the
proposed alignment-free method is able to outperform the current SOTA method
in terms of clustering results. Similarly, we are able to achieve higher
classification accuracy using well-known Node2Vec-based embedding compared to
other baseline embedding approaches.

###### Index Terms:

Sequence Classification, KNN Graph, Spike Sequence, COVID-19, Sequence
Clustering, $k$-mers

##  I Introduction

Because of the COVID-19 pandemic, a lot of sequence data is publicly available
on databases such as GISAID 111https://www.gisaid.org/. This motivated the
researchers to evaluate the data using multiple tools for the purpose of
advanced analysis. Networks are one tool to this end and have found many
applications in bioinformatics. More specifically, sequence similarity
networks are useful for the analysis of biological data using well-established
classification and clustering methods. Clustering the sequences based on the
lineages, locality, and time could help us understand the evolution and spread
of different lineages, which indeed helps the biologists and relevant
government authorities to take appropriate measures in advance.

Some effort has been made recently to classify and cluster sequences based on
hosts [25, 4] and lineages [9, 35, 7]. Clustering approaches can help in
identifying novel and rapidly growing lineages, while classification can
assist in keeping the track of existing ones. Since the spike (S protein)
sequence (see Figure 1) of a coronavirus is the point of contact to the host
cell, it is a vital characteristic of this type of virus [26, 37]. Therefore,
the spike the region is often the focus of the literature (rather than the
entire genome) for studying the behavior of coronaviruses, both in terms of
host and lineage specificity [25, 12]. However, for a spike sequence to be a
compatible input to machine learning (ML) models, it must be transformed into
a fixed-length numerical representation, known as the feature vector. There
are many methods proposed in the literature to produce such a representation
of a spike sequences, such as one-hot encoding [25], $k$-mers based encoding
[7], and position weight matrix based encoding [4].

Figure 1: The full-length genome sequence contains structural and non-
structural parts. We are interested in spike (s) region as mutations related
to coronavirus happen in that region.

Our contributions to this paper are the following:

  1. 1.

A new method to convert the biological sequences into a sequence similarity
graph is proposed, which can be used to study the sequences efficiently.

  2. 2.

Using the clustering algorithms, we show that the proposed model outperforms
the current SOTA approach using different internal clustering evaluation
metrics.

  3. 3.

Using different embedding representation methods for nodes, we show that node
classification using Node2Vec gives higher accuracy than other embedding
approaches.

The remaining paper is structured as follows. Section II contains related
work. Section III contains the proposed approach. Section IV contains
experimental setup and dataset statistics. Section V contains results and
discussion. Section VI concludes this paper.

##  II Related Work

In the biology domain, sequence analysis using the Phylogeny-based method [21]
is an important task. The method that uses substring (called mers) of length
$k$ (hence $k$-mers) counts was first explored in [15] for phylogenetic
applications. Authors in [15] proposed a phylogenetic tree-based method from
the sequences (non-coding and coding regions). In recent years, since we have
had a huge availability of sequencing data due to the coronavirus pandemic,
the sequence analysis problem attracts the attention of researchers [23, 18,
17]. Few methods proposed for supervised tasks are alignment-free. Similarly,
some approaches depend on the alignment of sequences for both supervised and
unsupervised tasks.

The conversion of sequences into embeddings is an important step in the
machine-learning pipeline while working with supervised and unsupervised
tasks. Numerical embedding generation is important in many fields like time
series forecasting [5, 2], analysis of graphs [11, 28], electromyography [36],
clinical data analysis [13], and network security [3]. A method to convert the
aligned sequences into fixed-length embeddings using a position weight matrix
is proposed in [4]. Their method shows improvement in terms of predictive
performance. However, the sequence alignment step in their pipeline makes it
difficult to use for a large number of sequences. A method based on $k$-mers
spectrum for supervised analysis of biological sequences is proposed in [7,
10]. For unsupervised analysis, many recent studies proposed embedding methods
for the biological sequences [12, 35, 8].

Using kernel matrix for supervised sequence analysis is another domain
explored in the literature. A kernel (gram) matrix is generated in this case
that corresponds to pairwise distances between sequences [8, 24]. This gram
matrix is then used with kernel-based classifiers for supervised analysis. A
kernel matrix-based method is proposed in [9] for biological sequence
classification. Although the proposed method shows promising results, it is
not applicable in real-world scenarios due to the expensive storage cost of
the kernel matrix.

##  III Proposed Approach

In this section, we describe the overall approach that we are using in this
paper. Our proposed method is divided into different steps described below.

###  III-A Generating Fixed-Length Numerical Representation

Given a set of spike sequences and attributes (coronavirus lineages) related
to them, we want to generate sequence similarity network (SSN). SSNs is a
network in which nodes are sequences and edges show the top $K$ nearest
neighbors to any given sequence. Since the KNN algorithm works in euclidean
space, the first step is to convert the biological sequences into a fixed-
length numerical representation. For this purpose, we use a method called
$k$-mers [9] (also called n-gram in the NLP domain) as given in Figure 2.

######  Definition III.1 ($k$-mers)

Given some sequence $\sigma$ on alphabet $\Sigma$, we generates substrings
(also called mers) of length $k$, _i.e._ , $k$-mers (using sliding window
approach. In $\Sigma$, we have the following 20 characters (amino acids)
“ACDEFGHIKLMNPQRSTVWY”.

For any sequence of length $N$, the total $k$-mers are:

$\text{total $k$-mers = }(N-k)+1$ (1)

Figure 2: $k$-mers example for $k=3$.

We use $k=3$ for the $k$-mers, computed using the validation set. After
generating the $k$-mers for a given sequence, we create the frequency vector
(numerical representation) of length $\Sigma^{k}$ (i.e., length comprised of
all possible $k$-mers in $\Sigma$ of length $k$), which contains the
count/frequencies of $k$-mers. The pseudocode to compute the frequency vectors
is given in Algorithm 1.

Algorithm 1 ComputeFrequencyVector

1:Input: Set $\mathcal{S}$ of $k$-mers on alphabet $\Sigma$, the size of mers
$n$

2:Output: Frequency Vector $V$

3:combos = GenerateAllKmersCombinations($\Sigma$, n)

4:$V$ = [0] * $|\Sigma|^{n}$ $\triangleright$ Total length of (zero) vector

5:for i $\leftarrow 1$ to $|\mathcal{S}|$ do

6: idx = combos.index($\mathcal{S}$[i]) $\triangleright$ Find index of
$i^{th}$ k-mer

7: $V$[idx] $\leftarrow$ $V$[idx] + 1 $\triangleright$ Increment bin by 1

8:return($V$)

###  III-B Generating Sequence Similarity Network (SSN)

After generating the $k$-mers-based numerical representation, we use the
$K$-Nearest Neighbors-based approach to generate the Sequence Similarity
Network (SSN), where nodes are the sequences and edges between them show the
similarity.

###### Remark 1

Note that we generate an unweighted SSN.

We use $K=20$ in this case (decided using standard validation set approach
[20]). The resultant SSN graph is given in Figure 3.

Figure 3: Input SSN Graph.

After generating the SSN, we use several supervised (classification) and
unsupervised (clustering) methods to perform sequence analysis.

###  III-C Supervised Analysis

Given the SSN, we perform lineage classification using different machine
learning (ML) algorithms in this setting. To apply the ML algorithms on
graphs, we need to get the feature embeddings for the nodes, which can be used
as input to the ML classifiers. To get the numerical embeddings, we use the
following methods:

####  III-C1 Locally Linear Embedding (LLE) [32]

Given a node $v$ from the graph along with its neighbors $N(v)$, this method
assumes that $v$ is a linear combination of $N(v)$ in the embedding space
(where this linear combination property must hold for all nodes in the graph).
It uses the objective function to minimize the distances between the
embeddings of two nodes that are neighbors to each other. The embeddings are
normalized in the end to get the final representation.

####  III-C2 Laplacian Eigenmaps [14]

Using the graph Laplacian, it generates the embedding for each node while
preserving the local neighborhood information for that node. Similar to LLE,
it also solves the objective function to keep the embeddings close to each
other for the nodes that are neighbors. However, it involves graph Laplacian.

####  III-C3 Higher-Order Proximity preserved Embedding (HOPE) [29]

Given the similarity matrix, it preserves the higher order proximity between
the embedding of nodes by minimizing the objective function in which we
minimize the squared differences between the similarity matrix and the pairs
of given embeddings.

####  III-C4 Graph Factorization (GF) [1]

Given the adjacency matrix of the graph, this method factorizes that matrix.
The GF minimizes the objective function (reducing the distance between the
neighboring nodes in the corresponding embeddings) to give an approximate
embedding representation for a given node. Since it does not give an exact
solution, the embeddings may contain some noise.

####  III-C5 DeepWalk [30]

Given a pair of nodes, the goal of this method is to preserve higher-order
proximity. This goal is fulfilled by maximizing the probability of observing
in the random walk (previous and next top $k$ nodes).

####  III-C6 Node2Vec [22]

It is similar to DeepWalk. The difference between the two is that this
approach adopts biased-random walks. This bias behavior gives a trade-off
between depth-first graph search and breadth-first graph search. Therefore it
results in efficient embeddings compared to DeepWalk.

###  III-D Unsupervised Analysis

In this setting, we perform clustering using different unsupervised clustering
algorithms. We use the following clustering algorithms.

####  III-D1 MiniBatch KMeans

It is based on the simple idea of KMeans with mini-batch optimization. The
idea of using a mini-batch is to reduce the runtime cost of the algorithm.

####  III-D2 Affinity Propagation

Given measures of similarity between data points, it uses the concept of
message passing to cluster the data. It does not require the number of
clusters as input. However, it has a quadratic runtime.

####  III-D3 Mean Shift

It is a centroid-based method that updates the candidates for centroids so
that they can be the mean of the points within a given region.

####  III-D4 Spectral Clustering

This method uses the spectrum (eigenvalues) of the kernel matrix (also called
gram/similarity matrix), which is computed from the input data to perform
clustering (after performing dimensionality reduction) in fewer dimensions.

####  III-D5 Ward

It uses agglomerative hierarchical clustering along with an objective function
so that merging two clusters at any iterations should satisfy the objective
function (minimizing error sum of squares ESS).

####  III-D6 Agglomerative Clustering

This method combines the pairs of clusters recursively from the data points.
It uses the dendrogram-based approach.

####  III-D7 Density-Based Spatial Clustering of Applications with Noise
(DBSCAN)

Given the input data, DBSCAN groups together data points that are close to
each other (nearest neighbors). It also keeps track of the outliers data
points, which lies alone in low-density regions (points that do not have
nearby neighbors).

####  III-D8 Ordering Points To Identify the Clustering Structure (OPTICS)

This method searches for the core sample of high density and expands clusters
from them. It is similar to DBSCAN. One of the main differences OPTICS has as
compared to the DBSCAN is that it keeps cluster hierarchy for a variable
neighborhood radius.

####  III-D9 BIRCH

It works by constructing a tree data structure. In that tree, the cluster
centroids behave as the leaf. Using the threshold value, clusters could be
extracted from the tree. Similarly, other algorithms can be used in
combination to get the final clusters, such as agglomerative clustering.

####  III-D10 Gaussian Mixture

It contains a mixture of multiple gaussian distributions. It is similar to the
KMeans algorithm. However, it is considered a Soft clustering approach
(overlapping allowed), while KMeans is considered a hard clustering algorithm
(no overlapping). Given the data, this approach assigns the data points to the
multivariate normal components. This assignment should maximize the component
posterior probability.

###### Remark 2

Note that all clustering methods except ward and agglomerative clustering take
$k$-mers-based numerical embeddings as input. The ward and agglomerative
clustering approaches take SSN as input.

##  IV Experimental Setup

We split the data into $70$-$30\%$ for training and testing (held out) sets,
respectively, for experimentation. We use $5$ fold cross-validation on the
training data for hyperparameter tuning. The final results are computed for
$30\%$ held-out testing set. We use $5$ random initialization of data to avoid
biases in the results. All experiments are performed on Windows 10 operating
system, having a Core i5 processor and 32 GB RAM. The code is written in
Python, which is available online
222https://github.com/sarwanpasha/Spike2Network.

###  IV-A Dataset Statistics

We extracted spike sequence data from GISAID 333https://www.gisaid.org/
website (as given in [6]), which consists of $7000$ sequences of length 1274
from $22$ lineages. The proportion of lineages in the dataset is given in
Table I.

TABLE I: Dataset statistics for $22$ lineages [6]. The character ‘-’ means
that information is not available.

Lineage |  Region of First Time Detection |  Variant Name |  No. Mut. S/Gen. |  No. of sequences
---|---|---|---|---
B.1.1.7 |  UK |  Alpha |  8/17 |  3369
B.1.617.2 |  India |  Delta |  8/17 |  875
AY.4 |  India |  Delta |  - |  593
B.1.2 |  USA |  - |  - |  333
B.1 |  USA |  |  |  292
B.1.177 |  Spain |  - |  - |  243
P.1 |  Brazil |  Gamma |  10/21 |  194
B.1.1 |  UK |  |  - |  163
B.1.429 |  California |  Epsilon |  3/5 |  107
B.1.526 |  New York |  Iota |  6/16 |  104
AY.12 |  India |  Delta |  - |  101
B.1.160 |  France |  - |  - |  92
B.1.351 |  South Africa |  Beta |  9/21 |  81
B.1.427 |  California |  Epsilon |  3/5 |  65
B.1.1.214 |  Japan |  - |  - |  64
B.1.1.519 |  USA |  - |  - |  56
D.2 |  Australia |  - |  - |  55
B.1.221 |  Netherlands |  - |  - |  52
B.1.177.21 |  Denmark |  - |  - |  47
B.1.258 |  Germany |  - |  - |  46
B.1.243 |  USA |  - |  - |  36
R.1 |  Japan |  - |  - |  32
Total |  - |  - |  - |  7000

###  IV-B Elbow Method for Optimal Number of Clusters

We use a data-driven method called Elbow, to compute the optimal number of
clusters [33, 12]. In this method, we evaluate the trade-off between runtime
and the sum of squared error and get the number of clusters where these two
metrics are smaller. In this paper, we selected $4$ as the optimal cluster
number (see Figure 4).

Figure 4: Elbow method using different value of $k$ to get optimal value. The
notation $k$ represents the number of clusters.

###  IV-C Baseline Model

As an unsupervised baseline, we use a method, called PWM2Vec, which is
proposed in [4]. This approach works by computing weights for each $k$-mers
using the position weight matrix[34] and then concat all those weights for all
$k$-mers within a sequence to get the final embedding.

###  IV-D Evaluation Metrics for Clustering

For the analysis of the unsupervised task, we use $k$-means clustering. The
following evaluation metrics are used to measure the quality of clusters:

Silhouette Coefficient [31]: This metric measures the similarity of a vector
with its own cluster (cohesion) and with vectors in other clusters
(separation). The value of this metric is between $-1$ (worst) to $1$ (best).

Calinski-Harabasz Score [16]: This metric works by analyzing inter-cluster
dispersion and the between-cluster dispersion (computing ratio between the
two). A lower score of this metric means bad clustering and vice versa.

Davies-Bouldin Score [19]: This metric computes the similarity between
clusters by measuring the ratio of distances within-cluster to between
clusters. In this case, a lower value corresponds to better clustering and
vice versa.

###  IV-E Evaluation Metrics For Classification

We use multiple classifiers to test the performance of embeddings. The
resultant feature vectors of length $200$ (using parameters tunning) are fed
to different classifiers as input. We use Multi-Layer Perceptron (MLP), SVM,
Decision Tree (DT), Naive Bayes (NB), Random Forest (RF), K-Nearest Neighbour
(KNN), and Logistic Regression (LR) algorithms. The evaluation metrics used to
test the performance are recall, precision, accuracy, weighted F1, macro F1,
and ROC area under the curve (AUC). To use the metrics designed for the binary
classification task, we use the one-vs-rest approach to use them for multi-
class classification.

##  V Results and Discussion

In this section, we show the results of classification and clustering methods.

###  V-A Clustering Results

####  V-A1 Subjective Evaluation

We use t-distributed stochastic neighbor embedding (t-SNE) [27] to represent
the data in 2-dimensions so that we can use scatter plots to visualize the
data and analyze the patterns. The t-SNE method maps input data to 2D real
vectors (low dimensions) while preserving the pair-wise distance between the
points from high dimensions. The subjective evaluation of different clustering
methods (using t-SNE) is shown in Figure 5. The t-SNE-based plot for the data
with original labels is also shown in Figure 5(k) for comparison purposes.
Compared with the original plot, we can conclude that the gaussian mixture and
ward give the clustering that is closely related to the original plot in some
sense.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 5: Subjective evaluation from different clustering approaches using
t-SNE.

####  V-A2 Objective Evaluation

The results for objective evaluation are given in Table II. In general, we can
observe that there is no clear winner for the goodness of clustering. For the
Silhouette coefficient, the Ward method performs the best. In the case of the
Calinski-Harabasz score, the Gaussian Mixture model shows the best
performance. Similarly, Agglomerative clustering shows the best performance in
the case of Davies Bouldin score. PWM2Vec (the SOTA method) performs best in
terms of clustering runtime because of lower dimensional embeddings.

TABLE II: Objective evaluation using internal clustering quality metrics for
different clustering algorithms. The best values are shown in bold.

| Evaluation Metrics |
---|---|---
Algorithm |  Silhouette Coefficient |  Calinski Harabasz Score |  Davies-Bouldin Score |  Clustering Runtime (Sec.)
PWM2Vec |  0.477 |  1762.983 |  1.007 |  1.45
MiniBatch KMeans |  0.133 |  1505.001 |  2.064 |  53.43
Affinity Propagation |  -0.600 |  0.111 |  1.901 |  511.47
Mean Shift |  -0.768 |  0.874 |  5.175 |  1326.34
Spectral Clustering |  -0.613 |  402.117 |  1.984 |  219.44
Ward |  0.593 |  2657.689 |  5.152 |  322.17
Agglomerative Clustering |  0.290 |  95.723 |  0.455 |  155.35
DBSCAN |  0.216 |  23.174 |  1.351 |  14.55
OPTICS |  -0.196 |  285.435 |  1.548 |  4229.51
BIRCH |  0.524 |  2242.030 |  4.757 |  227.01
Gaussian Mixture |  -0.418 |  2764.134 |  1.463 |  1030.63

###  V-B Classification Results

The classification results for different embedding methods and classification
algorithms are given in Table III. We can observe that Node2Vec significantly
outperforms all other embeddings for all evaluation metrics using SVM and KNN
classifiers. This behavior shows that the embeddings generated using Node2Vec
were able to preserve the structure of the network more efficiently than the
other embedding approaches.

TABLE III: Classification results (averaged over $5$ runs) using different
evaluation metrics. The best values are shown in bold.

Embed. |  |  Acc. |  Prec. |  Recall |  F1 (Weig.) |  F1 (Macro) |  ROC AUC |  Train Time (Sec.)
---|---|---|---|---|---|---|---|---
Laplacian Eigenmaps | SVM |  0.47 |  0.22 |  0.47 |  0.30 |  0.03 |  0.50 |  3.45
NB |  0.01 |  0.11 |  0.01 |  0.01 |  0.01 |  0.50 |  0.11
MLP |  0.40 |  0.25 |  0.40 |  0.30 |  0.04 |  0.50 |  12.98
KNN |  0.37 |  0.25 |  0.37 |  0.29 |  0.04 |  0.50 |  0.33
RF |  0.39 |  0.25 |  0.39 |  0.29 |  0.04 |  0.50 |  7.84
LR |  0.47 |  0.22 |  0.47 |  0.30 |  0.03 |  0.50 |  2.18
DT |  0.30 |  0.27 |  0.30 |  0.28 |  0.05 |  0.50 |  1.70
Locally Linear Embedding | SVM |  0.49 |  0.24 |  0.49 |  0.32 |  0.03 |  0.50 |  3.28
NB |  0.01 |  0.26 |  0.01 |  0.01 |  0.01 |  0.50 |  0.12
MLP |  0.42 |  0.28 |  0.42 |  0.32 |  0.04 |  0.50 |  20.98
KNN |  0.39 |  0.27 |  0.39 |  0.31 |  0.04 |  0.50 |  0.34
RF |  0.36 |  0.26 |  0.36 |  0.30 |  0.04 |  0.50 |  7.64
LR |  0.49 |  0.24 |  0.49 |  0.32 |  0.03 |  0.50 |  1.16
DT |  0.26 |  0.27 |  0.26 |  0.27 |  0.05 |  0.50 |  1.99
HOPE Embedding | SVM |  0.48 |  0.23 |  0.48 |  0.31 |  0.03 |  0.50 |  2.57
NB |  0.01 |  0.24 |  0.01 |  0.01 |  0.01 |  0.50 |  0.12
MLP |  0.45 |  0.27 |  0.45 |  0.32 |  0.04 |  0.50 |  18.77
KNN |  0.38 |  0.26 |  0.38 |  0.30 |  0.04 |  0.50 |  0.30
RF |  0.37 |  0.25 |  0.37 |  0.29 |  0.03 |  0.50 |  6.68
LR |  0.48 |  0.23 |  0.48 |  0.31 |  0.03 |  0.50 |  1.01
DT |  0.30 |  0.26 |  0.30 |  0.28 |  0.04 |  0.50 |  1.64
Graph Factorization | SVM |  0.46 |  0.21 |  0.46 |  0.29 |  0.03 |  0.50 |  2.92
NB |  0.37 |  0.51 |  0.37 |  0.40 |  0.31 |  0.64 |  0.11
MLP |  0.54 |  0.51 |  0.54 |  0.53 |  0.22 |  0.60 |  18.46
KNN |  0.69 |  0.69 |  0.69 |  0.69 |  0.51 |  0.73 |  0.39
RF |  0.54 |  0.53 |  0.54 |  0.44 |  0.14 |  0.54 |  5.47
LR |  0.46 |  0.21 |  0.46 |  0.29 |  0.03 |  0.50 |  1.26
DT |  0.40 |  0.40 |  0.40 |  0.40 |  0.14 |  0.56 |  1.28
DeepWalk | SVM |  0.48 |  0.23 |  0.48 |  0.31 |  0.03 |  0.50 |  2.73
NB |  0.46 |  0.26 |  0.46 |  0.32 |  0.03 |  0.50 |  0.13
MLP |  0.32 |  0.26 |  0.32 |  0.29 |  0.04 |  0.50 |  23.40
KNN |  0.38 |  0.26 |  0.38 |  0.30 |  0.04 |  0.50 |  0.33
RF |  0.48 |  0.23 |  0.48 |  0.31 |  0.03 |  0.50 |  9.27
LR |  0.48 |  0.23 |  0.48 |  0.31 |  0.03 |  0.50 |  1.15
DT |  0.25 |  0.25 |  0.25 |  0.25 |  0.04 |  0.50 |  2.37
Node2Vec | SVM |  0.79 |  0.74 |  0.79 |  0.77 |  0.52 |  0.75 |  1.40
NB |  0.63 |  0.70 |  0.63 |  0.65 |  0.48 |  0.75 |  0.10
MLP |  0.74 |  0.74 |  0.74 |  0.74 |  0.47 |  0.73 |  19.83
KNN |  0.76 |  0.76 |  0.76 |  0.76 |  0.55 |  0.77 |  0.32
RF |  0.78 |  0.73 |  0.78 |  0.75 |  0.51 |  0.74 |  4.97
LR |  0.78 |  0.72 |  0.78 |  0.75 |  0.49 |  0.73 |  3.81
DT |  0.68 |  0.68 |  0.68 |  0.68 |  0.42 |  0.71 |  1.06

To evaluate whether the computed results are statistically significant, we
used the student t-test and evaluated the $p$-values for the results using the
average and standard deviation results of $5$ runs. The standard deviation
results are reported in Table IV. We noted that the $p$-values were $<0.05$ in
the majority of the cases and for all embedding methods (because standard
deviation values are very low), hence confirming the statistical significance
of the reported results. Note that we have not reported the $p$-values in this
paper due to space constraints. However, we believe that they can easily be
computed by anyone using the reported average and standard deviation results.

TABLE IV: Standard Deviation values of 5 runs for Classification results on
the proposed and SOTA methods. The average results are reported in Table III.

Embed. Method |  ML Algo. |  Acc. |  Prec. |  Recall |  F1 weigh. | F1 Macro |  ROC- AUC |  Train. runtime (sec.)
---|---|---|---|---|---|---|---|---
Laplacian Eigenmaps |  SVM |  0.005339 |  0.005547 |  0.005339 |  0.004908 | 0.00894 |  0.003092 |  3.916014
NB |  0.165977 |  0.047923 |  0.165977 |  0.124314 | 0.028412 |  0.021885 |  0.234008
MLP |  0.010331 |  0.011387 |  0.010331 |  0.010000 | 0.024026 |  0.013581 |  5.913928
KNN |  0.015266 |  0.00962 |  0.015266 |  0.014351 | 0.016323 |  0.006342 |  2.590869
RF |  0.007615 |  0.010725 |  0.007615 |  0.008411 | 0.016425 |  0.007389 |  0.433674
LR |  0.006048 |  0.005895 |  0.006048 |  0.006385 | 0.01425 |  0.005388 |  2.752432
DT |  0.004937 |  0.005132 |  0.004937 |  0.004463 | 0.003758 |  0.002974 |  0.437479
Locally Linear Embedding |  SVM |  0.009732 |  0.014177 |  0.009732 |  0.011394 | 0.040983 |  0.022625 |  0.111355
NB |  0.04711 |  0.028739 |  0.04711 |  0.034415 | 0.021785 |  0.019178 |  0.047905
MLP |  0.004077 |  0.010426 |  0.004077 |  0.005091 | 0.020858 |  0.007274 |  2.914309
KNN |  0.011276 |  0.010159 |  0.011276 |  0.011393 | 0.017793 |  0.007173 |  0.048912
RF |  0.004304 |  0.01031 |  0.004304 |  0.007477 | 0.00978 |  0.007079 |  0.784592
LR |  0.003573 |  0.00751 |  0.003573 |  0.004969 | 0.010045 |  0.005569 |  5.316077
DT |  0.005904 |  0.00566 |  0.005904 |  0.004654 | 0.012841 |  0.005028 |  0.058772
HOPE Embedding |  SVM |  0.010442 |  0.013078 |  0.010442 |  0.011673 | 0.022631 |  0.007962 |  0.24096
NB |  0.008815 |  0.005557 |  0.008815 |  0.00703 | 0.011761 |  0.009719 |  0.006103
MLP |  0.007526 |  0.010713 |  0.007526 |  0.007969 | 0.020091 |  0.00948 |  2.020901
KNN |  0.006269 |  0.006129 |  0.006269 |  0.005595 | 0.014096 |  0.00755 |  0.050194
RF |  0.008558 |  0.010725 |  0.008558 |  0.011033 | 0.01859 |  0.006795 |  0.194081
LR |  0.006258 |  0.013144 |  0.006258 |  0.007598 | 0.009922 |  0.00391 |  0.285232
DT |  0.00958 |  0.011658 |  0.00958 |  0.01026 | 0.019321 |  0.007745 |  0.034749
Graph Factorization |  SVM |  0.006843 |  0.007175 |  0.006843 |  0.007533 | 0.018218 |  0.006985 |  0.127963
NB |  0.012983 |  0.014996 |  0.012983 |  0.013997 | 0.009849 |  0.005343 |  0.002441
MLP |  0.009597 |  0.016199 |  0.009597 |  0.008501 | 0.009785 |  0.003415 |  1.973822
KNN |  0.012222 |  0.006539 |  0.012222 |  0.011150 | 0.012703 |  0.007168 |  0.027236
RF |  0.004658 |  0.002715 |  0.004658 |  0.004319 | 0.014516 |  0.011376 |  0.128241
LR |  0.005999 |  0.006079 |  0.005999 |  0.005070 | 0.010821 |  0.008120 |  0.062362
DT |  0.002869 |  0.003766 |  0.002869 |  0.001588 | 0.017623 |  0.012498 |  0.001499
DeepWalk |  SVM |  0.014645 |  0.012340 |  0.014645 |  0.015661 | 0.024540 |  0.014596 |  1.106693
NB |  0.041749 |  0.023858 |  0.041749 |  0.033381 | 0.025358 |  0.008237 |  0.036856
MLP |  0.006472 |  0.008860 |  0.006472 |  0.003471 | 0.022909 |  0.010423 |  3.305528
KNN |  0.012241 |  0.018108 |  0.012241 |  0.013407 | 0.037266 |  0.019443 |  0.089361
RF |  0.015366 |  0.013036 |  0.015366 |  0.015552 | 0.021972 |  0.013453 |  0.217082
LR |  0.010999 |  0.012312 |  0.010999 |  0.012642 | 0.017529 |  0.008787 |  20.67326
DT |  0.012319 |  0.015275 |  0.012319 |  0.013563 | 0.020260 |  0.014356 |  0.346420
Node2Vec |  SVM |  0.007835 |  0.007516 |  0.007835 |  0.007297 | 0.022666 |  0.011886 |  0.215113
NB |  0.008105 |  0.023467 |  0.008105 |  0.012723 | 0.024591 |  0.010523 |  0.039378
MLP |  0.004398 |  0.005859 |  0.004398 |  0.005486 | 0.016899 |  0.006182 |  1.631292
KNN |  0.013624 |  0.010276 |  0.013624 |  0.014326 | 0.01157 |  0.006457 |  0.015902
RF |  0.004179 |  0.003099 |  0.004179 |  0.004322 | 0.010787 |  0.007188 |  0.118258
LR |  0.009168 |  0.011003 |  0.009168 |  0.01032 | 0.000823 |  0.000216 |  0.160578
DT |  0.008422 |  0.005511 |  0.008422 |  0.008296 | 0.017698 |  0.008793 |  0.041413

##  VI Conclusion

In this study, we propose a method to convert the protein (spike) sequences
into a graph to use well-established graph mining algorithms to analyze the
biological sequences in both a supervised and unsupervised manner. We show
that clustering the nodes gives better qualities of the clusters as compared
to the feature engineering-based approach. For the supervised analysis, we
show that Node2Vec could be used to generate the feature embeddings for the
nodes, which can give better classification accuracies compared to the other
embedding methods. In the future, we will work towards testing the proposed
model on more number of sequences to test the scalability. Using a neural
network-based model such as GNN for the classification is another interesting
future work. We will also consider other attributes, such as locality
information and calendar values for computing richer embeddings for sequences.

## References

  * [1] Ahmed, A., Shervashidze, N., Narayanamurthy, S., Josifovski, V., Smola, A.J.: Distributed large-scale natural graph factorization. In: International conference on World Wide Web. pp. 37–48 (2013)
  * [2] Ali, S., Mansoor, H., Khan, I., Arshad, N., Khan, M.A., Faizullah, S.: Short-term load forecasting using ami data. preprint, arXiv:1912.12479 (2019)
  * [3] Ali, S., Alvi, M.K., Faizullah, S., Khan, M.A., Alshanqiti, A., Khan, I.: Detecting ddos attack on sdn due to vulnerabilities in openflow. In: 2019 International Conference on Advances in the Emerging Computing Technologies (AECT). pp. 1–6 (2020)
  * [4] Ali, S., Bello, B., Chourasia, P., Punathil, R.T., Zhou, Y., Patterson, M.: Pwm2vec: An efficient embedding approach for viral host specification from coronavirus spike sequences. Biology 11(3), 418 (2022)
  * [5] Ali, S., Mansoor, H., Arshad, N., Khan, I.: Short term load forecasting using smart meter data. In: International Conference on Future Energy Systems. pp. 419–421 (2019)
  * [6] Ali, S., Murad, T., Chourasia, P., Patterson, M.: Spike2signal: Classifying coronavirus spike sequences with deep learning. In: International Conference on Big Data Computing Service and Applications (BigDataService). pp. 81–88 (2022)
  * [7] Ali, S., Patterson, M.: Spike2vec: An efficient and scalable embedding approach for covid-19 spike sequences. In: International Conference on Big Data (Big Data). pp. 1533–1540 (2021)
  * [8] Ali, S., Sahoo, B., Khan, M.A., Zelikovsky, A., Khan, I.U., Patterson, M.: Efficient approximate kernel based spike sequence classification. IEEE/ACM Transactions on Computational Biology and Bioinformatics (2022)
  * [9] Ali, S., Sahoo, B., Ullah, N., Zelikovskiy, A., Patterson, M., Khan, I.: A k-mer based approach for sars-cov-2 variant identification. In: International Symposium on Bioinformatics Research and Applications. pp. 153–164 (2021)
  * [10] Ali, S., Sahoo, B., Zelikovskiy, A., Chen, P.Y., Patterson, M.: Benchmarking machine learning robustness in covid-19 genome sequence classification. arXiv preprint arXiv:2207.08898 (2022)
  * [11] Ali, S., Shakeel, M.H., Khan, I., Faizullah, S., Khan, M.A.: Predicting attributes of nodes using network structure. ACM Transactions on Intelligent Systems and Technology 12(2), 1–23 (2021)
  * [12] Ali, S., Tamkanat-E-Ali, Khan, M.A., Khan, I., Patterson, M.: Effective and scalable clustering of sars-cov-2 sequences. In: International Conference on Big Data Research (ICBDR). pp. 1–8 (2021)
  * [13] Ali, S., Zhou, Y., Patterson, M.: Efficient analysis of covid-19 clinical data using machine learning models. Medical & Biological Engineering & Computing pp. 1–16 (2022)
  * [14] Belkin, M., Niyogi, P.: Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation 15(6), 1373–1396 (2003)
  * [15] Blaisdell, B.: A measure of the similarity of sets of sequences not requiring sequence alignment. Proceedings of the National Academy of Sciences 83, 5155–5159 (1986)
  * [16] Caliński, T., Harabasz, J.: A dendrite method for cluster analysis. Communications in Statistics-theory and Methods 3(1), 1–27 (1974)
  * [17] Chourasia, P., Ali, S., Ciccolella, S., Della Vedova, G., Patterson, M.: Clustering sars-cov-2 variants from raw high-throughput sequencing reads data. In: International Conference on Computational Advances in Bio and Medical Sciences. pp. 133–148 (2022)
  * [18] Chourasia, P., Ali, S., Ciccolella, S., Della Vedova, G., Patterson, M.: Reads2vec: Efficient embedding of raw high-throughput sequencing reads data. arXiv preprint arXiv:2211.08267 (2022)
  * [19] Davies, D.L., Bouldin, D.W.: A cluster separation measure. IEEE transactions on pattern analysis and machine intelligence (2), 224–227 (1979)
  * [20] Devijver, P., Kittler, J.: Pattern recognition: A statistical approach. In: London, GB: Prentice-Hall. pp. 1–448 (1982)
  * [21] Dhar, S., et al.: Tnet: Phylogeny-based inference of disease transmission networks using within-host strain diversity. In: International Symposium on Bioinformatics Research and Applications (ISBRA). pp. 203–216 (2020)
  * [22] Grover, A., Leskovec, J.: node2vec: Scalable feature learning for networks. In: International Conference on Knowledge Discovery & Data Mining (KDD). pp. 855–864 (2016)
  * [23] Krishnan, G., Kamath, S., Sugumaran, V.: Predicting vaccine hesitancy and vaccine sentiment using topic modeling and evolutionary optimization. In: International Conference on Applications of Natural Language to Information Systems (NLDB). pp. 255–263 (2021)
  * [24] Kuksa, P., Khan, I., Pavlovic, V.: Generalized similarity kernels for efficient sequence classification. In: SIAM International Conference on Data Mining (SDM). pp. 873–882 (2012)
  * [25] Kuzmin, K., Adeniyi, A.E., DaSouza Jr, A.K., Lim, D., Nguyen, H., Molina, N.R., Xiong, L., Weber, I.T., Harrison, R.W.: Machine learning methods accurately predict host specificity of coronaviruses based on spike sequences alone. Biochemical and Biophysical Research Communications 533(3), 553–558 (2020)
  * [26] Li, F.: Structure, function, and evolution of coronavirus spike proteins. Annual review of virology 3, 237–261 (2016)
  * [27] Van der Maaten, L., Hinton, G.: Visualizing data using t-sne. Journal of machine learning research 9(11) (2008)
  * [28] Mansoor, H., Ali, S., Alam, S., Khan, M.A., Khan, I., et al.: Impact of missing data imputation on the fairness and accuracy of graph node classifiers. arXiv preprint arXiv:2211.00783 (2022)
  * [29] Ou, M., Cui, P., Pei, J., Zhang, Z., Zhu, W.: Asymmetric transitivity preserving graph embedding. In: International conference on Knowledge discovery and data mining. pp. 1105–1114 (2016)
  * [30] Perozzi, B., Al-Rfou, R., Skiena, S.: Deepwalk: Online learning of social representations. In: International Conference on Knowledge Discovery & Data Mining (KDD). pp. 701–710 (2014)
  * [31] Rousseeuw, P.J.: Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. Journal of computational and applied mathematics 20, 53–65 (1987)
  * [32] Roweis, S., Saul, L.: Nonlinear dimensionality reduction by locally linear embedding. science 290(5500), 2323–2326 (2000)
  * [33] Satopaa, V., et al.: Finding a kneedle in a haystack: Detecting knee points in system behavior. In: International conference on distributed computing systems workshops. pp. 166–171 (2011)
  * [34] Stormo, G.D., Schneider, T.D., Gold, L., Ehrenfeucht, A.: Use of the ‘Perceptron’ algorithm to distinguish translational initiation sites in E. coli. Nucleic Acids Research 10(9), 2997–3011 (1982). https://doi.org/10.1093/nar/10.9.2997, https://doi.org/10.1093/nar/10.9.2997
  * [35] Tayebi, Z., Ali, S., Patterson, M.: Robust representation and efficient feature selection allows for effective clustering of sars-cov-2 variants. Algorithms 14(12), 348 (2021)
  * [36] Ullah, A., Ali, S., Khan, I., Khan, M.A., Faizullah, S.: Effect of analysis window and feature selection on classification of hand movements using emg signal. In: SAI Intelligent Systems Conference (IntelliSys). pp. 400–415 (2020)
  * [37] Walls, A.C., Park, Y.J., Tortorici, M.A., Wall, A., McGuire, A.T., Veesler, D.: Structure, function, and antigenicity of the sars-cov-2 spike glycoprotein. Cell 181(2), 281–292 (2020)