# Equation of state of Bose gases beyond the universal regime

Martí Planasdemunt  Departament de Física, Universitat Politècnica de
Catalunya, Campus Nord B4-B5, 08034 Barcelona, Spain  Jordi Pera
<EMAIL_ADDRESS>Departament de Física, Universitat Politècnica de
Catalunya, Campus Nord B4-B5, 08034 Barcelona, Spain  Jordi Boronat
<EMAIL_ADDRESS>Departament de Física, Universitat Politècnica de
Catalunya, Campus Nord B4-B5, 08034 Barcelona, Spain

(August 28, 2024)

###### Abstract

The equation of state of dilute Bose gases, in which the energy only depends
on the $s$-wave scattering length, is rather unknown beyond the universal
limit. We have carried out a bunch of diffusion Monte Carlo calculations up to
gas parameters of $10^{-2}$ to explore how the departure from the universality
emerges. Using different model potentials, we calculate the energies of the
gas in an exact way, within some statistical noise, and report the results as
a function of the three relevant scattering parameters: the $s$-wave
scattering length $a_{0}$, the $s$-wave effective range $r_{0}$, and the
$p$-wave scattering length $a_{1}$. If the effective range is not large we
observe universality in terms of $a_{0}$ and $r_{0}$ up to gas parameters of
$10^{-2}$. If $r_{0}$ grows the regime of universality in these two parameters
is reduced and effects of $a_{1}$ start to be observed. In the $(a_{0},r_{0})$
universal regime we propose an analytical law that reproduces fairly well the
exact energies.

The experimental realization of the Bose-Einstein condensed state (BEC)
Anderson  _et al._ (1995); Davis  _et al._ (1995) has allowed for the study of
a dilute Bose gas whose main properties were predicted theoretically many
years ago. The high tunability of the interaction strength in these systems
offers a unique platform to gain accurate insight on their properties.
However, this intriguing feature cannot be carried out as far as theory would
like due to the metastability of the dilute BEC in the laboratories Dalfovo
_et al._ (1999) . Three-body losses frequently hinder the achievement of large
densities $n$ and universality in terms of the gas parameter $x=na_{0}^{3}$,
with $a_{0}$ the $s$-wave scattering length, is what is normally observed.
Nevertheless, we already know of experimental results that go beyond
universality: the Fermi polaron in two dimensions Darkwah Oppong  _et al._
(2019); Bombín  _et al._ (2021), dilute Bose-Bose liquid drops Cabrera  _et
al._ (2018); Cikojevic  _et al._ (2020), and dipolar gases Böttcher  _et al._
(2019); Bombín  _et al._ (2024).

The first terms of the equation of state of Bose gases were calculated long
time ago and proved to be dependent only on the gas parameter,

$\frac{E_{\text{u}}}{N}=c_{\text{MF}}\,x+c_{\text{LHY}}\,x^{3/2}+c_{\text{WU}}\,x^{2}\ln
x\ ,$ (1)

with the energy in unities of $\hbar^{2}/(2ma_{0}^{2})$. The universal
constants in Eq. (1) are

$\displaystyle c_{\text{MF}}$ $\displaystyle=$ $\displaystyle 4\pi\ ,$ (2)
$\displaystyle c_{\text{LHY}}$ $\displaystyle=$
$\displaystyle\frac{512\sqrt{\pi}}{15}\ ,$ (3) $\displaystyle c_{\text{WU}}$
$\displaystyle=$ $\displaystyle\frac{32\pi(4\pi-3\sqrt{3})}{3}\ .$ (4)

The first term in Eq. (1) is the mean-field or Hartree-Fock term Bogoliubov
(1947). The second one was first calculated by Lee, Huang, and Yang Lee  _et
al._ (1957) and incorporates quantum fluctuations at the lowest order, and
finally the last one was obtained by Wu Wu (1959). The LHY and WU terms were
first derived by a hard-sphere interaction and, later on, it was proved that
they are also valid for repulsive potentials with the same $s$-wave scattering
length Brueckner and Sawada (1957a, b); Beliaev (1958); Lieb (1963). Beyond
$E_{\text{u}}/N$ (1), the energy of the Bose gas is no more a universal
expression depending only on the $s$-wave scattering length of the potential.
Hugenholtz and Pines Hugenholtz and Pines (1959) proposed a continuation of
the series (1) in the form

$\frac{E}{N}=\frac{E_{\text{u}}}{N}+d_{1}x^{2}+d_{2}x^{5/2}\ln x+\ldots\ ,$
(5)

with new coefficients $d_{i}$ which must depend on the shape of the
interatomic potential. These parameters are expected to be related to
scattering parameters corresponding to a larger momentum transfer, essentially
the $s$-wave effective range $r_{0}$ and the $p$-wave scattering length
$a_{1}$, and also to three-body $s$-wave contributions Braaten  _et al._
(2001); Braaten and Hammer (2006); Braaten and Nieto (1997); Hammer and
Furnstahl (2000). This perturbative series for a dilute Fermi gas is well
settled up to third order of the gas parameter, the terms beyond the universal
regime being known, some of them with analytic expressions and the rest with
integrals that are numerically estimated Bishop (1973); Pera  _et al._ (2024).
In this expansion for fermions, the analytic terms depend explicitly on
$r_{0}$ and $a_{1}$. The knowledge of the nonuniversal terms for the equation
of state of a Bose gas is very reduced up to now due to difficulties in the
calculation of the required diagrams. The most recent attempt to continue the
perturbative series was carried out by Braaten Braaten  _et al._ (2001);
Braaten and Hammer (2006). In this work, an expression beyond the LHY term is
proposed with terms that correspond to three-body contributions and another
one that depends on $r_{0}$. The weight of the three-body term was determined
by a fitting procedure to diffusion Monte Carlo (DMC) results available at
that time Giorgini  _et al._ (1999). As these data was focused to low
densities, this parameter was determined with a large uncertainty.

In this paper, we study the equation of state of Bose gases beyond the
universal limit using the DMC method Boronat and Casulleras (1994). Using
different model potentials (Fig. 1), from which we know the low-energy
scattering parameters, we explore an extended universality, but now in terms
of $a_{0}$, $r_{0}$, and $a_{1}$. Our results show that this hypothesis is
approximately correct within a certain domain of $r_{0}$ values and maximum
density and that results with different interactions, sharing common
scattering values, are always closer. We discuss the limitations of the
Braaten series Braaten  _et al._ (2001) to reproduce the DMC results, mainly
when the density grows. We propose an empirical form that fits well all the
DMC data for a limited set of $|r_{0}|\lesssim 2$ values and gas parameters up
to $1\times 10^{-2}$.

Figure 1: Model potentials used in the DMC calculations.

The DMC method solves stochastically the imaginary-time Schrödinger equation
providing exact results for the ground-state energy and other properties of
the system, within some statistical noise Boronat and Casulleras (1994). To
reduce the variance one introduces a trial wave function that acts as
importance sampling, leading the sampling to regions where the exact wave
function is expected to be large and importantly eliminating divergences when
the potential becomes singular. The model that we use is a standard Jastrow
wave function $\Psi_{\rm J}({\bf R})=\prod_{i<j}^{N}f(r_{ij})$, with $f(r)$ a
two-body correlation function. We build $f(r)$ as the two-body solution,
properly symmetrized to account well with the periodic boundary conditions
(see specific details in Ref. sup ). The time-step and the average number of
configurations (walkers) is adjusted to eliminate any systematic bias in the
results. It is worth noting that our Green’s function is accurate to second
order in the time step, reducing in this way the time-step dependence and
allowing for the use of larger values which reduce the correlation length. The
model potentials used in our simulations are shown in Fig. 1 and their
characteristic values are chosen to produce the desired scattering parameters.

Our simulations work with finite numbers of particles, confined in a cubic box
of the proper size to get the expected bulk density, and thus it is important
to study any finite-size effect. This is particularly important in our case
because we need to work with the largest possible accuracy. For a fixed
density and interaction, we carry on a set of simulations with increasing
number of particles, typically in the range $N=50$–$700$. Then, we get the
energy in the thermodynamic limit $E_{0}/N$ using the linear law
$E/N=E_{0}/N+\alpha\,1/N$ sup . The time step and asymptotic population of
walkers are chosen to eliminate any significant bias coming form these two
technical parameters.

Scattering theory at low energy, and up to momentum transfer
$\leavevmode\nobreak\ k^{2}$, can be developed using three scattering
parameters: the $a_{0}$ and $a_{1}$ scattering lengths Pera and Boronat
(2023),

$a_{l}^{2l+1}=\frac{1}{2l+1}\int_{0}^{\infty}V(r)r^{l+1}u_{l}^{(0)}(r)dr\ ,$
(6)

corresponding to $l=0$ and $l=1$, respectively, and the $s$-wave effective
range

$r_{0}=\frac{2}{a_{0}^{2}}\int_{0}^{\infty}\left[(r-a_{0})^{2}-{u_{0}^{(0)}}^{2}(r)\right]dr\
.$ (7)

In Eqs. (6,7), $V(r)$ is the interaction potential, $l$ is the angular
momentum quantum number, and $u_{l}^{(0)}$ is the reduced radial wave
function, solution of the radial two-body Schrödinger equation for that
potential, with zero energy and angular momentum $l$. For the model potentials
used in our work, these three parameters are analytically known as a function
of the parameters defining the interaction Pera and Boronat (2023). The
inverse problem, that is, to know the potential parameters for a given set of
scattering parameters is not trivial due to the complex character of the
equations to solve. To this end, we implemented a stochastic search that
proved to be very efficient (see Ref. sup ).

Figure 2: Left: DMC energies using different model potentials with a common
value of the $s$-wave effective range $r_{0}$ and $s$-wave scattering length
$a_{0}$. The scattering length is used as unite length everywhere. To show the
departure from the universal regime, we subtract to the energy per particle
the universal term (1). Right: DMC energies using different potentials for
values of $|r_{0}|<2$. Same colors stand for common $r_{0}$ values.

In Fig. 2 (left panel), we show the departure from the universal regime
plotting the DMC energies using the universal energy (1), with
$c_{\text{WU}}=0$, as a reference. From here on, we do not include the
logarithmic Wu contribution because it presents a pathological behavior when
$x$ grows, even below the universal limit Giorgini  _et al._ (1999). We have
studied densities from $na^{3}=10^{-4}$ up to $10^{-2}$ and the potentials
shown in Fig. 1. Below $\leavevmode\nobreak\ 10^{-3}$, the differences between
the different calculations are tiny and not visible at the scale of the plot.
This result agrees with previous estimations obtained using also DMC Giorgini
_et al._ (1999), but now the accuracy is improved by using modern CPU’s that
allow for a more accurate determination of finite-size effects. Increasing
more the density, one can appreciate as the energies spread and, at the larger
value analyzed in our work ($10^{-2}$), the effect of finite-range effects is
clear. A first significant finding of our study is that energies appear
ordered from low to high values of $r_{0}$. To better analyze a possible
universality in the ($a$,$r_{0}$) plane, we show in Fig. 2 (right panel) a
zoom in a regime of $r_{0}$ values around one. Organized in different colors,
we show the energies obtained using different potentials that share a common
$r_{0}$ value. As one can see, at least in this window of $r_{0}$ values, the
universality in terms of the two $s$-wave scattering parameters is clearly
observed up to the largest density here explored. However, as we will discuss
later on, the range of two-parameters universality is significantly reduced
when $r_{0}$ becomes larger. The role of the $p$-wave scattering length in the
cases reported in Fig. 2 (right panel) is quantitatively below our error bars,
since in the figure we plot energies with the same $r_{0}$ but distinct
$a_{1}$.

Braaten Braaten  _et al._ (2001) proposed a continuation of the perturbative
series beyond the universal regime (1) in the form

$\displaystyle\frac{E}{N}-\frac{E_{\text{u}}}{N}$ $\displaystyle=$
$\displaystyle 4\pi
x\left[\left(c_{0}+\frac{4\pi-3\sqrt{3}}{6\pi}\ln\bar{x}\right)\bar{x}\right.$
$\displaystyle+$
$\displaystyle\left.\left(\frac{8\left(4\pi-3\sqrt{3}\right)}{3\pi^{2}}\ln\bar{x}+c_{1}-\frac{16}{15\pi}r_{0}\right){\bar{x}}^{3/2}\right]\
,$

with $\bar{x}=16\pi x$. The parameter $c_{0}$, termed $c_{E}$ in Ref. Braaten
_et al._ (2001), is the three-body contact parameter. Braaten et al. obtained
the value of $c_{0}$ by fitting Eq. (Equation of state of Bose gases beyond
the universal regime) to DMC energies of Ref. Giorgini  _et al._ (1999), but
the scarce data available and their relatively large statistical errors made
difficult its estimation. The results found were dependent on the specific
interaction, for instance for a HS potential they obtained $c_{0}=-0.9\pm 0.5$
and for a SS one $c_{0}=-1.8\pm 0.8$. It is worth noticing that our results
for $c_{0}$ are compatible with the estimations in Ref. Braaten  _et al._
(2001), within the statistical errors.

