# Optimized Quantum Drude Oscillators for Atomic and Molecular Response
Properties

Szabolcs Góger  Almaz Khabibrakhmanov  Ornella Vaccarelli  Dmitry V. Fedorov
Alexandre Tkatchenko<EMAIL_ADDRESS>Department of Physics and
Materials Science, University of Luxembourg, L-1511 Luxembourg City,
Luxembourg

###### Abstract

The quantum Drude oscillator (QDO) is an efficient yet accurate coarse-grained
approach that has been widely used to model electronic and optical response
properties of atoms and molecules, as well as polarization and dispersion
interactions between them. Three effective parameters (frequency, mass,
charge) fully characterize the QDO Hamiltonian and are adjusted to reproduce
response properties. However, the soaring success of _coupled_ QDOs for many-
atom systems remains fundamentally unexplained and the optimal mapping between
atoms/molecules and oscillators has not been established. Here, we present an
optimized parametrization (OQDO) where the parameters are fixed by using only
dipolar properties. For the periodic table of elements as well as small
molecules, our OQDO model accurately reproduces atomic (spatial) polarization
potentials and multipolar dispersion coefficients, elucidating the high
promise of the presented model in the development of next-generation quantum-
mechanical force fields for (bio)molecular simulations.

The development of predictive model Hamiltonians that can describe various
properties of realistic molecules and materials is a cornerstone of modern
physics Cipcigan  _et al._ (2019) and chemistry Akimov and Prezhdo (2015). The
quantum Drude oscillator (QDO) is arguably the most powerful Hamiltonian for
accurate and efficient modeling of atomic and molecular response Hermann  _et
al._ (2017a); Cipcigan  _et al._ (2019); Wang and Jordan (2001); Sommerfeld
and Jordan (2005); Jones  _et al._ (2013a); Sadhukhan and Manby (2016); Lemkul
_et al._ (2016); Schröder and Steinhauser (2010); Anisimov  _et al._ (2005).
Within the coarse-grained QDO model, the response of valence electrons is
described _via_ a quasi-particle _drudon_ with a negative charge $-q$ and mass
$\mu$, harmonically bound to a positively-charged pseudo-nucleus of charge $q$
with a characteristic frequency $\omega$. The many-body extension of the QDO
model (the _coupled_ QDO model) has been widely employed to study both
molecules and materials, including their electronic Thole (1981); Hermann  _et
al._ (2017b) and optical Ambrosetti  _et al._ (2022) properties, polarization
Whitfield and Martyna (2006); Voora  _et al._ (2013), dispersion Bade (1957);
Bade and Kirkwood (1957); Whitfield and Martyna (2006); Tkatchenko  _et al._
(2012); Reilly and Tkatchenko (2015); Sadhukhan and Manby (2016); Sadhukhan
and Tkatchenko (2017); Ambrosetti  _et al._ (2014); DiStasio Jr.  _et al._
(2014); Stöhr  _et al._ (2021); Karimpour  _et al._ (2022a), and exchange
Vaccarelli  _et al._ (2021); Fedorov  _et al._ (2018); Tkatchenko  _et al._
(2021) interactions, as well as a wealth of non-additive field effects in
quantum mechanics Kleshchonok and Tkatchenko (2018); Stöhr  _et al._ (2021)
and quantum electrodynamics Karimpour  _et al._ (2022b, a). Coupled QDOs are
also extensively used in the development of van der Waals (vdW) density
functionals Tkatchenko and Scheffler (2009); Tkatchenko  _et al._ (2012);
Ambrosetti and Silvestrelli (2012), quantum mechanical Jones  _et al._
(2013a); Cipcigan  _et al._ (2019) and polarizable force fields Harder  _et
al._ (2006); Lopes  _et al._ (2013); Adluri  _et al._ (2015); Sokhan  _et al._
(2015); Poier  _et al._ (2022a) as well as recent machine learning force
fields Muhli  _et al._ (2021); Poier  _et al._ (2022b). Despite such a wide
applicability of the coupled QDO model, its success in describing real atoms
remains fundamentally unexplained and the optimal mapping between atoms and
oscillators has not been established. In this Letter, we develop an optimized
parametrization (OQDO) where the parameters are fixed by using only the well-
known atomic dipolar properties. Remarkably, OQDO reproduces spatial atomic
polarization potentials and atomic multipolar dispersion coefficients. Our
OQDO model for atoms and small molecules also paves the way to develop next-
generation quantum-mechanical force fields for (bio)molecular simulations.

Figure 1: Polarization potential curves $V_{\mathrm{pol}}(\mathbf{r})$
calculated with DFT-PBE0 and various QDO parametrizations for (a) hydrogen,
(b) carbon (no JQDO values are available), (c) sodium, and (d) argon atoms.
The FQDO and JQDO parametrization schemes are described by Eqs. (1) and (2),
respectively. OQDO(A) and OQDO(B) correspond to the two solutions of the
transcendental equation given by Eq. (3). In all cases the direction along the
applied field was chosen for the plots. The reference values for dipole
polarizability $\alpha_{1}$ are shown for each element. The numerical values
of the normalized root-mean-squared error ($\Delta$) are displayed for the
three QDO flavors. For hydrogen and carbon atoms, the unperturbed
electrostatic potentials (ESP) $V_{\mathrm{el}}(\mathbf{r})$ are shown as
insets, indicating that a QDO captures the response of atomic electron
density, but not the static potential itself.

The three parameters {$q,\mu,\omega$} fully define the QDO, and three atomic
response properties could be chosen to fix them, meaning that the choice of
QDO parameters is not unique. In addition, all QDO response properties –
multipolar polarizabilities and dispersion coefficients – are uniquely fixed
by the three parameters _via_ closed-form relations Jones  _et al._ (2013a).
The static dipole polarizability of a QDO, $\alpha_{1}=q^{2}/\mu\omega^{2}$,
conveniently combines all three parameters, and it is natural to set this
expression to the reference atomic $\alpha_{1}$. The QDO expression for the
dipole-dipole dispersion coefficient
$C_{6}=\frac{3}{4}\hbar\omega\alpha_{1}^{2}$ is identical to the London
formula and allows to fix $\omega$ if the reference atomic values of $C_{6}$
and $\alpha_{1}$ are given. Since $\alpha_{1}$ and $C_{6}$ are accurately
known for all elements in the periodic table Gobre (2016); Gould and Bučko
(2016); Schwerdtfeger and Nagle (2019), they form a baseline for the QDO
parametrization. However, one more condition is required to obtain
{$q,\mu,\omega$}, for which different constraints can be imposed. A reasonable
idea is to fix $q=1$ a.u., since a QDO should reproduce the response of
electrons. This results in the fixed-charge QDO (FQDO)

$q=1\;,~{}~{}\omega={4C_{6}}/{3\hbar\alpha_{1}^{2}}\;,~{}~{}\mu={9\hbar^{2}\alpha_{1}^{3}}/{16C_{6}^{2}}\quad.$
(1)

However, fixing $q$ and using QDO recursion relations for high-order response
usually yields large errors in the multipolar response properties (see Fig. 2
and Refs. Johnson and Becke (2006); Grimme  _et al._ (2010)). A more rigorous
approach was suggested by Jones _et al._ Jones  _et al._ (2013a) by employing
the dipole-quadrupole dispersion coefficient $C_{8}$. The mapping
$\\{\alpha_{1},C_{6},C_{8}\\}\rightarrow\\{q,\mu,\omega\\}$ yields the Jones
QDO (JQDO) parametrization scheme

$q=\sqrt{\mu\omega^{2}\alpha_{1}}\;,~{}~{}\omega={4C_{6}}/{3\hbar\alpha_{1}^{2}}\;,~{}~{}\mu={5\hbar
C_{6}}/{\omega C_{8}}\ .$ (2)

The JQDO approach improves the multipolar response over the FQDO model, while
simulations using the coupled JQDO model captured many remarkable properties
of bulk water and its surface Jones  _et al._ (2013b); Sokhan  _et al._
(2015). However, the $C_{8}$ dispersion coefficient is not directly
measurable, and accurate _ab initio_ calculations of quadrupole ($\alpha_{2}$)
and octupole ($\alpha_{3}$) polarizabilities and $C_{8}$ – $C_{10}$ dispersion
coefficients are currently technically feasible only for closed-shell species
(noble-gas atoms and small molecules) or alkali and alkaline-earth atoms with
$s$ valence shells Porsev and Derevianko (2003, 2006); Jiang  _et al._ (2015);
Tao and Rappe (2016). For other open-shell atoms (containing $p$, $d$, or $f$
valence shells), convergence of _ab initio_ response calculations becomes a
technical hurdle Rob . Thus, using higher-order atomic response properties
does not lead to a parametrization that would be universally applicable across
the periodic table as well as for small molecules.

Here, we introduce an optimized QDO parametrization (OQDO), where we
effectively map dipolar atomic quantities {$\alpha_{1},C_{6}$} to the
oscillator parameters. The third parameter is fixed by using the force balance
equation for vdW-bonded dimers derived recently Fedorov  _et al._ (2018);
Vaccarelli  _et al._ (2021); Tkatchenko  _et al._ (2021). Two equations for
$q$ and $\omega$ follow the JQDO scheme, whereas the third one is replaced
with a transcendental equation for a product $\mu\omega$ to be solved
numerically (_vide infra_)

$\mu=\frac{5\,\hbar\,{\rm C}_{6}}{\omega\,{\rm
C}_{8}}\rightarrow\exp(\frac{2\mu\omega R_{\rm
vdW}^{2}}{\hbar})=\frac{2^{7}\cdot(\alpha_{\rm
fsc}^{-1/3}a_{0})^{4}}{(3\hbar/\mu\omega)^{2}}\;,$ (3)

where $\alpha_{\rm fsc}=e^{2}/4\pi\varepsilon_{0}\hbar c$ is the fine-
structure constant, and the vdW radius ($R_{\rm vdW}$) is calculated _via_ the
universal formula connecting it with the dipole polarizability

$\alpha_{1}(R_{\rm vdW})={(4\pi\varepsilon_{0})\,R_{\rm
vdW}^{7}}/{(\alpha_{\rm fsc}^{-1/3}a_{0})^{4}}\ ,$ (4)

as obtained from the balance between exchange repulsion and vdW dispersion
attraction forces at the equilibrium distance in homonuclear atomic dimers
Fedorov  _et al._ (2018); Vaccarelli  _et al._ (2021).

Figure 2: (a) Periodic variations of the QDO length
$\sigma=\sqrt{\hbar/2\mu\omega}$ with the atomic number for the three
different parametrizations, as compared to the atomic vdW radii
($R_{\mathrm{vdW}}$) which are evaluated _via_ Eq. (4) using the reference
atomic polarizabilities Gobre (2016). (b) Correlation between
$R_{\mathrm{vdW}}$ and $\sigma$ (within the three QDO parametrizations) for
102 elements in the periodic table. The schemes OQDO(A) and OQDO(B) correspond
to two solutions of Eq. (10), as illustrated using the example of Ar in the
inset.

The predictive power of the QDO Hamiltonian arises from its ability to
effectively model integrated atomic response. However, when modeling a
molecule or a solid, coupled QDOs must capture interactions between atoms.
Considering two atoms $i$ and $j$ and using interatomic perturbation theory
Stone (2013); Jeziorski  _et al._ (1994), the interaction energy can be
written as the integrated product of the atomic electron density of $i$ with
the electric potential of atom $j$ Kaplan (2006); Stone (2013)

$E_{\mathrm{int}}=\int\rho_{\rm i}(\mathbf{r})V_{\rm
j}(\mathbf{r})\differential^{3}\mathbf{r}\ .$ (5)

This formula is valid for all dominant noncovalent interactions:
electrostatics, induction/polarization, exchange-repulsion, and dispersion.
The validity is evident for the former two cases Kaplan (2006); Stone (2013),
and it was shown that dispersion von Lilienfeld  _et al._ (2004) and exchange
Van Vleet  _et al._ (2016) interactions can also be _effectively_ represented
using the form of Eq. (5).

Response properties are given by variations of $E_{\rm int}$ as

$\delta E_{\mathrm{int}}=\int\left(\delta\rho_{\rm i}(\mathbf{r})V_{\rm
j}(\mathbf{r})+\rho_{\rm i}(\mathbf{r})\delta V_{\rm
j}(\mathbf{r})\right)\differential^{3}\mathbf{r}\;.$ (6)

For an external electric field $\mathbf{E}$, which can also model the effect
of environment,
$\delta\rho(\mathbf{r})=\rho_{E}(\mathbf{r})-\rho(\mathbf{r})$, where
$\rho_{E}(\mathbf{r})$ is the electron density under the external field. Then,
the dominant contribution to $\delta V_{\rm j}(\mathbf{r})$ is given by the
polarization potential Com

$V_{\mathrm{pol}}(\mathbf{r})=-\frac{1}{4\pi\varepsilon_{0}}\int\frac{\rho_{E}(\mathbf{r^{\prime}})-\rho(\mathbf{r^{\prime}})}{|\mathbf{r-r^{\prime}}|}\differential^{3}\mathbf{r^{\prime}}\;,$
(7)

which describes the change in the electrostatic potential of the system due to
the polarization of its charge density by the presence of another moiety (an
electric field in this case). For the QDO, the integral in Eq. (7) can be
evaluated analytically, yielding the following expression

$V_{\mathrm{pol}}^{{}_{\mathrm{QDO}}}(\mathbf{r})=\frac{-q}{4\pi\varepsilon_{0}}\left(\frac{\text{erf}(\tilde{r}/\sigma\sqrt{2})}{\tilde{r}}-\frac{\text{erf}(r/\sigma\sqrt{2})}{r}\right)\;,$
(8)

where $\tilde{\mathbf{r}}=\mathbf{r}+\alpha_{1}\mathbf{E}/q$ is the field-
induced oscillator coordinate and $\sigma=\sqrt{\hbar/2\mu\omega}$ is the QDO
spread Szabó  _et al._ (2022). In the Supplemental Material (SM) SM , we
present $V_{\mathrm{pol}}^{{}_{\mathrm{QDO}}}(\mathbf{r})$ in comparison to
$V_{\mathrm{pol}}(\mathbf{r})$ calculated for 21 atoms (between H – Ca and Kr)
within hybrid density-functional theory DFT-PBE0 Adamo and Barone (1999);
Perdew  _et al._ (1996); Epifanovsky  _et al._ (2021); Lu and Chen (2012);
Zhang and Lu (2021) shown to yield a highly accurate description of electronic
response Hait and Head-Gordon (2018) comparable to coupled-cluster
calculations SM . Here, we remark that the strength of the electric field was
chosen individually for each element depending on its reference static dipole
polarizability Gobre (2016); Schwerdtfeger and Nagle (2019) so that the field-
induced dipole moment is set as $\mathbf{d}=\alpha_{1}\mathbf{E}=0.01$ a.u.,
for all atoms CCC .