With the new data obtained in the present work we have repeated the Braaten
analysis. We selected energies for effective ranges $|r_{0}|<2$ and used both
the original model (Equation of state of Bose gases beyond the universal
regime) and the one in Eq. (Equation of state of Bose gases beyond the
universal regime) but without the logarithmic terms. In both cases, we get
statistically compatible fits ($\chi^{2}/\nu\lesssim 1$) only for gas
parameters $x\leq 4\times 10^{-3}$. With the complete Braaten model (Equation
of state of Bose gases beyond the universal regime), we obtained
$c_{0}=2.20\pm 0.06$ and $c_{1}=0.57\pm 0.16$ whereas, without the logarithmic
terms $c_{0}=-0.73\pm 0.06$ and $c_{1}=1.33\pm 0.15$, the latter showing a
slightly better $\chi^{2}/\nu$. The three-body contact parameter $c_{0}$ in
Eq. (Equation of state of Bose gases beyond the universal regime) keeps a
rather constant value $c_{0}=1.2\pm 0.1$ for energies corresponding to gas
parameters $x\leq 7\times 10^{-4}$. Its absolute value is close to the
original Braaten estimation but the sign that we obtain is changed. In Fig. 3,
we show the better Braaten equation (Equation of state of Bose gases beyond
the universal regime), restricted to the density regime where it is compatible
with our results, and the set of DMC points used for the fit (see the table of
all the DMC energies in Ref. sup ).

Table 1: Set of optimal parameters of Eq. (Equation of state of Bose gases beyond the universal regime).  $c_{0}$ | $-1.29\pm 0.07$
---|---
$c_{1}$ | $5.4\pm 0.4$
$c_{2}$ | $-9.9\pm 0.9$
$c_{3}$ | $6.2\pm 0.6$
$c_{4}$ | $2.5\pm 0.08$
$c_{5}$ | $-2.33\pm 0.12$

Figure 3: Left: Braaten equation of state (Equation of state of Bose gases
beyond the universal regime) and the set of DMC energies used for the fit.
Right: Our proposal for the equation of state (Equation of state of Bose gases
beyond the universal regime) and the DMC energies included in the fit. The
best parameters are reported in Table 1.

We extended our DMC calculations up to gas parameters $x=10^{-2}$, a value
which clearly goes beyond the validity of Eq. (Equation of state of Bose gases
beyond the universal regime). Inspired by the series expansion proposed by
Hugenholtz and Pines Hugenholtz and Pines (1959) (5) we found an empirical
equation of state in the form

$\displaystyle\frac{E}{N}-\frac{E_{\text{u}}}{N}$ $\displaystyle=$
$\displaystyle 4\pi
x\left[c_{0}\bar{x}+\left(c_{1}-\frac{16}{15\pi}r_{0}\right){\bar{x}}^{3/2}\right.$
$\displaystyle\left.+\left(c_{2}+c_{4}r_{0}\right){\bar{x}}^{2}+\left(c_{3}+c_{5}r_{0}\right){\bar{x}}^{5/2}\right]\
,$

which reproduces well our results. Compared with Eq (Equation of state of Bose
gases beyond the universal regime), we notice that we have suppressed the
logarithmic terms and increased the order of the expansion up to powers
$x^{7/2}$. We verified that we need to go this power in the series if energies
up to $x=1\times 10^{-2}$ are going to be reproduced. Regarding the
logarithmic terms, we did not observe any improvement in the law (Equation of
state of Bose gases beyond the universal regime) with their inclusion and
thus, we finally removed them. The best set of parameters are reported in
Table 1, with a value $\chi^{2}/\nu<1$. The equations of state for different
$r_{0}$ values in the range $|r_{0}|<2$ are shown, together with the DMC
energies, in Fig. 3. As one can appreciate, Eq. (Equation of state of Bose
gases beyond the universal regime) is able to reproduce quite accurately the
exact energies. For dilute Fermi gases, it is well known that series beyond
the universal terms include both the $s$-wave effective range and the $p$-wave
scattering length $a_{1}$ Bishop (1973); Pera  _et al._ (2024). In the range
of $r_{0}$ values used for the fit of Eq. (Equation of state of Bose gases
beyond the universal regime), we tried to include also the respective values
of $a_{1}$ that we also known (see Ref. sup ). We used different models to
include the $a_{1}$ values in the equation of state. However, we were not able
to discern any effect by introducing these values in the fit. This can also be
understood by looking at our results plotted in Fig. 3. There, we show results
that have in common $r_{0}$ but not $a_{1}$ and any meaningful difference is
not appreciated.

Figure 4: DMC energies for large $|r_{0}|$ not included in the fit of Eq.
(Equation of state of Bose gases beyond the universal regime). The error bars
are smaller than the size of the symbols. The solid lines stand for the
application of Eq. (Equation of state of Bose gases beyond the universal
regime) and the dashed lines incorporate the value of $a_{1}$ in the equation
of state, as commented in the text.

We carried out more DMC calculations covering larger values of the effective
range. However, we were not able to get a reliable equation of state, like Eq.
(Equation of state of Bose gases beyond the universal regime), including this
extended regime. In Fig. 4, we report the energies obtained for these
additional cases, with $r_{0}=19.7$ and $r_{0}=-4.96$. For the positive
$r_{0}$ case, we use two different potentials, HCSW and BW, and the negative
case corresponds to a SS model. The first observation on these results is the
breaking of universality, in terms of the pair $(a,r_{0})$, for the case of
large $r_{0}$ at gas parameters $x\gtrsim 3\times 10^{-3}$. Therefore, the
two-parameters universal range depends on the $|r_{0}|$ value, becoming
smaller when it grows. The case of the large $r_{0}$ is particularly
interesting. In the figure, we show three different calculations with a common
$r_{0}$ value, two of them sharing the same $p$-wave scattering length. Even
being it not perfect, we can see that the energies when the triplet
$(a,r_{0},a_{1})$ is the same are closer. This indicates that the role of
$a_{1}$, which was negligible for small $|r_{0}|$, becomes significant when
the effective range is larger.

We have analyzed if the equation of state (Equation of state of Bose gases
beyond the universal regime), fitted to a regime of small effective ranges, is
able to reproduce cases out of this domain, even in an approximate form. In
Fig. 4, we show with solid lines the predictions of Eq. (Equation of state of
Bose gases beyond the universal regime) with the parameters contained in Table
1 for the two new cases with $r_{0}=19.7$ and $r_{0}=-4.96$. The non-universal
corrections in both cases are significantly larger than the cases analyzed
before (see the different energy scales in Figs. 2 and 4) but still the
predictions are rather close to the exact results. Of course, Eq.(Equation of
state of Bose gases beyond the universal regime) can not reproduce the
dependence on $a_{1}$ observed for $r_{0}=19.7$ because this scattering
parameter is not included in the model. For the Bose equation of state, there
is not any perturbative model which includes $a_{1}$, to the best of our
knowledge. Noticeably, in the case of fermions $a_{1}$ appears explicitly in
the expansion, combined linearly with the effective range. As a first and
rather crude attempt, we summed $a_{1}^{3}$ to $r_{0}$ in the term
$\bar{x}^{3/2}$ of Eq. (Equation of state of Bose gases beyond the universal
regime) without any refitting procedure (in a second model we assumed
$r_{0}+\alpha a_{1}^{3}$ with $\alpha$ an adjusting parameter, but the
behavior of the energy with $a_{1}$ was worse than with a fixed $\alpha=1$
value). Notice that we add $a_{1}^{3}$ because the phase shift of $p$-wave
scattering is proportional to $k^{3}A_{1}$, $A_{1}$ being a volume that is
usually written as $a_{1}^{3}$ Pera and Boronat (2023). The results obtained
are shown in Fig. 4 with dashed lines. For the two small $|a_{1}|$ values, the
introduction of $a_{1}$ into the equation of state does not improve the
results. However, for the largest $|a_{1}|$ the result is surprisingly
accurate. It is also remarkable that the introduction of this scattering
parameter reproduces well the ordering of the DMC energies when $a_{1}$
changes.

Summarizing, we have explored the equation of state of dilute Bose gases by a
bunch of DMC calculations of energies, covering a large range of scattering
parameters and using different model potentials. The goal of this intensive
work was to provide a valuable set of results in order to explore the behavior
of the equation of state when universality in terms of the $s$-wave scattering
length ceases. In this regard, this is a continuation of work initiated some
years ago in Ref. Giorgini  _et al._ (1999) that allowed to determine the
approximated limit of the universality. We have carried out our calculations
up to gas parameters $x=10^{-2}$, quite far of the limit of universality
$x<10^{-3}$. Our results show sizeable deviations of the universal terms, with
both positive and negative departures. In the regime $|r_{0}|<2$, we clearly
observe universality in terms of $(a,r_{0})$, since energies calculated with
different model potentials, but sharing these two scattering parameters, are
statistically compatible up to $x\simeq 10^{-2}$. In this regime, the equation
of state proposed by Braaten Braaten  _et al._ (2001) reproduces our results
but in a smaller density domain, $x<4\times 10^{-3}$. We proposed a new
equation of state, with terms up to $x^{7/3}$, that reproduces quite
accurately the DMC energies up to $x=10^{-2}$. We used additional energies
corresponding to larger values of $r_{0}$ to check the validity of our law out
of the regime used for the fitting. Contrarily to the energies for small
effective range, we observed a sizable effect of the $p$-wave scattering
length. By fixing the set of three parameters $(a,r_{0},a_{1})$, we recovered
universality, although for not so large values of $x$. By introducing in the
equation of state the $a_{1}$ length, we improve our model, mainly for larger
values of $a_{1}$. We hope that our DMC calculations sup  can serve to
stimulate further theoretical work to improve the equation of state of Bose
gases beyond the universal limit. Also, from the experimental side, the
estimation of $r_{0}$ can serve to effectively improve the study of Bose gases
in regimes where the interactions are stronger. Finally, the introduction of
the new equation of state in the extended Gross-Pitaevskii formalism Dalfovo
_et al._ (1999) can enlarge significantly its domain of applicability.

We acknowledge Joaquim Casulleras for useful discussions to solve the inverse
scattering problem. We acknowledge financial support from Ministerio de
Ciencia e Innovación MCIN/AEI/10.13039/501100011033 (Spain) under Grant No.
PID2020-113565GB-C21 and from AGAUR-Generalitat de Catalunya Grant No.
2021-SGR-01411.

## References

  * Anderson  _et al._ (1995) M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell,  _Observation of Bose-Einstein Condensation in a Dilute Atomic Vapor_ , Science 269, 198 (1995), https://www.science.org/doi/pdf/10.1126/science.269.5221.198 .
  * Davis  _et al._ (1995) K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle,  _Bose-Einstein Condensation in a Gas of Sodium Atoms_ , Phys. Rev. Lett. 75, 3969 (1995).
  * Dalfovo  _et al._ (1999) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari,  _Theory of Bose-Einstein condensation in trapped gases_ , Rev. Mod. Phys. 71, 463 (1999).
  * Darkwah Oppong  _et al._ (2019) N. Darkwah Oppong, L. Riegger, O. Bettermann, M. Höfer, J. Levinsen, M. M. Parish, I. Bloch, and S. Fölling,  _Observation of Coherent Multiorbital Polarons in a Two-Dimensional Fermi Gas_ , Phys. Rev. Lett. 122, 193604 (2019).
  * Bombín  _et al._ (2021) R. Bombín, V. Cikojević, J. Sánchez-Baena, and J. Boronat,  _Finite-range effects in the two-dimensional repulsive Fermi polaron_ , Phys. Rev. A 103, L041302 (2021).
  * Cabrera  _et al._ (2018) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell,  _Quantum liquid droplets in a mixture of Bose-Einstein condensates_ , Science 359, 301 (2018), https://www.science.org/doi/pdf/10.1126/science.aao5686 .
  * Cikojevic  _et al._ (2020) V. Cikojevic, L. V. Markic, and J. Boronat,  _Finite-range effects in ultradilute quantum drops_ , New Journal of Physics 22, 053045 (2020).
  * Böttcher  _et al._ (2019) F. Böttcher, M. Wenzel, J.-N. Schmidt, M. Guo, T. Langen, I. Ferrier-Barbut, T. Pfau, R. Bombín, J. Sánchez-Baena, J. Boronat, and F. Mazzanti,  _Dilute dipolar quantum droplets beyond the extended Gross-Pitaevskii equation_ , Phys. Rev. Res. 1, 033088 (2019).
  * Bombín  _et al._ (2024) R. Bombín, V. Cikojević, F. Mazzanti, and J. Boronat,  _Quantum-Monte-Carlo–based functional for dysprosium dipolar systems_ , Phys. Rev. A 109, 033312 (2024).
  * Bogoliubov (1947) N. N. Bogoliubov,  _On the theory of superfluidity_ , J. Phys. (USSR) 11, 23 (1947).
  * Lee  _et al._ (1957) T. D. Lee, K. Huang, and C. N. Yang,  _Eigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature Properties_ , Phys. Rev. 106, 1135 (1957).
  * Wu (1959) T. T. Wu,  _Ground State of a Bose System of Hard Spheres_ , Phys. Rev. 115, 1390 (1959).
  * Brueckner and Sawada (1957a) K. A. Brueckner and K. Sawada,  _Bose-Einstein Gas with Repulsive Interactions: General Theory_ , Phys. Rev. 106, 1117 (1957a).
  * Brueckner and Sawada (1957b) K. A. Brueckner and K. Sawada,  _Bose-Einstein Gas with Repulsive Interactions: Hard Spheres at High Density_ , Phys. Rev. 106, 1128 (1957b).
  * Beliaev (1958) S. Beliaev,  _Energy spectrum of a non-ideal Bose gas_ , Sov. Phys. JETP 7, 299 (1958).
  * Lieb (1963) E. H. Lieb,  _Simplified Approach to the Ground-State Energy of an Imperfect Bose Gas_ , Phys. Rev. 130, 2518 (1963).
  * Hugenholtz and Pines (1959) N. M. Hugenholtz and D. Pines,  _Ground-State Energy and Excitation Spectrum of a System of Interacting Bosons_ , Phys. Rev. 116, 489 (1959).
  * Braaten  _et al._ (2001) E. Braaten, H.-W. Hammer, and S. Hermans,  _Nonuniversal effects in the homogeneous Bose gas_ , Phys. Rev. A 63, 063609 (2001).
  * Braaten and Hammer (2006) E. Braaten and H.-W. Hammer,  _Universality in few-body systems with large scattering length_ , Physics Reports 428, 259 (2006).
  * Braaten and Nieto (1997) E. Braaten and A. Nieto,  _Renormalization effects in a dilute Bose gas_ , Phys. Rev. B 55, 8090 (1997).
  * Hammer and Furnstahl (2000) H.-W. Hammer and R. Furnstahl,  _Effective field theory for dilute Fermi systems_ , Nuclear Physics A 678, 277 (2000).
  * Bishop (1973) R. Bishop,  _Ground-state energy of a dilute fermi gas_ , Annals of Physics 77, 106 (1973).
  * Pera  _et al._ (2024) J. Pera, J. Casulleras, and J. Boronat,  _Beyond universality in repulsive SU(N) Fermi gases_ (2024), arXiv:2206.06932,accepted for publication in Scipost Phys. [cond-mat.quant-gas] .
  * Giorgini  _et al._ (1999) S. Giorgini, J. Boronat, and J. Casulleras,  _Ground state of a homogeneous Bose gas: A diffusion Monte Carlo calculation_ , Phys. Rev. A 60, 5129 (1999).
  * Boronat and Casulleras (1994) J. Boronat and J. Casulleras,  _Monte Carlo analysis of an interatomic potential for He_ , Phys. Rev. B 49, 8920 (1994).
  * (26) _See Supplemental material_ , .
  * Pera and Boronat (2023) J. Pera and J. Boronat,  _Low-energy scattering parameters: A theoretical derivation of the effective range and scattering length for arbitrary angular momentum_ , American Journal of Physics 91, 90 (2023).