Before comparing $V_{\mathrm{pol}}({\bf r})$ for real atoms with different QDO
flavors, it is instructive to consider which atomic properties can be
faithfully captured by a QDO. Firstly, the QDO does not aim to describe static
properties of the atomic electron density, but rather its response under
applied static and fluctuating fields, as demonstrated also by the insets in
Fig. 1a,b. The electrostatic potential (ESP) of a QDO is given by $V_{\rm
el}^{{}_{\rm QDO}}=-q\cdot{\rm erf}(r/\sigma\sqrt{2})/r$, so the charge $q$
determines its magnitude. This explains why $V_{\rm el}^{{}_{\rm FQDO}}$
yields a good agreement with $V_{\rm el}^{\rm DFT}$ for hydrogen. However, the
QDO model does not describe $V_{\rm el}$ for many-electron atoms because
$q\sim 1$ a.u., while the ESP of atoms scales non-linearly with $Z$ (see the
example of carbon in the inset of Fig. 1b). Secondly, the harmonic response
captured by a QDO model should be sufficient to accurately describe integrated
electronic displacements induced by weak fields. However, it is much less
clear how well different QDO parameterizations perform for distributed
polarization potentials described by Eq. (7) for many-electron systems, given
the analytical form of $V_{\mathrm{pol}}^{{}_{\mathrm{QDO}}}(\mathbf{r})$ in
Eq. (8). To answer this question, in Fig. 1 we compare $V_{\mathrm{pol}}$ of
real atoms and $V_{\mathrm{pol}}^{{}_{\mathrm{QDO}}}(\mathbf{r})$ employing
the three QDO models discussed above. We used the accurate _ab initio_
reference data on $\alpha_{1}$ and $C_{6}$ Gobre (2016); Schwerdtfeger and
Nagle (2019); Derevianko  _et al._ (2010) to parametrize FQDO and OQDO. When
available, we also used the analogous data on $C_{8}$ Porsev and Derevianko
(2003, 2006); Jiang  _et al._ (2015) to parametrize the JQDO model. We observe
that the OQDO model is able to reproduce the full range of the polarization
potential of real atoms with a reasonable accuracy, showing a significantly
better agreement with the DFT-PBE0 results than FQDO and JQDO. To quantify
this, for each atom we calculated the root-mean-squared-error (RMSE) of the
three QDO curves with respect to the PBE0 reference curves and normalized the
RMSE using the equilibrium depth of the PBE0 curve. The OQDO flavor has an
error of 8.9% when averaged over 21 atoms, whereas JQDO and FQDO yield average
errors of 13.2% and 15.4%, respectively. We also emphasize that the
predictions of the OQDO model remain accurate for many-electron atoms such as
noble gases and alkali metals. It is especially reassuring that the OQDO model
reproduces the nonlinear $V_{\rm pol}(\mathbf{r})$ curves obtained from DFT
calculations without any adjustments. In fact, the OQDO performance is
sensitive to variations in the QDO parameters (solutions A or B in Fig. 1), so
the satisfactory agreement shows that the chosen OQDO(A) model accurately
describes real atoms. The significant differences between the predictions of
various parametrizations for $V_{\mathrm{pol}}(\mathbf{r})$ underline the
importance of the optimal mapping between atomic response properties and QDO
parameters.

Figure 3: Multipolar dispersion coefficients $C_{8}$ and $C_{10}$ as predicted
by FQDO ($q=1$) and OQDO models. Relative error
${\rm{RE}}=(C_{j}-C_{j}^{\rm{ref}})/C_{j}^{\rm{ref}}$ with respect to _ab
initio_ reference data Porsev and Derevianko (2003, 2006); Jiang  _et al._
(2015); Tao and Rappe (2016) is plotted. For the two models, numerical values
of mean absolute relative errors (MARE) are evaluated separately for atoms and
molecules. In the case of $\mathrm{O_{2}}$ and $\mathrm{CO_{2}}$, no reliable
_ab initio_ reference data for $C_{10}$ could be found.

We discuss now the technical aspects of deriving the two solutions of the OQDO
model (see OQDO(A) and OQDO(B) in Fig. 1), and their connection to real atoms.
The starting point is Eq. (4) that connects the atomic vdW radius and its
dipole polarizability. Within the QDO model, Eq. (4) can be written as follows
Fedorov  _et al._ (2018):

$\alpha_{1}(\mu\omega,R_{\rm
vdW})=\frac{2^{7}\cdot(4\pi\varepsilon_{0})\,R_{\rm
vdW}^{7}}{(3\hbar/\mu\omega)^{2}\,\exp(2\mu\omega R_{\rm vdW}^{2}/\hbar)}\ .$
(9)

The OQDO parametrization imposes that the product $\mu\omega$ in Eq. (9)
delivers the same $R_{\rm vdW}$ as from Eq. (4), for
$\alpha_{1}(\mu\omega,R_{\rm vdW})=\alpha_{1}(R_{\rm vdW})$. For simplicity,
we rewrite Eq. (9) in terms of the dimensionless variable $x$ as

$x=a\,e^{b\,x}\ \ ,\ \ \
x=\mu\,\omega\,a_{0}^{2}/\hbar={a_{0}^{2}}/{2\sigma^{2}}\ \ ,$ (10)

with the dimensionless coefficients $a$ and $b$ given by

$a=\frac{3\,\alpha_{\rm fsc}^{\nicefrac{{2}}{{3}}}}{8\sqrt{2}}\ ,\ \
b=\frac{R_{\rm
vdW}^{2}}{a_{0}^{2}}=\frac{(\alpha_{1}/4\pi\varepsilon_{0})^{\nicefrac{{2}}{{7}}}}{\alpha_{\rm
fsc}^{\nicefrac{{8}}{{21}}}a_{0}^{\nicefrac{{6}}{{7}}}}\ ,$ (11)

where we used Eq. (4) to express $R_{\rm vdW}$ in terms of $\alpha_{1}\,$. For
all elements in the periodic table, we found that Eq. (10) has two solutions,
A and B. This situation is illustrated by the inset of Fig. 2b for the case of
the argon atom. To rationalize this, it is instructive to consider when Eq.
(10) has only one solution, which happens when the polarizability of an atom
is equal to the critical value:

${\alpha_{1}^{\rm(c)}}/{4\pi\varepsilon_{0}}=\left({8\sqrt{2}}/{3e}\right)^{{}_{\nicefrac{{7}}{{2}}}}\\!\\!\alpha_{\rm
fsc}^{-1}\,a_{0}^{3}\,\approx\,431~{}{\rm a.u.}\,,$ (12)

which is greater than the largest known atomic dipole polarizability
($\alpha_{1}\approx 400$ a.u.) of Cs Schwerdtfeger and Nagle (2019). It is
worth mentioning that the existence of the two solutions extends beyond the
employed QDO model. To demonstrate this point, we obtained an analogous result
by using the Tang-Toennies potential Tang and Toennies (2003) with the
repulsive interaction treated by the Born-Mayer form (see SM SM ).

Since the OQDO frequency is fixed by the second condition of Eqs. (1) and (2),
the solutions A and B for the product $\mu\omega$ differ both in mass and
charge, yielding quite different results. First,
$V_{\mathrm{pol}}^{{}_{\mathrm{QDO}}}(\mathbf{r})$ constructed from solution B
do not resemble DFT potentials, while A is in good agreement with them (Fig.
1). Second, the overlap integral
$S=\exp\left({-\frac{\mu\omega}{2\hbar}R_{\mathrm{eq}}^{2}}\right)$ between
two QDOs at their equilibrium distance $R_{\mathrm{eq}}=2R_{\mathrm{vdW}}$ is
significantly larger for solution B, which violates the initial assumption
used to derive Eq. (9) that $S$ is small at $R_{\mathrm{eq}}$ Fedorov  _et
al._ (2018). Third, the QDO length $\sigma$ constructed from solution A
follows the same periodic trend as the atomic vdW radii, whereas solution B
does not seem to correlate well (see Fig. 2a). Therefore, throughout this work
we refer to solution A as the optimized and most promising parametrization.
For 102 atoms, the full set of QDO parameters corresponding to both solutions
A and B is given in the SM SM  together with the reference values of
{$\alpha_{1},C_{6}$}. Another noteworthy property of the OQDO model (see Fig.
2b) is a quasi-linear correlation between the QDO length (model quantity) and
the atomic vdW radius (physical observable). In fact, these quantities should
be connected _via_ the dipole polarizability Szabó  _et al._ (2022); Fedorov
_et al._ (2018). This property is not captured well by either the FQDO or the
JQDO models.