Supplemental Material for “Equation of state of Bose gases beyond the
universal regime”

##  S1 Trial wave functions

The DMC algorithm uses a trial or guiding wave function for importance
sampling. We use a Jastrow model,

$\Psi_{T}(\mathbf{R})=F_{J}(\mathbf{R})\Phi(\mathbf{R})\ .$ (S1)

The one-body term $\Phi(\mathbf{R})$, takes care of the symmetry or
antisymmetry of the trial wave function, that is, if particles are bosons or
fermions, respectively, as well as other global properties of the system. For
a bulk Bose system one just takes $\Phi(\mathbf{R})=1$. The second term is the
product of two-body correlation functions $f(r_{ij})$ for all pairs
$i,j=1,...,N$ in the system,

$F_{J}(\mathbf{R})=\prod_{i<j}f(r_{ij})\ .$ (S2)

For a central potential, the two-body factor $f(r_{ij})$ is solely determined
by the distance between the pair of particles. In this work, we use as a
Jastrow factor the exact solution of the two-body problem. In the following,
we discuss the explicit forms of $f(r)$ for the different model potentials
used in our study.

As our simulations are intended for a bulk system, that we model using
periodic boundary conditions, we need to specify the function $f(r)$ for any
interparticle distance $r\leq L/2$. We tried two different options: i) match
the two-body solution with a smooth function (exponential) going to a constant
at the end of that distance, and (ii) to use in all the domain the two-body
solution, properly symmetrized. We used both methods and obtained fully
compatible results. In the following, we describe the second method, which is
easier to build.

We impose that the two-body correlation function $f(r)$ satisfies

$\displaystyle f(L/2)=1,\qquad f^{\prime}(L/2)=0.$ (S3)

For practical reasons, the wave function $\mathcal{R}(r)$ needs to be
understood as an exponential, that means $\mathcal{R}(r)=e^{v(r)}$. Its
derivative, in terms of $v(r)$ is
$\mathcal{R}^{\prime}(r)=v^{\prime}(r)e^{v(r)}$ Then, the previous conditions
for $v(r)$, but already for the two-body correlation $f(r)=e^{\hat{v}(r)}$
result in

$\displaystyle e^{\hat{v}(L/2)}$ $\displaystyle=1\rightarrow\hat{v}(L/2)=0,$
(S4) $\displaystyle\hat{v}^{\prime}(L/2)e^{\hat{v}(L/2)}$
$\displaystyle=0\rightarrow\hat{v}^{\prime}(L/2)=0.$ (S5)

A solution for these two conditions is

$\hat{v}(r)=v(r)+v(L-r)-2v(L/2).$ (S6)

Actually, this is a general transformation for any function $v(r)$ to a
function $\hat{v}(r)$ that is symmetric with respect to $r=L/2$ in the
dominion $0<r<L$. Its first and second derivatives are

$\displaystyle\hat{v}^{\prime}(r)$
$\displaystyle=v^{\prime}(r)-v^{\prime}(L-r),$ (S7)
$\displaystyle\hat{v}^{\prime\prime}(r)$
$\displaystyle=v^{\prime\prime}(r)+v^{\prime\prime}(L-r).$ (S8)

One may verify that, apart from the symmetry, the latter satisfies (S4) and
(S5):

$\displaystyle\hat{v}(L/2)$ $\displaystyle=v(L/2)+v(L-L/2)-2v(L/2)=0,$ (S9)
$\displaystyle\hat{v}^{\prime}(L/2)$
$\displaystyle=v^{\prime}(L/2)-v^{\prime}(L-L/2)=0.$ (S10)

In the following, we report the detailed solution for each tested potential.

### Hard-core square-well (HCSW)

To find the solution $f(r)$, the starting point is the two-body solution
${\mathcal{R}(r)}$. In the case of the HCSW interaction potential, that is

$\mathcal{R}(r)=\begin{cases}0&0\leq r<R_{1}\\\
\frac{A}{r}\sin[k_{1}(r-R_{1})]&R_{1}\leq r<R_{2}\\\
\frac{B}{r}\sin[k_{2}(r-\delta)]&r\geq R_{2}\end{cases}$ (S11)

In terms of $v(r)=\ln{\mathcal{R}(r)}$,

$v(r)=\begin{cases}-\infty&0\leq r<R_{1}\\\
\ln{\frac{A}{r}\sin[k_{1}(r-R_{1})]}&R_{1}\leq r<R_{2}\\\
\ln{\frac{B}{r}\sin[k_{2}(r-\delta)]}&r\geq R_{2}\end{cases}$ (S12)

The parameters $A$, $B$, $\delta$ and $k_{2}$ are such that verify the
boundary conditions of both the function and its derivative. The resulting
equations are solved for $\delta$ and $k_{2}$, and $A$, $B$ are left as
variational parameters. In fact, only its relation in the form of the quotient
$A/B$ is left, as will be seen, since $A$ and $B$ are amplitudes and their
absolute value does not affect that $f(L/2)=1$. One gets

$\displaystyle A\sin[k_{1}(R_{2}-R_{1})]$
$\displaystyle=B\sin[k_{2}(R_{2}-\delta)]$ (S13) $\displaystyle
Ak_{1}\cos[k_{1}(R_{2}-R_{1})]$
$\displaystyle=Bk_{2}\cos[k_{2}(R_{2}-\delta)]$ (S14)

The resolution follows isolating $\delta$ from (S13),

$\delta=R_{2}-\frac{1}{k_{2}}\displaystyle\left\\{\arcsin{\displaystyle\left(\frac{A}{B}\sin[k_{1}(R_{2}-R_{1})]\right)}\right\\}.$
(S15)

Then, dividing (S13)/(S13), the transcendental equation for $k_{2}$ is
obtained

$\frac{\tan[k_{1}(R_{2}-R_{1})]}{k_{1}}=\frac{\tan[k_{2}(R_{2}-\delta)]}{k_{2}},$
(S16)

since $k_{1}(k_{2})=\sqrt{V_{1}+k_{2}^{2}}$ and $\delta=\delta(k_{2})$ in
(S15) (where the dependence on $A/B$ emerges). Defining $u=k_{2}$, the latter
is transformed into the zero of the function

$f(u)=\frac{\tan[k_{1}(u)(R_{2}-R_{1})]}{k_{1}(u)}-\frac{\tan[u(R_{2}-\delta(u))]}{u},$
(S17)

which can be calculated using the Newton-Raphson method, using that its
derivative is

$\displaystyle f^{\prime}(u)=$
$\displaystyle\frac{k_{1}^{\prime}(u)(R_{2}-R_{1})}{k_{1}(u)\cos^{2}[k_{1}(u)(R_{2}-R_{1})]}-\frac{k_{1}^{\prime}(u)\tan[k_{1}(u)(R_{2}-R_{1})]}{k_{1}(u)^{2}}-$
$\displaystyle-\frac{R_{2}-\delta(u)-u\delta^{\prime}(u)}{u\cos^{2}[u(R_{2}-\delta(u))]}+\frac{\tan[u(R_{2}-\delta(u))]}{u^{2}}$

where

$\displaystyle k_{1}^{\prime}(u)$
$\displaystyle=\frac{2u}{2\sqrt{V_{1}+u^{2}}}=\frac{u}{k_{1}(u)},$ (S18)
$\displaystyle\delta^{\prime}(u)$
$\displaystyle=\frac{\arcsin{\displaystyle\left(\frac{A}{B}\sin[k_{1}(u)(R_{2}-R_{1})]\right)}}{u^{2}}-\frac{\frac{A}{B}\cos[k_{1}(u)(R_{2}-R_{1})]k_{1}^{\prime}(u)(R_{2}-R_{1})}{u\sqrt{1-\displaystyle\left(\frac{A}{B}\sin[k_{1}(u)(R_{2}-R_{1})]\right)^{2}}}.$
(S19)

Once solved for $u=k_{2}$, $\delta$ is obtained from (S15). Thus, the analytic
expressions of $\mathcal{R}(r)$ and $v(r)$ are obtained.

From now on, the notation $v_{1}(r)$, $v_{2}(r)$ is employed to refer to the
solution $v(r)$ in the region $R_{1}\leq r<R_{2}$ and $r\geq R_{2}$
respectively. Since $v(r)$ is a piece-wise function, $\hat{v}(r)$ also is. To
construct it, (S6) is applied for each region.

$\hat{v}(r)=\begin{cases}-\infty+v_{2}(L-r)-2v_{2}(L/2)=-\infty&0\leq
r<R_{1}\\\ \hat{v_{1}}(r)=v_{1}(r)+v_{2}(L-r)-2v_{2}(L/2)&R_{1}\leq r<R_{2}\\\
\hat{v_{2}}(r)=v_{2}(r)+v_{2}(L-r)-2v_{2}(L/2)&R_{2}\leq r\leq L/2\end{cases}$
(S20)

Note that the term $v(L-r)$ always corresponds to the expression $v_{2}(L-r)$,
because $L-r>R_{2}$ for $0\leq r\leq L/2$ given that $R_{2}<L/2$. Also,
$v(L/2)$ corresponds to $v_{2}(L/2)$ for the same reason. The explicit
expressions for the new wave functions are

$\displaystyle\hat{v_{1}}(r)=$
$\displaystyle\ln{\frac{A}{r}\sin[k_{1}(r-R_{1})]}+\ln{\frac{B}{L-r}\sin[k_{2}(L-r-\delta)]}-2\ln{\frac{B}{L/2}\sin[k_{2}(L/2-\delta)]}$
$\displaystyle\hat{v_{2}}(r)=$
$\displaystyle\ln{\frac{B}{r}\sin[k_{2}(r-\delta)]}+\ln{\frac{B}{L-r}\sin[k_{2}(L-r-\delta)]}-2\ln{\frac{B}{L/2}\sin[k_{2}(L/2-\delta)]}$

and $f(r)=e^{\hat{v}(r)}$.

### Barrier-well (BW)

The two-body solution is

$\mathcal{R}(r)=\begin{cases}\frac{A}{r}\sinh[k_{1}r]&0\leq r<R_{1}\\\
\frac{B}{r}\sin[k_{2}(r-\delta)]&R_{1}\leq r<R_{2}\\\
\frac{C}{r}\sin[k_{3}(r-\xi)]&r\geq R_{2}\\\ \end{cases}$ (S21)

where $k_{1}=\sqrt{V_{1}-k_{3}^{2}}$ and $k_{2}=\sqrt{V_{2}+k_{3}^{2}}$. By
imposing the right boundary conditions, one gets

$\displaystyle A\sinh[k_{1}R_{1}]$ $\displaystyle=B\sin[k_{2}(R_{1}-\delta)]$
(S22) $\displaystyle Ak_{1}\cosh[k_{1}R_{1}]$
$\displaystyle=Bk_{2}\cos[k_{2}(R_{2}-\delta)]$ (S23) $\displaystyle
B\sin[k_{2}(R_{2}-\delta)]$ $\displaystyle=C\sin[k_{3}(R_{2}-\xi)]$ (S24)
$\displaystyle Bk_{2}\cos[k_{2}(R_{2}-\delta)]$
$\displaystyle=Ck_{3}\cos[k_{3}(R_{2}-\xi)]$ (S25)

This system of equations is solved for $C$, $\delta$, $\xi$, and $k_{3}$.
Isolating $\delta$ from (S22),

$\delta=R_{1}-\frac{1}{k_{2}}\displaystyle\left\\{\arcsin{\displaystyle\left(\frac{A}{B}\sinh[k_{1}R_{1}]\right)}\right\\}$
(S26)

and combining (S22)/(S23) and (S24)/(S25), the system is reduced to

$\displaystyle\frac{\tanh[k_{1}R_{1}]}{k_{1}}$
$\displaystyle=\frac{\tan[k_{2}(R_{1}-\delta)]}{k_{2}}$ (S27)
$\displaystyle\frac{\tan[k_{2}(R_{2}-\delta)]}{k_{2}}$
$\displaystyle=\frac{\tan[k_{3}(R_{2}-\xi)]}{k_{3}}$ (S28)

Given its form, setting $B=1$, and leaving $A$ as a variational parameter,
conveniently adjusted, the solution for $k_{3}$, $\xi$ is given by the zero of
the multivariate function

$f(k_{3},\xi)=\displaystyle\left(\frac{\tanh[k_{1}R_{1}]}{k_{1}}-\frac{\tan[k_{2}(R_{1}-\delta)]}{k_{2}},\frac{\tan[k_{2}(R_{2}-\delta)]}{k_{2}}-\frac{\tan[k_{3}(R_{2}-\xi)]}{k_{3}}\right).$
(S29)

Then, $\delta=\delta(k_{3})$ is recovered from (S26) and $C$ is determined
from (S24),

$C=\frac{B\sin[k_{2}(R_{2}-\delta)]}{\sin[k_{3}(R_{2}-\xi)]}.$ (S30)

The resulting transcendental equation is solved by means of a multidimensional
Newton-Raphson method. However, this time the derivative
$f^{\prime}(k_{3},\xi)$ is approximated numerically, due to the complexity of
its analytical expression. To explore its solution, $f(k_{3},\xi)$ is
understood as two surfaces in the $k_{3}$-$\xi$ plane, given by its first and
second components. Then, the initial guess and $A$ values are set provided
that both surfaces intersect with the plane $z=0$, meaning the solution
exists. From the two-body solution $\mathcal{R}(r)=e^{v(r)}$, one only needs
to combine it forming

$\hat{v}(r)=\begin{cases}\hat{v_{1}}(r)=v_{1}(r)+v_{3}(L-r)-2v_{3}(L/2)&0\leq
r<R_{1}\\\ \hat{v_{2}}(r)=v_{2}(r)+v_{3}(L-r)-2v_{3}(L/2)&R_{1}\leq r<R_{2}\\\
\hat{v_{3}}(r)=v_{3}(r)+v_{3}(L-r)-2v_{3}(L/2)&R_{2}\leq r\leq L/2\end{cases}$
(S31)

Again, the term $v(L-r)$ and $v(L/2)$ corresponds to expression for $v(r)$
that holds in the region $R_{2}\leq r\leq L/2$, since $L-r>R_{2}$ for $0\leq
r\leq L/2$ given that $R_{2}<L/2$. The explicit expressions are

$\displaystyle\hat{v_{1}}(r)=$
$\displaystyle\ln{\frac{A}{r}\sin[k_{1}r]}+\ln{\frac{C}{L-r}\sin[k_{3}(L-r-\xi)]}-2\ln{\frac{C}{L/2}\sin[k_{3}(L/2-\xi)]}$
$\displaystyle\hat{v_{2}}(r)=$
$\displaystyle\ln{\frac{B}{r}\sin[k_{2}(r-\delta)]}+\ln{\frac{C}{L-r}\sin[k_{3}(L-r-\xi)]}-2\ln{\frac{C}{L/2}\sin[k_{3}(L/2-\xi)]}$
$\displaystyle\hat{v_{3}}(r)=$
$\displaystyle\ln{\frac{C}{r}\sin[k_{3}(r-\xi)]}+\ln{\frac{C}{L-r}\sin[k_{3}(L-r-\xi)]}-2\ln{\frac{C}{L/2}\sin[k_{3}(L/2-\xi)]}$

### Soft sphere (SS)

The two-body solution is readily obtained,

$\mathcal{R}(r)=\begin{cases}\frac{A}{r}\sinh[k_{1}r]&0\leq r<R_{1}\\\
\frac{B}{r}\sin[k_{2}(r-\delta)]&r>R_{1}\\\ \end{cases}$ (S32)

where $k_{2}=\sqrt{mE}/\hbar$ and $k_{1}=\sqrt{mV_{1}/\hbar^{2}-k_{2}^{2}}$.
Imposing boundary conditions,

$\displaystyle A\sinh[k_{1}R_{1}]$ $\displaystyle=B\sin[k_{2}(R_{1}-\delta)]$
(S33) $\displaystyle Ak_{1}\cosh[k_{1}R_{1}]$
$\displaystyle=Bk_{2}\cos[k_{2}(R_{1}-\delta)]$ (S34)

From (S33), $\delta$ can be isolated

$\delta=R_{1}-\frac{1}{k_{2}}\displaystyle\left\\{\arcsin{\displaystyle\left(\frac{A}{B}\sinh[k_{1}R_{1}]\right)}\right\\}.$
(S35)

Then, as in the previous cases, the solution only depends on the ratio $A/B$,
which is left as a variational parameter that will be properly adjusted.
Dividing (S33)/(S34), the transcendental equation for $k_{2}$ is obtained,

$\frac{\tanh[k_{1}R_{1}]}{k_{1}}=\frac{\tan[k_{2}(R_{1}-\delta)]}{k_{2}}.$
(S36)

Then, defining $u=k_{2}$, it becomes the zero of the function

$f(u)=\frac{\tanh[k_{1}(u)R_{1}]}{k_{1}(u)}-\frac{\tan[u(R_{1}-\delta(u))]}{u},$
(S37)

which can be calculated using the Newton-Raphson method, knowing that its
derivative is

$\displaystyle f^{\prime}(u)=$
$\displaystyle\frac{k_{1}^{\prime}(u)R_{1}}{k_{1}(u)\cosh^{2}[k_{1}(u)R_{1}]}-\frac{k_{1}^{\prime}(u)\tanh[k_{1}(u)R_{1}]}{k_{1}(u)^{2}}-$
$\displaystyle-\frac{R_{1}-\delta(u)-u\delta^{\prime}(u)}{u\cos^{2}[u(R_{1}-\delta(u))]}+\frac{\tan[u(R_{1}-\delta(u))]}{u^{2}}$

where

$\displaystyle k_{1}^{\prime}(u)$
$\displaystyle=\frac{-2u}{2\sqrt{mV_{1}/\hbar^{2}-u^{2}}}=-\frac{u}{k_{1}(u)}$
(S38) $\displaystyle\delta^{\prime}(u)$
$\displaystyle=\frac{\arcsin{\displaystyle\left(\frac{A}{B}\sinh[k_{1}(u)R_{1}]\right)}}{u^{2}}-\frac{\frac{A}{B}\cosh[k_{1}(u)R_{1}]k_{1}^{\prime}(u)R_{1}}{u\sqrt{1-\displaystyle\left(\frac{A}{B}\sinh[k_{1}(u)R_{1}]\right)^{2}}}$
(S39)

Once solved for $u=k_{2}$, $\delta$ is obtained from (S35). Thus, the analytic
expressions of $\mathcal{R}(r)$ and $v(r)$ are known. Then, the expressions
for $\hat{v}(r)$ are

$\hat{v}(r)=\begin{cases}\hat{v_{1}}(r)=v_{1}(r)+v_{2}(L-r)-2v_{2}(L/2)&0\leq
r<R_{1}\\\ \hat{v_{2}}(r)=v_{2}(r)+v_{2}(L-r)-2v_{2}(L/2)&r\geq R_{1}\\\
\end{cases}$ (S40)

following the same reasoning as with the previous potentials.

##  S2 Solution of the inverse problem

Given an interaction potential $V(r)$, we are interested in its scattering
parameters $r_{0}$ ($s$-wave effective range), $a_{0}$ (s-wave scattering
length) and $a_{1}$ (p-wave scattering length). These are given by

$\displaystyle a_{l}^{2l+1}$
$\displaystyle=\frac{1}{2l+1}\int_{0}^{\infty}V(r)r^{l+1}u_{l}^{(0)}(r)dr,$
(S41) $\displaystyle r_{0}^{eff}$
$\displaystyle=\frac{2}{a_{0}^{2}}\int_{0}^{\infty}\displaystyle\left[(r-a_{0})^{2}-{u_{0}^{(0)}}^{2}(r)\right]dr,$
(S42)

where $l$ is the angular momentum quantum number and has been already
substituted by $l=0$ in Eq. (S42). The function $u_{l}^{(0)}(r)$ is the
reduced radial wave function, solution to the radial two-body Schrödinger
equation for that potential, with zero energy and angular momentum $l$. It
must verify the boundary condition $u_{l}^{(0)}(r=0)=0$ (the wave function
must vanish at the origin) and be properly normalized to behave as
$r^{l+1}-a_{l}^{2l+1}r^{-l}$ when $r\rightarrow\infty$.

Our interest lies in comparing the energy of systems described by different
interaction potentials, each of them with its scattering parameters. To do so,
one needs to be able to calculate the scattering parameters of the interaction
potentials employed. We refer to Ref. Pera and Boronat (2023) for their
explicit analytical expressions, which result from solving Eqs. ((S41),(S42))
for each potential. All the interaction potential models employed are simple,
piece-wise functions and can be tuned through a small number of parameters. As
an example, if the potential model presents a finite barrier, one can obtain
different potentials by changing its width (distance) and height (potential).
The referred equations explicitly relate the model’s parameters with the
scattering parameters. Thus, theoretically one can fix $n$ scattering
parameters for a potential model that presents $n$ tuning options. This is
exactly what we are interested in: to select any value of the scattering
parameters and find the exact different potentials that lead to them, to study
the dependencies at wish. The problem here lies in the fact that this
association leads to a non-trivial system of equations.

As an example, consider the barrier-well potential. This model consists of a
finite barrier and a well around it. The number of tuning options is four:
barrier length $R_{1}$, well length $R_{2}$, barrier height $V_{1}$ and well
height $V_{2}$. Then, one could solve a system of equations for up to four
scattering parameters. Since we are only interested in three of them, the
system of equations is given by: $a_{0}(R_{1},R_{2},V_{1},V_{2})$,
$r_{0}(R_{1},R_{2},V_{1},V_{2})$, $a_{1}(R_{1},R_{2},V_{1},V_{2})$. The
expressions can be found in Pera and Boronat (2023), by switching $k_{1}$,
$k_{2}$ to $ik_{1}$, $ik_{2}$ respectively from the well-barrier equations.
This can be done because the barrier-well model is the inverse of the well-
barrier, meaning the only difference is a change of sign of the potential
values, which affects the equations by turning $k_{1}$ and $k_{2}$ (the square
root of the potential values) imaginary. Then, one could attempt to fix a
triad $a_{0}$, $r_{0}$, $a_{1}$ and $R_{1}$ and find the values of $R_{2}$,
$V_{1}$ and $V_{2}$ that verify the system of equations. It is not possible to
isolate analytically the unknowns from the system of equations, so a numerical
method is the only way to proceed. A usual approach is the multidimensional
Newton-Raphson method. However, only that does not solve the problem in
general; previously it is necessary that one: (i) ensures that there exists a
solution for the set of parameters fixed and (ii) provides an initial guess
for the unknown variables that is close to the solution of the system of
equations.

We use a stochastic strategy to explore the scattering equations, by solving
the inverse problem. That means, trying different combinations
$\\{R_{1},R_{2},V_{1},V_{2}\\}$ and calculating the scattering parameters for
each set, employing the equations provided. Then, each combination leads to a
point in the $(a_{0},a_{1},r_{0})$ space. Repeating this process for random
combinations (all of them feasible for that potential model) provides an idea
of the images of the scattering equations. In the stochastic sense, each
combination is a choice and provides a triad which is solution of the
scattering parameters. Then, the regions with more density of points are those
where it is more likely to find a triad solution. Conversely, the regions with
no points are not likely reporting a potential based on that model with those
scattering parameters. By inspecting the whole picture, one has an idea about
what scattering parameters’ triads can be obtained, and, through a simple
algorithm, recovers the parameters of the potential model that lead to that
triad. Finally, those potential parameters serve as the initial condition for
the Newton-Raphson method to find triads similar to the one that is obtained
by them, up to the desired accuracy.

This strategy has an even major interest in our scope: it allows us to compare
the images of the scattering parameters’ equations of different potential
models (Fig. S1). Thus, it can be extremely useful to find coincident triads
for different potential models and, if possible, which are approximately its
values. This case corresponds to a region where two points from different
potentials are close in distance.

(a) Zoom out.

(b) Zoom in.

Figure S1: Space of solutions $(a_{0},a_{1},r_{0})$ of the scattering
equations for different potentials (see legend) for $10^{6}$ points.

(a) Zoom out.

(b) Zoom in.

Figure S2: Slice $a_{0}\approx 1$ in the space of solutions of the scattering
equations for different potentials (see legend) for $10^{8}$ points.

We choose $a_{0}=1$ for all the interaction potentials, such that the distance
is in units of $a_{0}$ in the DMC calculations. Hence, we are interested in
the points contained in the plane $a_{0}=1$. In fact, Fig. S2 shows the points
contained in $a_{0}\in[1-\delta,1+\delta]$. Recall that due to the stochastic
nature of this strategy, it is not possible to obtain any point with $a_{0}=1$
exactly, thus we select those that are sufficiently close such that they serve
as an initial condition for the Newton-Raphson method to find a solution with
$a_{0}=1$ exactly.