For practical calculations of the vdW energy and constructing predictive force
fields, the multipolar contributions associated with the $C_{8}$ and $C_{10}$
coefficients can become relevant Johnson and Becke (2006); Grimme  _et al._
(2010); Tao and Rappe (2016). The available reference data for higher-order
molecular dispersion coefficients have significant uncertainties. Our careful
examination of the literature reporting the reference values of $C_{8}$ and
$C_{10}$ (see Porsev and Derevianko (2003, 2006); Jiang  _et al._ (2015); Tao
and Rappe (2016) and references therein) identifies uncertainties of up to 20%
for the reference $C_{8}$ and $C_{10}$ values. Within the QDO formalism, it is
straightforward to evaluate these coefficients using closed-form expressions
derived by Jones _et al._ Jones  _et al._ (2013a). In Fig. 3, we present the
predictions of $C_{8}$ and $C_{10}$ by FQDO and OQDO models as compared to
accurate reference values compiled from the literature Porsev and Derevianko
(2003, 2006); Jiang  _et al._ (2015); Tao and Rappe (2016) for a set of 16
atoms (including alkali and alkaline-earth metals and noble gases) and 12
small molecules. We do not include JQDO results in this comparison, since it
explicitly uses $C_{8}$ as an input parameter. Overall, our results show that
the OQDO parametrization improves the dispersion coefficients compared to the
FQDO one, reducing the MARE from 31% to 25% for $C_{8}$ and from 68% to 33%
for $C_{10}$ when averaged over all 28 (26 for $C_{10}$) systems. The OQDO
model consistently surpasses the FQDO in accuracy for all systems, except for
alkaline-earth metals where the FQDO gives more accurate results. Moreover,
Fig. 3 shows that the deviations of OQDO dispersion coefficients from the
reference values are consistent in terms of their sign and magnitude. Namely,
for majority of systems FQDO underestimates $C_{8}$ and $C_{10}$, but roughly
for one third of them the dispersion coefficients are overestimated. The
maximal errors of FQDO are observed for Xe both in case of $C_{8}$ (66%) and
$C_{10}$ (190%). In contrast, OQDO consistently underestimates both dispersion
coefficients for all systems, except for $C_{8}$ of Li as well as $C_{10}$ of
Li and Cs, where there is overestimation. The maximum errors of OQDO are
observed for ${\rm CO_{2}}$ (46%) in case of $C_{8}$ and CO (58%) in case of
$C_{10}$, which are significantly smaller than maximum errors of FQDO. The
consistency of OQDO errors allows a straightforward rescaling of dispersion
coefficients: with our best rescaling factors, 1.3 for $C_{8}$ and 1.5 for
$C_{10}$, one can decrease the MARE of OQDO to 15% and 22%, respectively,
which is not far from the uncertainty of the reference molecular $C_{8}$ and
$C_{10}$ values. A similar analysis is carried out for multipolar
polarizabilities $\alpha_{2}$ and $\alpha_{3}$ in the SM SM . In addition, the
OQDO model reproduces spatially-distributed polarization potentials rather
well for most elements in the periodic table. This highlights the promise of
OQDO for constructing _coupled_ OQDO models for complex molecules and
materials.

Our robust OQDO parametrization for all elements in the periodic table
substantially advances the ability to model a wide range of response
properties of molecules and materials based on the coupled QDO framework, also
paving the way to develop next-generation quantum-mechanical force fields for
(bio)molecular simulations.

S.G., A.K., and O.V. contributed equally to this work. We acknowledge
financial support from the European Research Council (ERC Consolidator Grant
“BeStMo”), the Luxembourg National Research Fund _via_ FNR CORE Jr project
“PINTA(C17/MS/11686718)”, “DRIVEN (PRIDE17/12252781)”, and “ACTIVE
(PRIDE19/14063202)”.

## References

  * Cipcigan  _et al._ (2019) F. S. Cipcigan, J. Crain, V. P. Sokhan, and G. J. Martyna, Rev. Mod. Phys. 91, 025003 (2019).
  * Akimov and Prezhdo (2015) A. V. Akimov and O. V. Prezhdo, Chem. Rev. 115, 5797 (2015).
  * Hermann  _et al._ (2017a) J. Hermann, R. A. DiStasio Jr., and A. Tkatchenko, Chem. Rev. 117, 4714 (2017a).
  * Wang and Jordan (2001) F. Wang and K. D. Jordan, J. Chem. Phys. 114, 10717 (2001).
  * Sommerfeld and Jordan (2005) T. Sommerfeld and K. D. Jordan, J. Phys. Chem. A 109, 11531 (2005).
  * Jones  _et al._ (2013a) A. P. Jones, J. Crain, V. P. Sokhan, T. W. Whitfield, and G. J. Martyna, Phys. Rev. B 87, 144103 (2013a).
  * Sadhukhan and Manby (2016) M. Sadhukhan and F. R. Manby, Phys. Rev. B 94, 115106 (2016).
  * Lemkul  _et al._ (2016) J. A. Lemkul, J. Huang, B. Roux, and A. D. MacKerell, Chem. Rev. 116, 4983 (2016).
  * Schröder and Steinhauser (2010) C. Schröder and O. Steinhauser, J. Chem. Phys. 133, 154511 (2010).
  * Anisimov  _et al._ (2005) V. M. Anisimov, G. Lamoureux, I. V. Vorobyov, N. Huang, B. Roux, and A. D. MacKerell, J. Chem. Theory Comput. 1, 153 (2005).
  * Thole (1981) B. Thole, Chem. Phys. 59, 341 (1981).
  * Hermann  _et al._ (2017b) J. Hermann, D. Alfè, and A. Tkatchenko, Nat. Commun. 8, 14052 (2017b).
  * Ambrosetti  _et al._ (2022) A. Ambrosetti, P. Umari, P. L. Silvestrelli, J. Elliott, and A. Tkatchenko, Nat. Commun. 13, 813 (2022).
  * Whitfield and Martyna (2006) T. W. Whitfield and G. J. Martyna, Chem. Phys. Lett. 424, 409 (2006).
  * Voora  _et al._ (2013) V. K. Voora, J. Ding, T. Sommerfeld, and K. D. Jordan, J. Phys. Chem. B 117, 4365 (2013).
  * Bade (1957) W. L. Bade, J. Chem. Phys. 27, 1280 (1957).
  * Bade and Kirkwood (1957) W. L. Bade and J. G. Kirkwood, J. Chem. Phys. 27, 1284 (1957).
  * Tkatchenko  _et al._ (2012) A. Tkatchenko, R. A. DiStasio Jr., R. Car, and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  * Reilly and Tkatchenko (2015) A. M. Reilly and A. Tkatchenko, Chem. Sci. 6, 3289 (2015).
  * Sadhukhan and Tkatchenko (2017) M. Sadhukhan and A. Tkatchenko, Phys. Rev. Lett. 118, 210402 (2017).
  * Ambrosetti  _et al._ (2014) A. Ambrosetti, D. Alfè, R. A. DiStasio Jr., and A. Tkatchenko, J. Phys. Chem. Lett. 5, 849 (2014).
  * DiStasio Jr.  _et al._ (2014) R. A. DiStasio Jr., V. V. Gobre, and A. Tkatchenko, J. Phys.: Condens. Matter 26, 213202 (2014).
  * Stöhr  _et al._ (2021) M. Stöhr, M. Sadhukhan, Y. S. Al-Hamdani, J. Hermann, and A. Tkatchenko, Nat. Commun. 12, 137 (2021).
  * Karimpour  _et al._ (2022a) M. R. Karimpour, D. V. Fedorov, and A. Tkatchenko, Phys. Rev. Research 4, 013011 (2022a).
  * Vaccarelli  _et al._ (2021) O. Vaccarelli, D. V. Fedorov, M. Stöhr, and A. Tkatchenko, Phys. Rev. Research 3, 033181 (2021).
  * Fedorov  _et al._ (2018) D. V. Fedorov, M. Sadhukhan, M. Stöhr, and A. Tkatchenko, Phys. Rev. Lett. 121, 183401 (2018).
  * Tkatchenko  _et al._ (2021) A. Tkatchenko, D. V. Fedorov, and M. Gori, J. Phys. Chem. Lett. 12, 9488 (2021).
  * Kleshchonok and Tkatchenko (2018) A. Kleshchonok and A. Tkatchenko, Nat. Commun. 9, 3017 (2018).
  * Karimpour  _et al._ (2022b) M. R. Karimpour, D. V. Fedorov, and A. Tkatchenko, J. Phys. Chem. Lett. 13, 2197 (2022b).
  * Tkatchenko and Scheffler (2009) A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  * Ambrosetti and Silvestrelli (2012) A. Ambrosetti and P. L. Silvestrelli, Phys. Rev. B 85, 073101 (2012).
  * Harder  _et al._ (2006) E. Harder, V. M. Anisimov, I. V. Vorobyov, P. E. M. Lopes, S. Y. Noskov, A. D. MacKerell, and B. Roux, J. Chem. Theory Comput. 2, 1587 (2006).
  * Lopes  _et al._ (2013) P. E. M. Lopes, J. Huang, J. Shim, Y. Luo, H. Li, B. Roux, and A. D. MacKerell, J. Chem. Theory Comput. 9, 5430 (2013).
  * Adluri  _et al._ (2015) A. N. S. Adluri, J. N. Murphy, T. Tozer, and C. N. Rowley, J. Phys. Chem. B 119, 13422 (2015).
  * Sokhan  _et al._ (2015) V. P. Sokhan, A. P. Jones, F. S. Cipcigan, J. Crain, and G. J. Martyna, Proc. Natl. Acad. Sci. 112, 6341 (2015).
  * Poier  _et al._ (2022a) P. P. Poier, L. Lagardère, and J.-P. Piquemal, J. Chem. Theory Comput. 18, 1633 (2022a).
  * Muhli  _et al._ (2021) H. Muhli, X. Chen, A. P. Bartók, P. Hernández-León, G. Csányi, T. Ala-Nissila, and M. A. Caro, Phys. Rev. B 104, 054106 (2021).
  * Poier  _et al._ (2022b) P. P. Poier, T. Jaffrelot Inizan, O. Adjoua, L. Lagardère, and J.-P. Piquemal, J. Phys. Chem. Lett. 13, 4381 (2022b).
  * Gobre (2016) V. V. Gobre,  _Efficient modelling of linear electronic polarization in materials using atomic response functions_ , Ph.D. thesis, Fritz Haber Institute Berlin (2016).
  * Gould and Bučko (2016) T. Gould and T. Bučko, J. Chem. Theory Comput. 12, 3603 (2016).
  * Schwerdtfeger and Nagle (2019) P. Schwerdtfeger and J. K. Nagle, Mol. Phys. 117, 1200 (2019).
  * Johnson and Becke (2006) E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Jones  _et al._ (2013b) A. Jones, F. Cipcigan, V. P. Sokhan, J. Crain, and G. J. Martyna, Phys. Rev. Lett. 110, 227801 (2013b).
  * Porsev and Derevianko (2003) S. G. Porsev and A. Derevianko, J. Chem. Phys. 119, 844 (2003).
  * Porsev and Derevianko (2006) S. G. Porsev and A. Derevianko, J. Exp. Theor. Phys. 102, 195 (2006).
  * Jiang  _et al._ (2015) J. Jiang, J. Mitroy, Y. Cheng, and M. W. J. Bromley, Atom. Data Nucl. Data 101, 158 (2015).
  * Tao and Rappe (2016) J. Tao and A. M. Rappe, J. Chem. Phys. 144, 031102 (2016).
  * (49) Coupled-cluster response calculations for frequency-dependent multipole polarizabilities of open-shell atoms containing $p$, $d$, or $f$ valence orbitals undergo serious convergence problems when employing Gaussian basis sets (private communication with Prof. Dr. Robert A. DiStasio Jr., Cornell University). See also Ref. Che  and the related references there, for a discussion of difficulties in theoretical and experimental determination of multipole polarizabilities.
  * Stone (2013) A. Stone,  _The theory of intermolecular forces_ (OUP oxford, 2013).
  * Jeziorski  _et al._ (1994) B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994).
  * Kaplan (2006) I. G. Kaplan,  _Intermolecular interactions: physical picture, computational methods and model potentials_ (John Wiley & Sons, 2006).
  * von Lilienfeld  _et al._ (2004) O. A. von Lilienfeld, I. Tavernelli, and U. Rothlisberger, Phys. Rev. Lett. 93, 153004 (2004).
  * Van Vleet  _et al._ (2016) M. J. Van Vleet, A. J. Misquitta, A. J. Stone, and J. R. Schmidt, J. Chem. Theory Comput. 12, 3851 (2016).
  * (55) Often the term _polarization potential_ refers to a potential given by an expansion in terms of multipole polarizabilities Zhang  _et al._ (2012) to describe the long-range interaction between an atom and a charged particle. Here, we use this term for the function given by Eq. (7) which is another kind of polarization potential caused by an external field.
  * Szabó  _et al._ (2022) P. Szabó, S. Góger, J. Charry, M. R. Karimpour, D. V. Fedorov, and A. Tkatchenko, Phys. Rev. Lett. 128, 070602 (2022).
  * (57) See Supplemental Material at …
  * Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Epifanovsky  _et al._ (2021) E. Epifanovsky _et al._ , J. Chem. Phys. 155, 084801 (2021).
  * Lu and Chen (2012) T. Lu and F. Chen, J. Comput. Chem. 33, 580 (2012).
  * Zhang and Lu (2021) J. Zhang and T. Lu, Phys. Chem. Chem. Phys. 23, 20323 (2021).
  * Hait and Head-Gordon (2018) D. Hait and M. Head-Gordon, Phys. Chem. Chem. Phys. 20, 19800 (2018).
  * (64) The field strength defined this way is larger than the one arising in vdW-bonded systems – the field strengths generated by the electron density of free atoms at equilibrium vdW distance is on the order of $2\times 10^{-6}$ a.u. However, we also assessed a wider range of field strengths for three atoms, between $10^{-6}$ and $5\times 10^{-2}$ a.u., finding that our conclusions hold for this range as well SM .
  * Derevianko  _et al._ (2010) A. Derevianko, S. G. Porsev, and J. F. Babb, Atom. Data Nucl. Data 96, 323 (2010).
  * Tang and Toennies (2003) K. T. Tang and J. P. Toennies, J. Chem. Phys. 118, 4976 (2003).
  * (67) K. U. Lao, Y. Yang, and R. A. DiStasio, Jr., On the higher-order static polarizabilities and dispersion coefficients of the fullerenes: An ab initio study, doi:10.26434/chemrxiv.12940652.v1.
  * Zhang  _et al._ (2012) Y.-H. Zhang, L.-Y. Tang, X.-Z. Zhang, J. Jiang, and J. Mitroy, J. Chem. Phys. 136, 174107 (2012).


11institutetext:  Department of Metabolism, Digestion &Reproduction, Imperial
College London, UK 22institutetext: Department of Computing, Imperial College
London, UK 33institutetext: National Heart & Lung Institute, Imperial College
London, UK
33email<EMAIL_ADDRESS>
# PAT-CNN: Automatic Segmentation and Quantification of Pericardial Adipose
Tissue from T2-Weighted Cardiac Magnetic Resonance Images

Zhuoyu Li 11 Camille Petri 22 3 3 James Howard 33 Graham Cole 33 Marta Varela
33

###### Abstract

Background: Increased pericardial adipose tissue (PAT) is associated with many
types of cardiovascular disease (CVD). Although cardiac magnetic resonance
images (CMRI) are often acquired in patients with CVD, there are currently no
tools to automatically identify and quantify PAT from CMRI. The aim of this
study was to create a neural network to segment PAT from T2-weighted CMRI and
explore the correlations between PAT volumes (PATV) and CVD outcomes and
mortality.

Methods: We trained and tested a deep learning model, PAT-CNN, to segment PAT
on T2-weighted cardiac MR images. Using the segmentations from PAT-CNN, we
automatically calculated PATV on images from 391 patients. We analysed
correlations between PATV and CVD diagnosis and 1-year mortality post-imaging.

Results: PAT-CNN was able to accurately segment PAT with Dice score/ Hausdorff
distances of $0.74\pm 0.03/27.1\pm 10.9~{}mm$, similar to the values obtained
when comparing the segmentations of two independent human observers ($0.76\pm
0.06/21.2\pm 10.3~{}mm$). Regression models showed that, independently of sex
and body-mass index, PATV is significantly positively correlated with a
diagnosis of CVD and with 1-year all cause mortality (p-value <0.01).