Once a solution triad is found, with its associated potential model
parameters, we use them to calculate the scattering parameters differently
again (solving the integral forms (S41), (S42)), such that if they coincide
the solution is validated. Also, it is checked that a bound state does not
appear, the $a_{0}$ value of the triad solution being always at the left of
the first Feshbach resonance.

##  S3 Correction of finite-size effects in DMC

A DMC computation is set for a specific number of particles $N$ at a given
density $n$, defining a cubic box of side $L=\sqrt[3]{N/n}$. Obviously, $N$ is
limited by the increase of computer time, that grows as $\sim N^{2}$. We need
to estimate finite-size effects because our calculations are intended for a
bulk phase in the thermodynamic limit.

To overcome this issue, we perform several DMC calculations with different
numbers of particles and, as a check, we compare our results with $E_{u}/N$ at
a density $na_{0}^{3}=10^{-5}$ (Fig. S3) where universality holds. A similar
analysis is shown in Fig. S4 where $E_{u}/N$ is no longer valid. Notice that
for efficiency reasons in the DMC branching process, the product $E\,\Delta t$
needs to remain approximately constant, hence the time-step $\Delta t$ should
be reduced when $N$, and thus $E$, increases. As an approximation, we deal
with this by setting the product $N\,\Delta t$ constant along DMC calculations
with different numbers of particles, at each density and interatomic
potential.

We observe that the difference between the energy per particle in the DMC
calculation $E/N$ and the universal energy $E_{u}/N$ is linear with the
inverse number of particles $1/N$. At low density, this statement holds in the
light of Figure S3(a); while for larger density, Fig. S4(a) the linear law is
less accurate. As one can see, at low density the agreement between the
extrapolation to $1/N\to 0$ and the expected result is perfect, improving
slightly previous DMC estimations on the same system Giorgini  _et al._
(1999).

Based on the check at low density, all our DMC results consist of calculating
the energy with $N=150$ and $N=300$, and then extrapolate linearly to the
limit $1/N\to 0$. The errors are safely propagated accordingly. At low density
(Fig. S3(b)), this procedure leads to reproducing $E_{u}/N$.

(a) Linear fit to the DMC energies in terms of $1/N$.

(b) Extrapolation using data obtained with $N=150$ and 300\.

Figure S3: DMC energies for different numbers of particles interacting through
a Hard-Sphere potential at density $na_{0}^{3}=10^{-5}$.

(a) Linear fit to the DMC energies in terms of $1/N$.

(b) Extrapolation using data obtained with $N=150$ and 300\.

Figure S4: DMC energies for different numbers of particles interacting through
a Soft-Sphere potential (with scattering parameters $r_{0}=-4.96$,
$a_{1}=1.83$) at density $na^{3}=10^{-3}$.

##  S4 DMC results

Model | $a_{0}$ | $r_{0}$ | $a_{1}$ | $V_{1}$ | $V_{2}$ | $R_{1}$ | $R_{2}$
---|---|---|---|---|---|---|---
SS | $1$ | $-4.963717$ | $1.828551$ | $3.154281E-2$ | – | 5 | –
SS | $1$ | $0.246789895603954$ | $1.192367656777861$ | $0.553699098716986$ | – | $2.25305223642389$ | –
SS | $1$ | $0$ | $1.257306682690198$ | $0.352264068783665$ | – | $2.52455873378522$ | –
HCSW | $1$ | $19.7493$ | $-2.44558$ | $2.816293678E-2$ | – | $2.23871367134203$ | $6.87815732919500$
HCSW | $1$ | $19.7493$ | $-3.32211$ | $8.879685093760E-2$ | – | $2.65721897814958$ | $5.91964724378216$
HCSW | $1$ | $0.951316$ | $-1.6636$ | $0.212152394371251$ | – | $1.24099358420681$ | $2.65518281462878$
HCSW | $1$ | $0.951316$ | $2.00307$ | $7.4019719854E-3$ | – | $1.07208359841336$ | $4.12341010567992$
HCSW | $1$ | $0.6842384$ | $1.05598013574085$ | $3.54339300812E-3$ | – | $1.01028825874101$ | $3.06371922192216$
HCSW | $1$ | $0.6842384$ | $-1.4545$ | $7.14004450115E-2$ | – | $1.03334693697237$ | $2.13914846375276$
HCSW | $1$ | $1.8227872$ | $0.932469981150490$ | $3.82475925030E-2$ | – | $1.29735440984730$ | $4.04020841908368$
HCSW | $1$ | $0.7$ | $0$ | $1.01776107899E-2$ | – | $1.02287211424451$ | $2.90289035827995$
BW | $1$ | $19.7493$ | $-2.44558$ | $0.393919355038301$ | $2.088440731E-2$ | $3.73355462692938$ | $7.312800$
BW | $1$ | $0.6842384$ | $1.05598013574085$ | $0.8005690524$ | $0.101545616787338$ | $2.28829062303712$ | $3.09329071235972$
BW | $1$ | $1.192368$ | $0.246789895603954$ | $0.55130064452$ | $5.7315830258E-3$ | $2.26041796351839$ | $2.62247347275118$
Table 2: Parameters of the two-body interaction potentials. | HCSW | HCSW | HCSW | HCSW
---|---|---|---|---
$na_{0}^{3}$ | $r_{0}=0.95$ | $r_{0}=0.95$ | $r_{0}=0.68$ | $r_{0}=0.68$
| $a_{1}=-1.66$ | $a_{1}=2.00$ | $a_{1}=1.05$ | $a_{1}=-1.45$
$1.0E-5$ | $1.2633E-4\pm 8.7E-7$ | $1.266E-4\pm 1.1E-6$ | $1.2695E-4\pm 6.5E-7$ | $1.2683E-4\pm 8.5E-7$
$1.0E-4$ | $1.306E-3\pm 1.1E-5$ | $1.307E-3\pm 1.2E-5$ | $1.306E-3\pm 1.2E-5$ | $1.307E-3\pm 1.1E-5$
$2.0E-4$ | $2.667E-3\pm 1.7E-5$ | $2.658E-3\pm 2.0E-5$ | $2.658E-3\pm 2.0E-5$ | $2.654E-3\pm 1.7E-5$
$3.0E-4$ | $4.050E-3\pm 3.2E-5$ | $4.040E-3\pm 2.4E-5$ | $4.037E-3\pm 2.6E-5$ | $4.038E-3\pm 2.6E-5$
$5.0E-4$ | $6.857E-3\pm 6.5E-5$ | $6.867E-3\pm 6.8E-5$ | $6.862E-3\pm 5.8E-5$ | $6.865E-3\pm 6.1E-5$
$7.0E-4$ | $9.773E-3\pm 7.0E-5$ | $9.764E-3\pm 7.5E-5$ | $9.756E-3\pm 6.6E-5$ | $9.758E-3\pm 6.2E-5$
$1.0E-3$ | $1.424E-2\pm 1.4E-4$ | $1.422E-2\pm 1.2E-4$ | $1.420E-2\pm 1.2E-4$ | $1.420E-2\pm 1.2E-4$
$2.0E-3$ | $2.989E-2\pm 2.6E-4$ | $2.984E-2\pm 2.9E-4$ | $2.976E-2\pm 2.1E-4$ | $2.977E-2\pm 1.8E-4$
$3.0E-3$ | $4.650E-2\pm 3.0E-4$ | $4.641E-2\pm 3.2E-4$ | $4.623E-2\pm 4.4E-4$ | $4.622E-2\pm 3.5E-4$
$4.0E-3$ | $6.375E-2\pm 3.5E-4$ | $6.382E-2\pm 3.3E-4$ | $6.339E-2\pm 4.5E-4$ | $6.342E-2\pm 4.6E-4$
$5.0E-3$ | $8.211E-2\pm 5.3E-4$ | $8.189E-2\pm 6.5E-4$ | $8.127E-2\pm 7.4E-4$ | $8.131E-2\pm 4.8E-4$
$6.0E-3$ | $1.0086E-1\pm 6.1E-4$ | $1.0062E-1\pm 9.4E-4$ | $9.972E-2\pm 6.3E-4$ | $9.971E-2\pm 6.8E-4$
$7.0E-3$ | $1.2040E-1\pm 7.1E-4$ | $1.202E-1\pm 1.1E-3$ | $1.1878E-1\pm 6.5E-4$ | $1.188E-1\pm 1.1E-3$
$8.0E-3$ | $1.407E-1\pm 1.3E-3$ | $1.4008E-1\pm 9.0E-4$ | $1.3844E-1\pm 9.5E-4$ | $1.3837E-1\pm 7.6E-4$
$9.0E-3$ | $1.613E-1\pm 1.3E-3$ | $1.608E-1\pm 1.2E-3$ | $1.5844E-1\pm 9.8E-4$ | $1.585E-1\pm 1.5E-3$
$1.0E-2$ | $1.815E-1\pm 1.7E-3$ | $1.813E-1\pm 1.2E-3$ | $1.790E-1\pm 1.2E-3$ | $1.790E-1\pm 1.7E-3$
Table 3: DMC energies. | BW | BW | SS | HCSW
---|---|---|---|---
$na_{0}^{3}$ | $r_{0}=0.68$ | $r_{0}=0.25$ | $r_{0}=0.25$ | $r_{0}=1.82$
| $a_{1}=1.05$ | $a_{1}=1.19$ | $a_{1}=1.19$ | $a_{1}=0.93$
$1.0E-5$ | $1.2671E-4\pm 9.9E-7$ | $1.2684E-4\pm 9.6E-7$ | $1.269E-4\pm 1.1E-6$ | $1.257E-4\pm 1.2E-6$
$1.0E-4$ | $1.303E-3\pm 1.2E-5$ | $1.308E-3\pm 1.1E-5$ | $1.307E-3\pm 1.3E-5$ | $1.308E-3\pm 1.0E-5$
$2.0E-4$ | $2.649E-3\pm 1.7E-5$ | $2.656E-3\pm 2.1E-5$ | $2.657E-3\pm 2.2E-5$ | $2.660E-3\pm 1.4E-5$
$3.0E-4$ | $4.032E-3\pm 3.1E-5$ | $4.030E-3\pm 3.5E-5$ | $4.033E-3\pm 3.7E-5$ | $4.045E-3\pm 4.0E-5$
$5.0E-4$ | $6.836E-3\pm 5.7E-5$ | $6.847E-3\pm 5.9E-5$ | $6.847E-3\pm 6.0E-5$ | $6.882E-3\pm 4.3E-5$
$7.0E-4$ | $9.716E-3\pm 5.5E-5$ | $9.731E-3\pm 6.4E-5$ | $9.728E-3\pm 6.3E-5$ | $9.802E-3\pm 5.7E-5$
$1.0E-3$ | $1.403E-2\pm 1.4E-4$ | $1.416E-2\pm 1.3E-5$ | $1.4157E-2\pm 7.5E-5$ | $1.428E-2\pm 1.1E-4$
$2.0E-3$ | $2.985E-2\pm 1.6E-4$ | $2.961E-2\pm 1.6E-4$ | $2.964E-2\pm 2.5E-4$ | $3.014E-2\pm 2.8E-4$
$3.0E-3$ | $4.609E-2\pm 3.3E-4$ | $4.579E-2\pm 3.6E-4$ | $4.585E-2\pm 4.1E-4$ | $4.721E-2\pm 2.9E-4$
$4.0E-3$ | $6.331E-2\pm 3.5E-4$ | $6.263E-2\pm 3.5E-4$ | $6.267E-2\pm 3.1E-4$ | $6.518E-2\pm 5.3E-4$
$5.0E-3$ | $8.122E-2\pm 7.3E-4$ | $8.005E-2\pm 4.4E-4$ | $7.985E-2\pm 5.8E-4$ | $8.413E-2\pm 7.6E-4$
$6.0E-3$ | $9.982E-2\pm 7.7E-4$ | $9.785E-2\pm 8.2E-4$ | $9.794E-2\pm 9.7E-4$ | $1.0416E-1\pm 9.3E-4$
$7.0E-3$ | $1.1839E-1\pm 9.0E-4$ | $1.1618E-1\pm 6.8E-4$ | $1.158E-1\pm 1.1E-3$ | $1.2484E-1\pm 7.7E-4$
$8.0E-3$ | $1.3761E-1\pm 9.5E-4$ | $1.3465E-1\pm 1.2E-3$ | $1.347E-1\pm 1.3E-3$ | $1.4646E-1\pm 8.8E-4$
$9.0E-3$ | $1.5730E-1\pm 8.9E-4$ | $1.5386E-1\pm 1.5E-3$ | $1.538E-1\pm 1.4E-3$ | $1.691E-1\pm 1.5E-3$
$1.0E-2$ | $1.775E-1\pm 1.6E-3$ | $1.730E-1\pm 1.6E-3$ | $1.729E-1\pm 1.7E-3$ | $1.930E-1\pm 1.9E-3$
Table 4: DMC energies. | HCSW | SS | HCSW | BW
---|---|---|---|---
$na_{0}^{3}$ | $r_{0}=0.7$ | $r_{0}=0$ | $r_{0}=19.7493$ | $r_{0}=19.7493$
| $a_{1}=0$ | $a_{1}=1.2573$ | $a_{1}=-3.32211$ | $a_{1}=-2.44558$
$1.0E-5$ | $1.268E-4\pm 1.2E-6$ | $1.269E-4\pm 1.1E-6$ | $1.3092E-4\pm 7.2E-7$ | $1.270E-4\pm 1.1E-6$
$1.0E-4$ | $1.306E-3\pm 1.2E-5$ | $1.309E-3\pm 1.3E-5$ | $1.3244E-3\pm 8.5E-6$ | $1.3363E-3\pm 9.2E-6$
$2.0E-4$ | $2.658E-3\pm 1.4E-5$ | $2.642E-3\pm 2.4E-5$ | $2.742E-3\pm 1.7E-5$ | $2.688E-3\pm 2.2E-5$
$3.0E-4$ | $4.038E-3\pm 3.3E-5$ | $4.034E-3\pm 2.5E-5$ | $4.290E-3\pm 3.8E-5$ | $4.154E-3\pm 2.3E-5$
$5.0E-4$ | $6.863E-3\pm 6.2E-5$ | $6.846E-3\pm 4.0E-5$ | $7.569E-3\pm 4.5E-5$ | $7.295E-3\pm 4.7E-5$
$7.0E-4$ | $9.757E-3\pm 5.6E-5$ | $9.734E-3\pm 7.4E-5$ | $1.124E-2\pm 1.1E-4$ | $1.0536E-2\pm 7.5E-5$
$1.0E-3$ | $1.4201E-2\pm 7.1E-5$ | $1.4149E-2\pm 9.7E-5$ | $1.756E-2\pm 1.5E-4$ | $1.5733E-2\pm 8.7E-5$
$2.0E-3$ | $2.978E-2\pm 2.0E-4$ | $2.955E-2\pm 1.6E-4$ | $4.605E-2\pm 4.1E-4$ | $3.595E-2\pm 3.1E-4$
$3.0E-3$ | $4.625E-2\pm 3.5E-4$ | $4.567E-2\pm 2.6E-4$ | $8.922E-2\pm 6.2E-4$ | $6.083E-2\pm 6.0E-4$
$4.0E-3$ | $6.342E-2\pm 4.7E-4$ | $6.231E-2\pm 3.9E-4$ | $1.515E-1\pm 1.3E-3$ | $8.939E-2\pm 6.0E-4$
$5.0E-3$ | $8.136E-2\pm 5.6E-4$ | $7.936E-2\pm 4.6E-4$ | $2.369E-1\pm 1.7E-3$ | $1.2179E-1\pm 7.6E-4$
$6.0E-3$ | $9.977E-2\pm 8.0E-4$ | $9.681E-2\pm 5.9E-4$ | $3.512E-1\pm 2.3E-3$ | $1.568E-1\pm 1.5E-3$
$7.0E-3$ | $1.190E-1\pm 1.1E-3$ | $1.146E-1\pm 1.0E-3$ | $5.017E-1\pm 2.6E-3$ | $1.950E-1\pm 1.0E-3$
$8.0E-3$ | $1.3850E-1\pm 7.4E-4$ | $1.330E-1\pm 1.2E-3$ | $6.886E-1\pm 3.8E-3$ | $2.357E-1\pm 1.7E-3$
$9.0E-3$ | $1.586E-1\pm 1.3E-3$ | $1.5155E-1\pm 8.5E-4$ | $9.284E-1\pm 6.2E-3$ | $2.763E-1\pm 1.7E-3$
$1.0E-2$ | $1.791E-1\pm 1.3E-3$ | $1.703E-1\pm 1.4E-3$ | $1.2158E-0\pm 7.3E-3$ | $3.204E-1\pm 1.9E-3$
Table 5: DMC energies. | HCSW | SS
---|---|---
$na_{0}^{3}$ | $r_{0}=19.7493$ | $r_{0}=-4.96$
| $a_{1}=-2.44558$ | $a_{1}=1.83$
$1.0E-5$ | $1.2717E-4\pm 9.2E-7$ | $1.2712E-4\pm 6.8E-7$
$1.0E-4$ | $1.329E-3\pm 1.2E-5$ | $1.304E-3\pm 1.1E-5$
$2.0E-4$ | $2.724E-3\pm 2.5E-5$ | $2.650E-3\pm 2.4E-5$
$3.0E-4$ | $4.209E-3\pm 3.8E-5$ | $4.025E-3\pm 3.7E-5$
$5.0E-4$ | $7.211E-3\pm 6.4E-5$ | $6.792E-3\pm 5.4E-5$
$7.0E-4$ | $1.0579E-2\pm 9.5E-5$ | $9.622E-3\pm 6.8E-5$
$1.0E-3$ | $1.590E-2\pm 1.5E-4$ | $1.3927E-2\pm 7.3E-5$
$2.0E-3$ | $3.742E-2\pm 2.9E-4$ | $2.855E-2\pm 2.4E-4$
$3.0E-3$ | $6.510E-2\pm 6.3E-4$ | $4.366E-2\pm 3.1E-4$
$4.0E-3$ | $1.0024E-1\pm 5.4E-4$ | $5.898E-2\pm 5.6E-4$
$5.0E-3$ | $1.4327E-1\pm 7.3E-4$ | $7.428E-2\pm 4.4E-4$
$6.0E-3$ | $1.954E-1\pm 1.6E-3$ | –
$7.0E-3$ | $2.560E-1\pm 2.4E-3$ | –
$8.0E-3$ | $3.272E-1\pm 1.7E-3$ | –
$9.0E-3$ | $4.092E-1\pm 3.5E-3$ | –
$1.0E-2$ | $5.014E-1\pm 4.1E-3$ | $1.527E-1\pm 1.1E-3$
Table 6: DMC energies.