Conclusions: PAT-CNN can segment PAT from T2-weighted CMR images automatically
and accurately. Increased PATV as measured automatically from CMRI is
significantly associated with the presence of CVD and can independently
predict 1-year mortality.

###### Keywords:

Cardiac MRI Segmentation Pericardial Fat Volume Cardiovascular Disease Risk
Factors Convolutional Neural Network

##  1 Introduction

Large volumes of pericardial adipose tissue, PAT, and endocardial adipose
tissue, EAT, predispose to a number of cardiovascular conditions [24] such as
atrial fibrillation [8]. This is due to the high metabolic activity of adipose
tissue which can stimulate remodelling in cardiac muscle, leading to fibrosis,
chamber enlargement and inflammation, among others. PAT can also regulate
cardiac tissue innervation [16] likely pro-arrhythmically. Obese patients
usually have large amounts of PAT, but correlations between cardiac fat volume
and body mass index (BMI) are relatively weak [8, 10].

So far, PAT has been predominantly segmented and quantified using computed
tomography (CT) [1]. In these images, fat’s low Hounsfield number allows for
comparatively straightforward identification and segmentation using simple
intensity threshold-based methods such as region growing and active contours
[22] or more complicated approaches like neural networks [5]. Fat is in
general more difficult to segment in Cardiac Magnetic Resonance Images (CMRI),
where its intensity is much more variable and similar to other tissues’.
Quantification of PAT volume (PATV) from MRI has so far only been reported
using time-consuming manual segmentation of PAT [6] or dedicated CMRI
sequences [7] such as Dixon imaging [10, 12]. This prevents widespread
quantification of PATV from more common types of CMRI, despite CMRI’s
increasing prominence in the diagnosis and assessment of cardiovascular
disease (CVD). A tool to rapidly and accurately quantify PAT from CMRI would
allow rapid quantification of PATV in the clinic, to aid CVD risk assessment.

MRI segmentation Convolutional neural networks (CNNs) have shown great ability
to segment MR images, enabling easy quantification of several imaging
biomarkers whose analysis would otherwise be too time consuming. They
consistently outperform non-learned techniques such as thresholding [20],
region growing [19] and probability graph models [9]. Once trained, CNNs can
perform image analysis in real time in the clinic. To date, and to the best of
our knowledge, there are no CNN-based tools for the segmentation of PAT from
T2-weighted CMR images.

Aims In this project, we aim to estimate PATV from clinical CMRI. We focus on
axial T2-weighted rapid gradient-echo scans which are routinely acquired in
the clinic and in which CNNs have been used to accurately segment several
cardiac structures [11]. We initially create a non-learned semi-automatic
image processing protocol to segment PAT and generate training data for the
proposed CNN. We then train and test a dedicated CNN to segment PAT in these
images. Finally, using data from 391 patients, we investigate correlations
between: a) PATV and the presence of CVD and b) PATV and 1-year all-cause
mortality.

##  2 Methods

Data We uses fully anonymized T2-prepared spoiled GRE axial stacks acquired at
1.5T (TE: 1.56 ms, TR: 3.1 ms, recon voxel size: 1.4 x 1.4 x 6.0 $mm^{3}$) in
391 patients with known or suspected CVD, under ethical approval. The number
of slices in each volumetric scan ranged from 23 to 55.

Preprocessing Ground truth PAT segmentation was semi-automatically performed
in 70 training images (Fig 1). As the first step, we identified the cardiac
chambers and left ventricular (LV) myocardium using HRNet, an existing CNN
trained in T2-weighted images with a similar contrast [11]. Using these
segmentations, we created a bounding box around the heart. We then performed a
2-class Otsu segmentation [21] on the cropped images. In most cases, one of
the classes roughly coincided with subcutaneous fat, allowing the rapid
segmentation of PAT after some further manual corrections. These semi-
automatic PAT segmentations were reviewed by a second independent observer.

Figure 1: Semi-automatic pipeline for ground truth PAT segmentation. ROI:
Region of Interest

Model: Architecture To allow fully automatic PAT labelling, we trained PAT-
CNN, a 3D Res-UNet CNN [13]. The network (shown in Fig 2) included a
contracting path with four convolutional layers, an expanding path with four
upsampling layers and concatenate connections. A residual unit was added to
each layer. In residual units, batch normalization and parametric rectified
linear unit (PReLU) were implemented to improve the feature selection results.
The output layer used a sigmoid function as the activation function and was
trained to classify each pixel as: PAT (label 1) or background/other tissues
(label 0).

Figure 2: Structure diagram of PAT-CNN. The kernel sizes of convolutions and
transpose convolutions are 3 x 3 x 3. The numbers represent filters numbers in
each layer. The triangle arrows represents the residual units.

Model: Dataset Setup 70 scans with manual PAT labels were randomly divided
into a training+validation dataset (60 images) and a testing dataset (10
images). A second independent observer independently segmented PAT in the test
dataset to assess inter-observer segmentation metrics.

Model: Data Augmentation To enhance PAT-CNN’s invariance and generalization,
as well as expand the size of the dataset, we applied data augmentation during
training using the image analysis library SimpleITK [25]. We randomly applied
rotation, crop, flip, blur and Rayleigh noise addition to each input image of
the network in each epoch. The probability of each transformation was 0.1.

Model: Loss Function The loss function combined cross-entropy loss and Dice
loss terms [26]. The Dice loss term help improves PAT-CNN’s accuracy in the
presence of the large imbalance between the volumes of PAT and voxels labelled
as background/other organs. The cross-entropy loss term minimises the likely
instability of the gradients of a pure Dice-score-based loss function. The
total loss function can be represented as below, where $x_{n}$ is the ground
truth label of the $n^{th}$ sample, $y_{n}$ is its prediction, N is the total
number of samples (and the small term $\epsilon$ in the Dice coefficient is
used to avoid the denominator becoming 0):

$Loss(X,Y)=1-\frac{2\sum_{n=1}^{N}x_{n}y_{n}+\epsilon}{\sum_{n=1}^{N}x_{n}+\sum_{n=1}^{N}y_{n}+\epsilon}-\frac{1}{N}\sum_{n=1}^{N}x_{n}log(y_{n})$
(1)

Model: Hyperparameters and Performance Assessment Through our experimental
results, the optimal batch size was 8. We used Adam [15] and set the learning
rate to $10^{-4}$, training for 60 epochs. This combination performed best on
the validation dataset. We used the Dice Score and Hausdorff Distance to
evaluate PAT-CNNs performance on the test dataset. Using these metrics, we
additionally compared PAT-CNN to the vanilla U-Net [23] commonly using for
segmenting medical images.

Statistical Analysis: Data We used PAT-CNN to estimate PATV (in $cm^{3}$) in
all 391 volumetric T2-weighted axial CMR images. From the anonymized imaging
metadata and radiology reports, we extracted the following demographic
information ($mean\pm std$): age ($55\pm 18$ years), sex ($42\%$ female), BMI
($27.7\pm 5.9~{}kg/m^{2}$) and 1-year post-scan mortality (‘Deceased’,
$8.2\%$). We also compiled diagnostic information, which was synthesised into
73 binary clinical diagnostic labels (presence/absence of: dilated
cardiomyopathy, left/right ventricular failure, left/right ventricular
hypertrophy, myocardial infarction, ischaemia, etc). The sum of these
diagnostic labels was compiled into a single variable ‘CVD Diagnosis’ ($2.9\pm
1.9$).

Statistical Analysis: Methods We investigated the relationship between the
regressor variables: PATV, sex, age and BMI and the dependent variables: ‘CVD
diagnosis’, ‘Deceased’ and PATV, using univariate Pearson/phi correlation
analysis. Those regressors that were significantly correlated ($\alpha$ =
0.01) with the independent variables were then used as inputs to a
multivariate analysis using ordinary least squares (for the ‘CVD Diagnosis’
and ‘PATV’ dependent variables) and a logistic regression model (for the
dichotomous variable ‘Deceased’). We note that the clinical dataset used did
not include time-to-mortality information and survival analyses were therefore
not performed.

##  3 Results

PAT-CNN Performance On the test dataset, PAT was segmented with an average
Dice score/Hausdorff distance of $0.74\pm 0.03/27.1\pm 10.9~{}mm$. The model
coped well with variations in PAT shape, size and location (Fig 3) and the
presence of pathological changes. PAT-CNN performed better than U-Net and its
results compared well with the agreement between two human operators on the
same dataset (Table 1).

| PAT-CNN | U-Net | Observer 2
---|---|---|---
Dice Score | $0.74\pm 0.03$ | $0.71\pm 0.03$ | $0.76\pm 0.06$
Hausdorff Distance (mm) | $27.1\pm 10.9$ | $32.8\pm 11.2$ | $21.2\pm 10.3$

Table 1: Performance of PAT-CNN, U-Net and a second human observer on the test
images. (The segmentations from human Observer 1 were treated as the ground
truth for this analysis.)

Statistical Analysis Across the 391 analysed patients, PATV had a mean of
$139.60~{}cm^{3}$ and a standard deviation of $80.24~{}cm^{3}$, in excellent
agreement with PAT quantification from CT [5]. PATV is significantly
correlated with patients’ age ($\rho_{Pearson}=0.38$), sex and BMI
($\rho_{Pearson}=0.34$), but in a multivariate analysis, only sex and BMI
remain significant (Table 2). PATV and age were the only regressor variables
found to be significantly correlated to both CVD diagnostic labels and 1-year
mortality post-imaging in both univariate and multivariate analyses (Table 2).
This is in contrast to BMI which was not an independent predictor of either
mortality or CVD diagnosis in univariate analysis (Table 2).

| PATV ($cm^{3}$) |  Deceased | CVD Diagnosis
---|---|---|---
Sex |  -50.2* |  |  -0.30
Age ($years$) |  1.70 |  0.03* |  0.02*
BMI ($kg/m^{2}$) |  3.92* |  |
PATV ($cm^{3}$) |  |  0.01* |  0.01*
Table 2: Regression coefficients from 3 multivariate regression analyses (n =
391), for the dependent variables listed in each column.
Each row shows each regressor considered in the model. (For the dichotomous
variable ‘Sex’, females are represented by 1 and males by 0.) Only regressors
found to be statistically significant in univariate analyses were considered.
Entries are the regression coefficients in the specified units, with the
asterix (*) denoting significant correlations at $\alpha=0.01$. PATV:
pericardial adipose tissue volume; BMI: body-mass index; CVD: cardiovascular
disease.

##  4 Discussion

We present PAT-CNN, a CNN able to identify PAT accurately in T2-weighted CMRI
across subjects of various sizes and in the presence of different types of
pathology. PATV estimated using PAT-CNN’s segmentation of CMRI is
significantly correlated to the presence of cardiac pathology and is an
independent predictor of 1-year mortality, as seen in studies from the
literature in which PAT was segmented manually [8].

Figure 3: Representative examples of T2-weighted slices from four different
patients overlaid with the contours for the ground truth (red) and PAT-CNN
(green) PAT labels.

PAT Segmentation PAT is distributed around in the heart in several non-
contiguous patches with various shapes, making PAT segmentation more complex
than the segmentation of contiguous and less variable structures such as the
cardiac chambers themselves. Moreover, when present, pericardial fluid can
easily be mistakenly labelled as PAT, as fluid and fat have a similar
intensity and appearance on T2-weighted CMR images.

Despite the difficulty of the task, PAT-CNN was able to segment PAT with an
average Dice score of 0.74, outperforming U-Net and performing only marginally
worse than the average Dice score of 0.76 obtained when comparing the
segmentations from different human observers (see Table 1).

PAT-CNN uses a 3D Res-UNet previously used to segment heart chambers [14, 18]
, with a modified combined Dice-Cross-Entropy loss that outperformed Dice-only
and Cross-Entropy-only losses for this task. We therefore show that the 3D
Res-UNet has a good performance in CMRI even when segmenting irregular, non-
connected regions such as PAT.

Statistical Analysis of PATV

PATV estimated automatically using PAT-CNN shows similar correlations to
clinical variables as PATV estimated from CT and echocardiography [2, 17].
Females, elderly people and those with higher BMI are more likely to have more
PATV. Statistically significant correlations were found between PATV and a
diagnosis of cardiovascular disease at the time of scan, as has been reported
in other studies [2]. PATV was also an independent predictor of one-year all-
cause mortality in the analysed 319 subjects, as also found on CT [4].
Interestingly, BMI was not a predictor for mortality or the presence of CVD,
although it was significantly but weakly correlated with PAT
($\rho_{Pearson}=0.34$), also as seen in other studies [2, 10].

Our analysis used 1-year all-cause mortality, as we could not obtain
information about CVD-related deaths or temporal information about the time to
death. Body surface area (BSA) is more commonly used than BMI for the indexing
of cardiac function variables. BSA would have been an interesting variable to
include in our regression models, but it was not available in our dataset.

Future Work PAT segmentations, such as those provided by PAT-CNN, provide
information beyond that of PATV alone. In the current study, we did not
analyse the clinical importance of the location or spatial distribution of
PAT, which is likely to have clinical value beyond PAT volume only. We also
did not attempt to index PATV by body area or heart size, which may yield
biomarkers with further clinical value [3].

Alongside PAT, epicardial adipose tissue (EAT) is also likely to play an
important role in cardiovascular pathophysiology and it should be identifiable
in CMRI. EAT, however, is expected to deposit in regions of smaller volume and
its identification is limited by partial volume effects at MRI’s typical
spatial resolution. EAT’s reliable identification would require MRI at much
higher spatial resolution or the use of specialised CMRI sequences such as
Dixon techniques [10, 12]. Dixon technique CMR images are not usually acquired
as part of clinical examinations, stressing the importance of our approach to
automatically segment PAT from commonly acquired T2-weighted clinical scans.

Pending training and testing across images from different sites and scanners,
PAT-CNN can be deployed as a scanner-side analysis tool to automatically
estimate PATV. It can also be used to explore the mechanisms of different CVDs
by highlighting potential relationships between PAT volume and location and
CVD manifestations.

##  5 Conclusions

We propose PAT-CNN, a neural network to automatically segment PAT from CMR
images. We show that PAT segmentations are reliable, making PAT-CNN ready to
be deployed on clinical images to yield biomarkers of potential clinical
interest. Using data from 391 patients, we find that PAT volume, estimated
using PAT-CNN, is significantly correlated with both a diagnosis of CVD and
1-year all-cause mortality, independently of sex and BMI.

## Acknowledgements

This work was supported by the British Heart Foundation Centre of Research
Excellence at Imperial College London (RE/18/4/34215).