# ViralFP: A webserver of viral fusion proteins

Pedro Moreira  Escola de Engenharia da Universidade do MinhoRua da
UniversidadeBragaPortugal4710-057<EMAIL_ADDRESS>,  Sara Pereira
Escola de Engenharia da Universidade do MinhoRua da
UniversidadeBragaPortugal4710-057 ,  Miguel Rocha  Escola de Engenharia da
Universidade do MinhoRua da UniversidadeBragaPortugal4710-057
<EMAIL_ADDRESS>and  Diana Lousa  ITQB NOVA, Universidade Nova de
LisboaAv. da RepúblicaOeirasPortugal2780-157<EMAIL_ADDRESS>
###### Abstract.

Viral fusion proteins are attached to the membrane of enveloped viruses (a
group that includes Coronaviruses, Dengue, HIV and Influenza) and catalyze
fusion between the viral envelope and host membrane, enabling the virus to
insert its genetic material into the host cell. Given the importance of these
biomolecules, this work presents a centralized database containing the most
relevant information on viral fusion proteins, available through a free-to-use
web server accessible accessible by the URL https://viralfp.bio.di.uminho.pt/.
This server also contains several bioinformatic tools, such as Clustal
sequence alignment and Weblogo, and includes a machine learning-based tool
capable of predicting the location of fusion peptides (the component of fusion
proteins that inserts into the host’s cell membrane) within the fusion protein
sequence. Given the crucial role of these proteins in viral infection, their
importance as natural targets of our immune system and their potential as
therapeutic targets, this webserver aims to foster our ability to fight viral
pandemics.

Key words: Fusion Proteins, Fusion Peptides, Web Server, Machine Learning.

##  1\. Introduction

Enveloped viruses, a group of viruses that includes SARS-CoV-2, Dengue,
Influenza and HIV, are characterized by having an exterior lipid envelope that
contains at least one type of glycoprotein on its surface (Marsh, 2005,
Preface I) (Marsh and Helenius, 2006) (Weissenhorn et al., 2007). These
proteins are essential for the entry of the virus into the host cell, because
they enable the binding of the virus to receptors on the host cell and
catalyse the fusion between the viral envelope and the host cell membrane.

In some viruses, one protein can induce both the binding to the host receptor
and the fusion reaction, whereas in other viruses these tasks are executed by
different proteins. Viral fusion proteins (VFP) are key players in a process
called membrane fusion, in which the virus inserts its genetic material into
their host cell, enabling it to replicate (Apellániz et al., 2014) (Harrison,
2015).

During membrane fusion, the two initially separate membrane bilayers, the
viral envelope and the host cellular membrane, are brought into close contact
so to they can merge themselves into one single membrane (Baquero et al.,
2013) (White et al., 2008). According to the currently most accepted model,
after the initial apposition step, membrane fusion passes through a “hemi-
fusion intermediate” phase, in which the outer leaflets of the two membranes
are merged, but the distal ones are still not in contact. This hemi-fused step
forms a stalk-like structure, diminishing the contact area to a minimum,
lowering the work required to overcome hydration force repulsion. Finally,
fusion of the unattached inner leaflets occurs, which starts to form a
hydrophilic pore that will eventually enlarge, connecting the 2 initially
separated intermembrane volumes (Apellániz et al., 2014) (Harrison, 2015)
(Hughson, 1995) (Teissier and Pécheur, 2007).

An enveloped virus has multiple copies of at least one VFP in its pre-fusion
conformation on its outer membrane, which needs to be converted (through a
environmental stimulus trigger) into a membrane embedded pre-hairpin and then
to the final fusion-competent conformation, which is referred as a trimer-of-
hairpins, so that the fusion reaction can happen (Marsh, 2005, pp. 27-28)
(White et al., 2008). VFPs are fundamental for membrane fusion to occur due to
the high kinetic barrier that must be overcome during the process. VFPs are
the catalyst that enables the virus to overcome this energetic barrier
(Harrison, 2015) (Teissier and Pécheur, 2007).

Figure 1. Structure of the app: within the server there are the 3 Docker
containers: the MySQL, with the database dump; the Django (back-end) and the
Angular (front-end). The latter two can be accessed from a client through
those browser links, even though it is the latter that is intended as the one
the client would use.

VFPs can assume multiple possible forms, depending on their corresponding
virus. Currently, the consensus is that VFPs can be grouped into 3 classes: I,
II and III, each with their own characteristics, such their prefusion, pre-
hairpin and post-fusion structures, orientation within virus’ membrane, types
of possible triggering mechanisms and location of their fusion peptide (FPep)
(Apellániz et al., 2014) (Epand, 2003) (Harrison, 2015) (White et al., 2008).

The FPep is a conserved sequence within VFP of the same viral family that can
insert into the host membrane, allowing the fusion of the viral envelope with
the host cell’s membrane. This peptide can promote the reaction of lipid
mixing of membrane vesicles by itself, even if it is detached from the rest of
the metabolite (Apellániz et al., 2014) (Epand, 2003).

This importance of these proteins (and their FPep) within the membrane fusion
process makes them potential therapeutic targets the diseases related to their
respective viruses (Apellániz et al., 2014) (White et al., 2008). One example
of this importance can be found in the defining pandemic of 2020: SARS-CoV-2.
Being an enveloped virus, one of the key therapeutic targets that were
analysed is its VFP, which is called spike (S) protein (Li et al., 2020).
Several of the most promising and effective vaccines that are being developed
(or already administered) for this virus are based on the S protein’s mRNA.
This group of vaccines include the Moderna’s mRNA-1273 (Jackson et al., 2020)
and the Pfizer-BioNTech’s BNT162b2 mRNA (Polack et al., 2020). In both of
them, the mRNA molecule is inserted into a lipid nanoparticle, that stabilizes
the sequence until it can enter our cells. Once inside the cells, they will
produce the protein in its prefusion conformation. These molecules induce an
immune response on the organism, which will produce antibodies to target those
proteins.

Given the known importance of these proteins, a lot of information for each of
these proteins can be found in biological databases, like UniProt (Consortium,
2015), NCBI Protein (Geer et al., 2010) and PDB (Burley et al., 2017).
However, there was not a centralized repository that contained all this
information, forcing researchers to search several databases when they need to
have comprehensive characterization of one or several of these proteins.

To better understand how VFPs and their components are characterized, some
bioinformatic tools can be used, like sequence alignment tools (like BLAST
(Johnson et al., 2008b), Clustal (Sievers and Higgins, 2014) and HMMER (Finn
et al., 2011)), that allow researchers to discover a peptide’s conserved
residues, as well other pieces of information regarding their features. This
allows to obtain valuable insights about the way they operate.

Another way by which these type of proteins are being analysed is through the
use of machine learning (ML) models. Several type of tasks are being developed
using this type of algorithms for VFP, but some of the most promising ones has
been predicting the location of FPeps within VFP sequences (Pereira et al.,
2019) (Sequeira et al., 2019) (Wu et al., 2016) (Wu et al., 2019) and
predicting inhibitors for different VFPs to potentially be used as potential
therapies for some diseases (Xu et al., 2015).

With all that in mind, the main objective of the current work is to develop a
CRUD (Create, Read, Update and Delete) web server that is capable to store,
retrieve, display and analyse the gathered data on enveloped viruses (through
information gathered from NCBI Taxonomy (Federhen, 2012) entries) and their
VFPs, as well to manipulate the server’s information, all through user-
friendly interfaces. This data is stored in a relational database, which
contains relevant data regarding several hundred VFPs. The back-end component
of the application is built using Django (and Django REST API) library, whilst
its front-end is built using Angular, a Typescript framework to build single
page applications. This web service will be freely available and be useful for
the all elements of scientific community who are working with VFPs and FPep.

This application includes selected bioinformatics tools, such as sequence
alignment tools (like BLAST and CLUSTAL), and ML models capable to detect
FPeps within the VFP sequence. Finally, manual curation of the database
information is performed by reviewing the original protein and taxonomy
repositories’ data entries, as well researching through scientific
publications.

##  2\. MATERIALS AND METHODS

The structure of the application developed in this work can be seen in the
Figure 1. There are three distinct components, each with its own Docker
container:

  * •

The MySQL database (to store all the data that the server requires and will
access to display).

  * •

Django back-end, with Django REST framework (to access the database data and
serialize it, besides containing several views that will be called by the
front-end).

  * •

Angular front-end (the application intended for the end user that displays the
database information and contains several tools).

.

In the next topics, a deeper explanation for each of those components will be
made.

###  2.1. Database

Figure 2. Structure of the database

The VFP / FPep data present in the MySQL database was retrieved from UniProt,
NCBI Protein and PDB entries, whilst the taxonomy information was originally
from NCBI Taxonomy. Its relational structure, names of the tables and
attributes, as well of the type of data stored in each attribute can be seen
on the scheme present on figure 2.

The ”main” table of this database is called Protein and contains the
information regarding each VFP’s name, class, type of activation, fusogenic
peptide description (its name, location within the protein and sequence), as
well as having the UniProt and NCBI Protein entries for the protein.

Another important table is the Fusion Peptide, which describes the location of
this peptide within the VFP sequence (attribute “Residues”), sequence and some
attributes regarding annotation method and experimental evidence.

The last key table is the Taxonomy, which has all the relevant taxonomy of the
virus to which the protein belongs (family, genus, species and subspecies /
strain), as well as having the virus’ NCBI Taxonomy entry.

Additional tables within the database include a table with the possible
structures of the VFP (through PDB entries); a table for the virus’ hosts,
which has the common name and NCBI Taxonomy of each virus’ host organism; a
inhibitor / antibody table, that contains information regarding those types of
protein that are associated with a particular VFP, and, finally, tables for
additional bibliographical references.

The data present in the database was manually curated by verifying the
original repositories of the data entries. Finally, the inhibitors and
antibodies with each VFP was filled using Biopython for the Pubmed automated
searches, and pypdb (Gilpin, 2015) for the PDB searches.

###  2.2. Back-end

The main purpose of the back-end is to access the data of the MySQL database.
This component serializes the query data into JSON objects, that can easily be
used by client apps (in this case the front-end), so that later it can be
displayed or saved. Besides that, the back-end allows client apps to add or
modify data, through deserialization methods.

One view for each MySQL database table was built, usually with the same name
of the respective table. Some data was included in the serializers from other
tables accessible through foreign keys whenever that extra information was
relevant. Whenever a query request is made, the returned data will be returned
as response (also in a JSON format). All these pages will allow to receive
requests to insert, update or delete data entries. Besides that, all of those
will contain features to allow to filter and / or search data, besides
containing pagination methods.

Since bioinformatics packages are more commonly written in Python that in
other programming language used in this project (Typescript), those tools are
inserted here, that later will be accessed by the front-end by requests.

Clustal Omega analysis was implemented using a Clustal console. This view will
receive the sequences and the additional parameters as part of the request,
perform the Clustal alignment and provide that result as a text response.

To build Weblogos the packages ${}^{\prime}weblogo^{\prime}$
(https://weblogo.readthedocs.io/) and ${}^{\prime}logomaker^{\prime}$ (Tareen
and Kinney, 2019) are used. If the text output is required, the view sends the
result of ${}^{\prime}weblogo^{\prime}$ only. To obtain the image output, the
weblogo matrix (obtained with the ${}^{\prime}weblogo^{\prime}$ package given
a set of protein sequences) was required to be split into several sub-
matrices, one for each line of the weblogo. In order to use
${}^{\prime}logomaker^{\prime}$, special attention was needed to assure that
the same y-scale was used throughout all weblogo, as well as to make sure that
the last line does not change in size relative to the rest of the lines
(because the package will correct the size of the columns to fill all the
line), something that was achieved by adding columns of zeros in the last sub-
matrix. A base64-string (a string that can be converted into a image object)
will be the output of this view for that option.

The ML view allows to perform the FPep prediction for a given sequence, using
models trained using Propythia package (Sequeira et al., 2021). This is a
Python framework which allows bioinformatic researchers to write code to build
ML models for classification of peptides from features retrieved from their
physicochemical properties. Its key purpose is to facilitate the development
of all the necessary pipelines to build ML models. It achieves this by
providing functions to preprocess data, perform feature selection, run
clustering algorithms, train and optimize ML models, as well as to deal with
protein data and obtain related descriptors from their sequences (Sequeira et
al., 2019) (Sequeira et al., 2021). The models used in this application were
trained from a dataset amiably provided by Sequeira. This file was obtained by
extracting the features (through the functions available in Descriptor module
of Propythia) from a dataset where the non-positive results are composed with
both transmembrane domains sequences and non-FPep and non-transmembrane
domains sequences. Those features went through a feature reduction process
using an SVC model, which, according to Sequeira’s work, is the feature
selection process that tends to produce better ML models for this dataset.
Finally, ML models were trained with the Propythia’s $train\\_best\\_model$
function. Those models are stored in Pickle files. If a certain model is
selected, the file containing that model will be loaded and used in an
instance of the Propythia’s ML class (which uses the already mentioned
dataset). This view also requires selecting a window size and a size for the
gap. After defining those parameters, the $predicted\\_window$ function of the
Propythia is called to perform the predict the likelihood of the subpeptides
of that window size to be FPeps. Those probabilities will be returned in a
JSON object.

A view was added to obtain conservation scores of a peptide against a virus
family’s weblogo. This view receives a request which must contain a protein
sequence, a window size and a taxonomy family, and for each subpeptide of the
size of the window, the function computes the conservation value for the
selected alignment. For each position of the sequence, the view obtains the
highest value for all the conservation features from the peptides that contain
it. A JSON object with those values is returned.

To perform this function, VFP sequences present in the database used in the
web server were used to create those weblogos. After retrieving all the
sequences and the family of each respective virus from the database, the
families that contained at least four sequences in the database were
retrieved. This allowed to select the following 13 families: Arenaviridae,
that have 4 sequences; Coronaviridae (16 sequences); Filoviridae (16);
Flaviviridae (71); Herpesviridae (22); Orthomyxoviridae (170); Paramyxoviridae
(55); Peribunyaviridae (4); Phenuiviridae (4); Pneumoviridae (9); Retroviridae
(164); Rhabdoviridae (6) and Togaviridae (32).

For each of those list of sequences, a Clustal Omega (Sievers et al., 2011)
alignment was performed. The weblogos were created using the ’weblogo’ package
with those alignments.

For each of the weblogos, the ’final’ scores for each position of the PWM
matrix were precalculated, by applying the conservation equation described by
Crooks (Crooks et al., 2004), in which this metric is obtained by the
following equation:

(1) $\ R_{seq}=S_{max}-S_{obs}=log_{2}N-(-\sum_{n=1}^{N}p_{n}*log2(p_{n}))$

where $R_{seq}$ is the conservation value of a given position of the PWM;
$S_{max}$ is the maximal entropy possible in a position, equal to
$log_{2}(N)$, $N$ being the number of possible distinct symbols, which means
that in a protein sequence this value is 20, corresponding to the number of
possible amino acids; and $S_{obs}$ is the observed entropy in a given
position. Because these weblogos are the result of alignments of families with
differing number of protein sequences, that number of sequences was also taken
into the calculations.

So, before obtaining the conservation score for the sequence, the weblogo
matrix was preprocessed by calculating in each position of the PMW the result
of the expression:

(2) $\ PositionScore=\frac{(log_{2}(20)-S_{obs})*bits}{n_{seqs}}$

where the Sobs is the entropy value for that specific row, bits the number of
bits on each position of the PWM and $n_{seqs}$ the number of the sequences
present in the database for the selected family.

After those calculations, there is one Pandas matrix for each family with only
the scores per position. So, after transforming the Pandas dataframe to Numpy
arrays, it is possible to retrieve the local alignment scores for each
sequence (or subsequence) by verifying the scores of the sequence against all
the possible positions in the PWM and save the best of those as a feature.

###  2.3. Front-end