## References

  * [1] Benčević, M., Galić, I., Habijan, M., Pižurica, A.: Recent Progress in Epicardial and Pericardial Adipose Tissue Segmentation and Quantification Based on Deep Learning: A Systematic Review. Applied Sciences 12(10), 5217 (2022). https://doi.org/10.3390/app12105217
  * [2] Britton, K.A., Massaro, J.M., Murabito, J.M., Kreger, B.E., Hoffmann, U., Fox, C.S.: Body Fat Distribution, Incident Cardiovascular Disease, Cancer, and All-Cause Mortality. Journal of the American College of Cardiology 62(10), 921–925 (9 2013). https://doi.org/10.1016/J.JACC.2013.06.027
  * [3] Cai, S., Wald, R., Deva, D.P., Kiaii, M., Ng, M.Y., Karur, G.R., Bello, O., Li, Z.J., Leipsic, J., Jimenez-Juan, L., Kirpalani, A., Connelly, K.A., Yan, A.T.: Cardiac MRI measurements of pericardial adipose tissue volumes in patients on in-centre nocturnal hemodialysis. Journal of Nephrology 33(2), 355–363 (4 2020). https://doi.org/10.1007/s40620-019-00665-4
  * [4] Cheng, V.Y., Dey, D., Tamarappoo, B., Nakazato, R., Gransar, H., Miranda-Peats, R., Ramesh, A., Wong, N.D., Shaw, L.J., Slomka, P.J., Berman, D.S.: Pericardial Fat Burden on ECG-Gated Noncontrast CT in Asymptomatic Patients Who Subsequently Experience Adverse Cardiovascular Events. JACC: Cardiovascular Imaging 3(4), 352–360 (4 2010). https://doi.org/10.1016/J.JCMG.2009.12.013
  * [5] Commandeur, F., Goeller, M., Betancur, J., Cadet, S., Doris, M., Chen, X., Berman, D.S., Slomka, P.J., Tamarappoo, B.K., Dey, D.: Deep Learning for Quantification of Epicardial and Thoracic Adipose Tissue from Non-Contrast CT. IEEE Transactions on Medical Imaging 37(8), 1835–1846 (8 2018). https://doi.org/10.1109/TMI.2018.2804799
  * [6] Davidovich, D., Gastaldelli, A., Sicari, R.: Imaging cardiac fat. European Heart Journal - Cardiovascular Imaging 14(7), 625–630 (7 2013). https://doi.org/10.1093/EHJCI/JET045
  * [7] Ding, X., Pang, J., Ren, Z., Diaz-Zamudio, M., Jiang, C., Fan, Z., Berman, D.S., Li, D., Terzopoulos, D., Slomka, P.J., Dey, D.: Automated pericardial fat quantification from coronary magnetic resonance angiography: feasibility study. Journal of Medical Imaging 3(1), 014002 (2 2016). https://doi.org/10.1117/1.jmi.3.1.014002, https://doi.org/10.1117/1.jmi.3.1.014002
  * [8] Fitzgibbons, T.P., Czech, M.P.: Epicardial and Perivascular Adipose Tissues and Their Influence on Cardiovascular Disease: Basic Mechanisms and Clinical Associations. Journal of the American Heart Association 3(2) (3 2014). https://doi.org/10.1161/JAHA.113.000582
  * [9] Grosgeorge, D., Petitjean, C., Dacher, J.N., Ruan, S.: Graph cut segmentation with a statistical shape model in cardiac MRI. Computer Vision and Image Understanding 117(9), 1027–1035 (2013). https://doi.org/https://doi.org/10.1016/j.cviu.2013.01.014
  * [10] Henningsson, M., Brundin, M., Scheffel, T., Edin, C., Viola, F., Carlhäll, C.J.: Quantification of epicardial fat using 3D cine Dixon MRI. BMC Medical Imaging 20(1), 1–9 (7 2020). https://doi.org/10.1186/s12880-020-00478-z
  * [11] Howard, J., Zaman, S., Ragavan, A., Hall, K., Leonard, G., Sutanto, S., Ramadoss, V., Razvi, Y., Linton, N.F., Bharath, A., Shun-Shin, M.J., Rueckert, D., Francis, D., Cole, G.: Automated analysis and detection of abnormalities in transaxial anatomical cardiovascular magnetic resonance images: a proof of concept study with potential to optimize image acquisition. International Journal of Cardiovascular Imaging 37(3), 1033–1042 (2021). https://doi.org/10.1007/s10554-020-02050-w
  * [12] Kellman, P., Hernando, D., Shah, S., Zuehlsdorff, S., Jerecic, R., Mancini, C., Liang, Z.P., Arai, A.E.: Multiecho dixon fat and water separation method for detecting fibrofatty infiltration in the myocardium. Magnetic Resonance in Medicine 61(1), 215–221 (1 2009). https://doi.org/10.1002/mrm.21657
  * [13] Kerfoot, E., Clough, J., Oksuz, I., Lee, J., King, A.P., Schnabel, J.A.: Left-Ventricle Quantification Using Residual U-Net. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 11395 LNCS, 371–380 (2019). https://doi.org/10.1007/978-3-030-12029-0_40
  * [14] Kerfoot, E., Puyol Anton, E., Ruijsink, B., Clough, J., King, A.P., Schnabel, J.A.: Automated CNN-based reconstruction of short-axis cardiac MR sequence from real-time image data. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). vol. 11040 LNCS, pp. 32–41. Springer Verlag (2018). https://doi.org/10.1007/978-3-030-00946-5_4
  * [15] Kingma, D., Ba, J.: Adam: A Method for Stochastic Optimization. Computer Science (2014)
  * [16] Lavie, C.J., Pandey, A., Lau, D.H., Alpert, M.A., Sanders, P.: THE PRESENT AND FUTURE Obesity and Atrial Fibrillation Prevalence, Pathogenesis, and Prognosis Effects of Weight Loss and Exercise. Tech. rep. (10 2017). https://doi.org/10.1016/J.JACC.2017.09.002
  * [17] Liu, J., Fox, C.S., Hickson, D.M., Sarpong, D., Ekunwe, L., May, W.D., Hundley, G.W., Carr, J.J., Taylor, H.A.: Pericardial Adipose Tissue, Atherosclerosis, and Cardiovascular Disease Risk FactorsThe Jackson Heart Study. Diabetes Care 33(7), 1635–1639 (7 2010). https://doi.org/10.2337/DC10-0245
  * [18] Lourenço, A., Kerfoot, E., Dibblin, C., Alskaf, E., Anjari, M., Bharath, A.A., King, A.P., Chubb, H., Correia, T.M., Varela, M.: Left Atrial Ejection Fraction Estimation Using SEGANet for Fully Automated Segmentation of CINE MRI. In: Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics). vol. 12592 LNCS, pp. 137–145. Springer Science and Business Media Deutschland GmbH (2021). https://doi.org/10.1007/978-3-030-68107-4_14, https://link.springer.com/chapter/10.1007/978-3-030-68107-4_14
  * [19] Mancas, M., Gosselin, B., Macq, B.: Segmentation using a region-growing thresholding. In: Dougherty, E.R., Astola, J.T., Egiazarian, K.O. (eds.) Image Processing: Algorithms and Systems IV. vol. 5672, pp. 388–398. International Society for Optics and Photonics, SPIE (2005). https://doi.org/10.1117/12.587995
  * [20] Norouzi, A., Rahim, M.S.M., Altameem, A., Saba, T., Rad, A.E., Rehman, A., Uddin, M.: Medical Image Segmentation Methods, Algorithms, and Applications. IETE Technical Review 31(3), 199–213 (2014). https://doi.org/10.1080/02564602.2014.906861
  * [21] Otsu, N.: A Threshold Selection Method from Gray-Level Histograms. IEEE Transactions on Systems, Man, and Cybernetics 9(1), 62–66 (1979). https://doi.org/10.1109/TSMC.1979.4310076
  * [22] Rodrigues, O., Morais, F.F., Morais, N.A., Conci, L.S., Neto, L.V., Conci, A.: A novel approach for the automated segmentation and volume quantification of cardiac fats on computed tomography. Computer Methods and Programs in Biomedicine 123, 109–128 (1 2016). https://doi.org/10.1016/J.CMPB.2015.09.017
  * [23] Ronneberger, O., Fischer, P., Brox, T.: U-Net: Convolutional Networks for Biomedical Image Segmentation. pp. 234–241. Springer, Cham (2015). https://doi.org/10.1007/978-3-319-24574-4_28, http://link.springer.com/10.1007/978-3-319-24574-4_28
  * [24] Weschenfelder, C., Quadros, A., Santos, J., Garofallo, S.B., Marcadenti, A.: Adipokines and Adipose Tissue-Related Metabolites, Nuts and Cardiovascular Disease. Metabolites 10(1), 32– (2020)
  * [25] Yaniv, Z., Lowekamp, B.C., Johnson, H.J., Beare, R.: {SimpleITK} Image-Analysis Notebooks: a Collaborative Environment for Education and Reproducible Research. Journal of Digital Imaging 31(3), 290–303 (2018). https://doi.org/10.1007/s10278-017-0037-8
  * [26] Yeung, M., Sala, E., Schönlieb, C.B., Rundo, L.: Unified Focal loss: Generalising Dice and cross entropy-based losses to handle class imbalanced medical image segmentation. Computerized Medical Imaging and Graphics 95, 102026 (1 2022). https://doi.org/10.1016/J.COMPMEDIMAG.2021.102026