This front-end was built using a framework called Nebular
(https://akveo.github.io/nebular/), which is an open-source Angular User
Interface Library that consists in a pre-built modular, customizable and
configurable application. This enabled to speed up the development process by
already providing a several useful features like tables and card formats that
are user-friendly and visually appealing.

####  2.3.1. Tables’ Pages

The information of the database is found on the “Tables & Data” option of the
menu. This list contains links for the Taxonomy, Fusion Protein and Fusion
Peptides tables. Those pages’ tables (and all the others that can be accessed
through them) are organised in a Tree Grid format: each data entry will
correspond to a table line; upon clicking on one, its related data, which in
this app corresponds to links to related data (the data from other tables
connected by foreign keys), will appear below.

The app’s tools can be grouped together in two categories: sequence alignment
and sequence prediction.

####  2.3.2. Sequence Alignment’ Pages

In the Sequence Alignment tools, a page that allows to send a selected
sequence to the NCBI BLAST was added; this page included methods that allow
the user to select the type of Blast to perform and select the database to run
the BLAST against (Johnson et al., 2008a). Similarly, the HMMER page allows to
select a sequence and the database against which the search is made, and then
send these parameters into the EMBL-EBI portal, showing the results in their
web portal (Potter et al., 2018).

The Weblogo page contains two possible ways to build these graphs. The first
uses the weblogo.threeplusone.com portal (Crooks et al., 2004), in which the
app will ask for at least three sequences, which must be in a FASTA format,
send those to the back-end to perform a Clustal Omega alignment and from that
the web portal builds a weblogo (it can be in a PNG or raw data formats). The
second method will use exclusively the back-end (using the already described
methods): the PNG output will appear in the same page, on which the user can
select the type of colouring provided (from the ones allowed by the
’logomaker’ package), as well as the number of stacks per line; the text one
will appear in a pop-up page.

The Clustal page alignment will use the console from the back-end. The user
can select parameters such as the weight matrix and format of output; upon
clicking the button, one pop-up will open in the browser to display the
alignment in the required format; a second one can appear if the the guide
tree data is requested.

####  2.3.3. Sequence Prediction’ Pages

The app contains two tools relative to sequence prediction: epitope prediction
and FPep prediction.

Regarding the first, the site will access to the IEDB portal (Vita et al.,
2014) in order to obtain an antibody epitope prediction (an epitope is,
according to the IEDB, “molecular targets of adaptive immune responses”); the
app will send a request with the sequence, selected analysis method and window
size to the back-end that will redirect to the IEDB API. The results will be
displayed in a Smart Table format, which is a Nebular feature; both Bepipred
methods’ results will appear into 2 different result tables, the rest of them
just require one (Vita et al., 2015).

Regarding ML prediction of the location of FPeps within the VFP sequence, this
page allows to define a sequence, as well as specifying which model to use,
from the set of available models trained in Propythia: SVM, RF, Gradient
Boosting, KNN, Linear Regression, Gaussian NB and ANN models. Finally, the
user can select the window and gaps sizes. The page can display those
predictions in two types of output: a table which contains every possible
peptide probability of being a FPep; and a graphical interface that, for each
amino acid, displays a color that represents the max probability of all
peptides that contain that particular character of being a FPep, with the
respective color legend (like shown in Figure 3); also, for each amino acid, a
popup box will show up, informing the user about the position within the VFP
sequence and the probability value of belong to a FPep of each character.
These results can be downloaded into a text file.

Figure 3. Example of a Machine Learning prediction of the fusion peptide
location.

Within the FPep prediction page the user can find the conservation tool. This
will display those results in a very similar manner as the one described for
the graphical representation for the FPep prediction. Upon receiving the JSON
object with the output of conservation methods, for each sequence’s position,
the max score of the peptides that include it is provided and from that value
the colour corresponding to that score is attributed to it. The family can be
selected from a list, which upon clicking the button, updates the graphical
output to the desired family.

##  3\. Results

Figure 4. Fusion Protein Page Scheme.

There are two key components of the front-end application that can require an
explanation of their elements: the tables and the tools.

The figure 4 shows the representation of the Fusion Protein page. It will
allow to explain all the tables pages’ different elements. Below there is a
list that describes them, from the left to the right and top to down:

  * •

On the dashed green rectangle on the left, there is the side menu of the
application, where one can access the three main data tables, as well as the
bioinformatic tools present in the server.

  * •

On the top right, in the black box, the user can find the link to the global
search tool; upon entering a search term, a pop-up windows will display
whether the main pages have results for that query; if they do, the page will
present examples of some specific results within those pages.

  * •

The red rectangle contains the search methods also previously described, that
will (partially) determine what is shown in the table; whenever available,
after 3 characters, some auto-complete for the search term suggestions will
appear.

  * •

The blue rectangle allows to download the query result into a CSV file.

  * •

The orange rectangle redirects to the back-end’s administrator page.

  * •

The green rectangle is the table with results from the query request
determined by the search term and the current page; within each line, there
are links to related data (if available).

  * •

The yellow rectangles shows all the necessary tools to allow to use the
selected sequences, chosen in the checkboxes on the yellow dashed rectangle,
or through the “Select All Sequences” button, in the Tools; upon choosing at
least one sequence, the bottom yellow rectangle will show up, allowing the
user can delete some (in the table) or all (in the “Clear All Sequences”) of
the selected sequences and redirect the chosen sequences to the Tools,
accessible by the buttons on the bottom of the page.

  * •

Some columns like the one inside the dashed grey rectangle will contain links
to redirect to the original source of the data; the links can be UniProt, NCBI
Protein, PDB or NCBI Taxonomy data entries.

  * •

The purple box buttons corresponds to the pagination methods that will also
determine the content of the table.

Figure 5. Machine Learning Page.

All the tool pages are similar to the one present at figure 5. Regarding those
tools, all pages contain a text box to insert the sequence(s) input (it can be
auto filled if protein sequences were selected in the Fusion Protein page):
some, mainly the multiple alignment ones, require a minimum of 3 protein
sequences in a FASTA format; the rest just need a protein sequence (the text
boxes contains placeholders with a small description for the type of expected
input). Besides that, the rest of the parameters can be chosen within the
pages, with clear labels for each of them, in the form of buttons, check boxes
or select lists. The results of the tools can appear in new tabs to show text
results, like in the case of Clustal and the text result of Weblogo, or to
access other web domains, like BLAST or HMMER. The results can also appear
within the same page, like in the case of ML prediction and the image result
of Weblogo (built from the back-end’s functions).

##  4\. Contribution

This work provides a significant contribution for the scientific community
working with viruses and their proteins, since it presents a tool that can be
useful to any researcher interested in obtaining several types of information
regarding VFPs.

Currently, the application implements several useful functionalities: it
allows a easy access to the MySQL database through a user-friendly web
interface. All the desired bioinformatics tools are implemented in the
application, through web interfaces. ML models capable of predicting FPeps
within the VFP sequence were trained as described by Sequeira’s work (Sequeira
et al., 2021) and included into the server. The application also contains a
tool to provide conservation scores of a peptide against a virus’ taxonomy
family multiple alignment that can complement the ML prediction. Finally, a
full manual curation of the database information was performed.

##  5\. Future Work

Future tasks within the website will include fixing bugs that may come up and
find and correct unaccounted situations within the different components of
application; include features to facilitate the usage for researchers; improve
the initial loading time of the website and the app’s aesthetic appearance, as
well as removing all the web links present in the data tables that redirect to
empty tables .

Alongside these technical tasks, in order to improve the utility of this
application for the scientific community we will resolve known issues in the
tools, mainly within the Weblogo, Clustal and IEDB’s pages, as well as to add
more options in those pages to increase the possibilities for the output of
those tools; looking into possible solutions to make getting the outputs of
the ML faster and try to implement alternative ways to consider conservation
features so that they can be used as features in ML pipelines.

Regarding the database’s information, it is imperative to add all the relevant
new entries published after the submission of this manuscript, hopefully by
implementing methodologies capable of performing this type of search
automatically and periodically. Finally it is important to redo the automated
searches of known inhibitors / antibodies within Pubmed and PDB in a way that
can take into consideration the virus’ strain.

###### Acknowledgements.

Special thanks to Rúben Rodrigues by being the main help in the development of
the webserver and the deployment into the server. Another special appreciation
to Ana Marta Sequeira, for providing much needed support to understand the
Propythia package and to provide machine learning models and datasets.

## Funding

This work is funded by COMPETE 2020, Portugal 2020 and FCT - Fundação para a
Ciência e a Tecnologia, under the project ”Using computational and
experimental methods to provide a global characterization of viral fusion
peptides”, through the funding program ”02/SAICT/2017 - Projetos de
Investigação Científica e Desenvolvimento Tecnológico (IC&DT)”, with the
reference ”PTDC/CCI-BIO/28200/2017”. This work was also financially supported
by Project LISBOA-01-0145-FEDER-007660 (Microbiologia Molecular, Estrutural e
Celular) funded by FEDER funds through COMPETE2020 — Programa Operacional
Competitividade e Internacionalização (POCI) and by national funds through FCT
— Fundação para a Ciência e a Tecnologia.

## References

  * (1)
  * Apellániz et al. (2014) Beatriz Apellániz, Nerea Huarte, Eneko Largo, and José L. Nieva. 2014.  The three lives of viral fusion peptides.  _Chemistry and Physics of Lipids_ 181 (2014), 40–55.  https://doi.org/10.1016/j.chemphyslip.2014.03.003
  * Baquero et al. (2013) Eduard Baquero, Aurélie A. Albertini, Patrice Vachette, Jean Lepault, Stéphane Bressanelli, and Yves Gaudin. 2013\.  Intermediate conformations during viral fusion glycoprotein structural transition.  _Current Opinion in Virology_ 3, 2 (2013), 143–150.  https://doi.org/10.1016/j.coviro.2013.03.006
  * Burley et al. (2017) Stephen K Burley, Helen M Berman, Gerard J Kleywegt, John L Markley, Haruki Nakamura, and Sameer Velankar. 2017.  Protein Data Bank (PDB): the single global macromolecular structure archive.  _Protein Crystallography_ (2017), 627–641.
  * Consortium (2015) UniProt Consortium. 2015\.  UniProt: a hub for protein information.  _Nucleic acids research_ 43, D1 (2015), D204–D212.
  * Crooks et al. (2004) G Crooks, G Hon, J Chandonia, and S Brenner. 2004\.  WebLogo: A Sequence Logo Generator.  _Genome Research_ 14 (2004), 1188–1190.  https://doi.org/10.1101/gr.849004.1
  * Epand (2003) Richard M. Epand. 2003\.  Fusion peptides and the mechanism of viral fusion.  _Biochimica et Biophysica Acta - Biomembranes_ 1614, 1 (2003), 116–121.  https://doi.org/10.1016/S0005-2736(03)00169-X
  * Federhen (2012) Scott Federhen. 2012\.  The NCBI taxonomy database.  _Nucleic acids research_ 40, D1 (2012), D136–D143.
  * Finn et al. (2011) Robert D Finn, Jody Clements, and Sean R Eddy. 2011\.  HMMER web server: interactive sequence similarity searching.  _Nucleic acids research_ 39, suppl_2 (2011), W29–W37.
  * Geer et al. (2010) Lewis Y Geer, Aron Marchler-Bauer, Renata C Geer, Lianyi Han, Jane He, Siqian He, Chunlei Liu, Wenyao Shi, and Stephen H Bryant. 2010.  The NCBI biosystems database.  _Nucleic acids research_ 38, suppl_1 (2010), D492–D496.
  * Gilpin (2015) William Gilpin. 2015\.  PyPDB: a Python API for the Protein Data Bank.  _Bioinformatics_ 32, 1 (09 2015), 159–160.  https://doi.org/10.1093/bioinformatics/btv543 arXiv:https://academic.oup.com/bioinformatics/article-pdf/32/1/159/6687793/btv543.pdf
  * Harrison (2015) Stephen C. Harrison. 2015\.  Viral membrane fusion.  _Virology_ 479-480 (2015), 498–507.  https://doi.org/10.1016/j.virol.2015.03.043 arXiv:NIHMS150003
  * Hughson (1995) Frederick M. Hughson. 1995\.  Structural characterization of viral fusion proteins.  _Current Biology_ 5, 3 (1995), 265–274.  https://doi.org/10.1016/S0960-9822(95)00057-1
  * Jackson et al. (2020) Lisa A Jackson, Evan J Anderson, Nadine G Rouphael, Paul C Roberts, Mamodikoe Makhene, Rhea N Coler, Michele P McCullough, James D Chappell, Mark R Denison, Laura J Stevens, et al. 2020\.  An mRNA vaccine against SARS-CoV-2—preliminary report.  _New England Journal of Medicine_ (2020).
  * Johnson et al. (2008a) Mark Johnson, Irena Zaretskaya, Yan Raytselis, Yuri Merezhuk, Scott Mcginnis, and Thomas L Madden. 2008a.  NCBI BLAST : a better web interface.  _Nucleic Acids Research_ 36, April (2008), 5–9.  https://doi.org/10.1093/nar/gkn201
  * Johnson et al. (2008b) Mark Johnson, Irena Zaretskaya, Yan Raytselis, Yuri Merezhuk, Scott McGinnis, and Thomas L Madden. 2008b.  NCBI BLAST: a better web interface.  _Nucleic acids research_ 36, suppl_2 (2008), W5–W9.
  * Li et al. (2020) Jing Li, Peng Zhan, and Xinyong Liu. 2020.  Targeting the entry step of SARS-CoV-2: a promising therapeutic approach.  _Signal transduction and targeted therapy_ 5, 1 (2020), 1–2.
  * Marsh (2005) M. Marsh. 2005.  _Current topics in microbiology and immunology_. Vol. 6.  Springer. 207–213 pages.  https://doi.org/10.3109/08820137709055812
  * Marsh and Helenius (2006) Mark Marsh and Ari Helenius. 2006.  Virus entry: Open sesame.  _Cell_ 124, 4 (2006), 729–740.  https://doi.org/10.1016/j.cell.2006.02.007
  * Pereira et al. (2019) Sara Pereira, Diana Lousa, and Isabel Rocha. 2019\.  _Building a database and development of a Machine Learning algorithm to identify and characterize viral Fusion Peptides_.  Master’s thesis. Universidade do Minho.
  * Polack et al. (2020) Fernando P Polack, Stephen J Thomas, Nicholas Kitchin, Judith Absalon, Alejandra Gurtman, Stephen Lockhart, John L Perez, Gonzalo Pérez Marc, Edson D Moreira, Cristiano Zerbini, et al. 2020\.  Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine.  _New England Journal of Medicine_ 383, 27 (2020), 2603–2615.
  * Potter et al. (2018) Simon C Potter, Sean R Eddy, Youngmi Park, Rodrigo Lopez, Robert D Finn, and The Hmmer. 2018\.  HMMER web server : 2018 update.  _Nucleic Acids Research_ 46, June (2018), 200–204.  https://doi.org/10.1093/nar/gky448
  * Sequeira et al. (2019) Ana Sequeira, Miguel Rocha, and Diana Lousa. 2019\.  _Building an automated platform for the classification of peptides/proteins using machine learning_.  Master’s thesis. Universidade do Minho.
  * Sequeira et al. (2021) Ana Marta Sequeira, Diana Lousa, and Miguel Rocha. 2021\.  ProPythia: A Python Automated Platform for the Classification of Proteins Using Machine Learning. In _Practical Applications of Computational Biology & Bioinformatics, 14th International Conference (PACBB 2020)_. Springer, 32–41.
  * Sievers and Higgins (2014) Fabian Sievers and Desmond G Higgins. 2014.  Clustal omega.  _Current protocols in bioinformatics_ 48, 1 (2014), 3–13.
  * Sievers et al. (2011) Fabian Sievers, Andreas Wilm, David Dineen, Toby J. Gibson, Kevin Karplus, Weizhong Li, Rodrigo Lopez, Hamish McWilliam, Michael Remmert, Johannes Söding, Julie D. Thompson, and Desmond G. Higgins. 2011.  Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega.  _Molecular Systems Biology_ 7, 539 (2011).  https://doi.org/10.1038/msb.2011.75
  * Tareen and Kinney (2019) Ammar Tareen and Justin B. Kinney. 2019.  Logomaker: Beautiful sequence logos in python.  _bioRxiv_ (2019).  https://doi.org/10.1101/635029 arXiv:https://www.biorxiv.org/content/early/2019/05/13/635029.full.pdf
  * Teissier and Pécheur (2007) Élodie Teissier and Eve Isabelle Pécheur. 2007\.  Lipids as modulators of membrane fusion mediated by viral fusion proteins.  _European Biophysics Journal_ 36, 8 (2007), 887–899.  https://doi.org/10.1007/s00249-007-0201-z
  * Vita et al. (2015) Randi Vita, James A Overton, Jason A Greenbaum, Julia Ponomarenko, D Clark, Jason R Cantrell, Daniel K Wheeler, Joseph L Gabbard, Deborah Hix, Alessandro Sette, and Bjoern Peters. 2015.  The immune epitope database (IEDB) 3.0.  _Nucleic Acids Research_ 43, October 2014 (2015), 405–412.  https://doi.org/10.1093/nar/gku938
  * Vita et al. (2014) Randi Vita, James A. Overton, Jason A. Greenbaum, Julia Ponomarenko, Jason D. Clark, Jason R. Cantrell, Daniel K. Wheeler, Joseph L. Gabbard, Deborah Hix, Alessandro Sette, and Bjoern Peters. 2014.  The immune epitope database (IEDB) 3.0.  _Nucleic Acids Research_ 43, D1 (10 2014), D405–D412.  https://doi.org/10.1093/nar/gku938 arXiv:https://academic.oup.com/nar/article-pdf/43/D1/D405/7330533/gku938.pdf
  * Weissenhorn et al. (2007) Winfried Weissenhorn, Andreas Hinz, and Yves Gaudin. 2007\.  Virus membrane fusion.  _FEBS Letters_ 581, 11 (2007), 2150–2155.  https://doi.org/10.1016/j.febslet.2007.01.093
  * White et al. (2008) Judith M. White, Sue E. Delos, Matthew Brecher, and Kathryn Schornberg. 2008.  Structures and mechanisms of viral membrane fusion proteins: Multiple variations on a common theme.  _Critical Reviews in Biochemistry and Molecular Biology_ 43, 3 (2008), 189–219.  https://doi.org/10.1080/10409230802058320
  * Wu et al. (2016) Sijia Wu, Jiuqiang Han, Ruiling Liu, Jun Liu, and Hongqiang Lv. 2016.  A computational model for predicting fusion peptide of retroviruses.  _Computational Biology and Chemistry_ 61 (2016), 245–250.  https://doi.org/10.1016/j.compbiolchem.2016.02.013
  * Wu et al. (2019) Sijia Wu, Xiaoming Wu, Jie Tian, Xiaobo Zhou, and Liyu Huang. 2019.  PredictFP2: a new computational model to predict fusion peptide domain in all retroviruses.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ (2019).  https://doi.org/10.1109/tcbb.2019.2898943
  * Xu et al. (2015) Yongtao Xu, Shui Yu, Jian Wei Zou, Guixiang Hu, Noorsaadah A.B.D. Rahman, Rozana Binti Othman, Xia Tao, and Meilan Huang. 2015\.  Identification of peptide inhibitors of enveloped viruses using support vector machine.  _PLoS ONE_ 10, 12 (2015), 1–15.  https://doi.org/10.1371/journal.pone.0144171