# A Generalized Lerche-Newberger Formula

Parker Kuklinski  Michael Warnock  David Hague

###### Abstract

The Lerche-Newberger formula simplifies harmonic sums of Bessel functions and
has seen application in plasma physics and frequency modulated quantum
systems. In this paper, we rigorously prove the formula and extend the
classical result to a family of multi-dimensional extensions of the single
variable Bessel functions called generalized Bessel functions. Since
prevailing definitions of these functions do not accommodate arbitrary complex
order, we use an auxiliary family of functions called generalized Anger
functions and show that the single-variable result holds in multiple
dimensions for a certain selection of parameters. We conclude by applying
these results to physical systems.

##  1 Introduction

Sums involving Bessel functions are pervasive in physics [1] [2] [3]. A
Jacobi-Anger expansion, for instance, allows for a Fourier representation of a
sinusoidal frequency modulated signal as a sum of harmonic components weighted
by Bessel functions of varying order [4]. Other sums like the Kapteyn series
arise in astrophysics applications [5], while Neumann and Schlömilch-type
series also also appear in the literature [6] [7]. The Lerche-Newberger
formula [8] gives a closed form representation to an infinite sum of Bessel
function products with a transposed harmonic term (i.e. $1/(n-a)$). These sums
first arose in 3-dimensional plasma systems with an oscillating ambient
magnetic field from the associated $3\times 3$ plasma wave dispersion relation
[9] [10]. In the present context we are interested in applying the Lerch-
Newberger formula to quantum systems under frequency modulation [11]. Berns
et. al. [12] investigate an intermediate ”quasiclassical” regime in a driven
two level quantum system, and the transition rate of the qubit is modeled by a
modified Lerche-Newberger summation.

In this paper we aim to adapt the methods of Kuklinski and Hague [13] to
extend the traditional Lerche-Newberger formula to a class of multi-
dimensional Bessel functions called generalized Bessel functions [14].
Generalized Bessel functions introduce higher order harmonics in the
modulation function of the associated integral. In a signal processing
context, the generalized Bessel functions are the Fourier coefficients of a
multi-tone sinusoidal frequency modulated (MTSFM) signal; these signals are
useful for radar and sonar applications due to their constant amplitude,
spectral efficiency, and tunability [15]. The parameters of these MTSFM are
fed through to the arguments of the GBFs in the Fourier transformed space. The
generalized Bessel functions have found applications in laser physics [16],
crystallography [17], and astrophysics [18]. In connection to frequency
modulated quantum systems, a version of the Lerche-Newberger sum found in
Berns et. al. with Bessel functions replaced by GBFs arises in a more general
oscillatory system.

A naïve extension of the Lerche-Newberger formula by a total substitution of
one-dimensional Bessel functions with their multi-dimensional counterparts
fails since the generalized Bessel functions are not well defined for
fractional order. To remedy this, we instead work with generalized Anger
functions which extend the usual integral representation of integer-order
generalized Bessel functions to arbitrary order. These Anger functions,
although agreeing with Bessel functions on integer order, do not satisfy the
same collection of identities and thus certain unresolvable terms emerge in
the corresponding Lerche-Newberger extension.

The rest of this paper is organized as follows: In section 2 we rigorously
derive the one-dimensional Lerche-Newberger formula over the usual domain of
parameters. In section 3 these results are extended by replacing the one-
dimensional Bessel functions with generalized Anger functions. Section 4
treats an application to a multiply frequency modulated two-level quantum
system. An appendix handles some details of the proof of the one-dimensional
Lerche-Newberger formula.

##  2 Derivation of Lerche-Newberger formula

In this section we will derive the following sum presented by Newberger [19]
[8]:

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}J_{\alpha-\gamma n}(z)J_{\beta+\gamma
n}(z)}{n+\mu}=\frac{\pi}{\sin\mu\pi}J_{\alpha+\gamma\mu}(z)J_{\beta-\gamma\mu}(z)$
(1)

Here, we restrict $\text{Re}(\alpha+\beta)>-1$,
$\mu\in\mathbb{C}\backslash\mathbb{Z}$, and $\gamma\in(0,1]$. To proceed, we
follow an argument from Lerche [20]; recall from Watson [3] that for
$\text{Re}(\alpha+\beta)>-1$ we have

$J_{\alpha}(z)J_{\beta}(z)=\frac{2}{\pi}\int_{0}^{\pi/2}J_{\alpha+\beta}(2z\cos\theta)\cos(\alpha-\beta)\theta
d\theta.$ (2)

Substituting this into the summation on the left side of (1), we have

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}J_{\alpha-\gamma n}(z)J_{\beta+\gamma
n}(z)}{n+\mu}=\frac{2}{\pi}\sum_{n=-\infty}^{\infty}\left[\int_{0}^{\pi/2}\frac{(-1)^{n}}{n+\mu}J_{\alpha+\beta}(2z\cos\theta)\cos{(\alpha-\beta-2\gamma
n)\theta}d\theta\right].$ (3)

Ultimately we would like to interchange summation and integration, and we do
this with dominated convergence theorem [21]. Since the doubly-infinite
summation is a composition of two separate limiting operations, one over the
negative indices and one over the positive indices, we fmust split the double
summation into two one-sided summations. one we label $S_{+}$ with
$n\in\\{1,2,...\\}$ and the other we call $S_{-}$ with $n\in\\{-1,-2,...\\}$
such that

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}J_{\alpha-\gamma n}(z)J_{\beta+\gamma
n}(z)}{n+\mu}=\frac{2}{\pi}\left(S_{+}+S_{-}+\frac{1}{\mu}\int_{0}^{\pi/2}J_{\alpha+\beta}(2z\cos{\theta})\cos{(\alpha-\beta)\theta}d\theta\right).$
(4)

To apply dominated convergence theorem to the summations $S_{\pm}$, we first
state a result for two simpler sequence of functions. Let $f_{n}(\theta)$ and
$h_{n}(\theta)$ be the sequences of partial sums

$f_{n}(\theta)=\sum_{k=1}^{n}\frac{(-1)^{k}}{k+\mu}\cos(k\theta),\hskip
28.45274ptg_{n}(\theta)=\sum_{k=1}^{n}\frac{(-1)^{k}}{k+\mu}\sin(k\theta)$ (5)

for $\theta$ in some finite interval in $\mathbb{R}$. These partial sums are
related to Lerch zeta functions [22]. We are searching for functions
$F(\theta)$ and $G(\theta)$ independent of $n$ such that
$|f_{n}(\theta)|<F(\theta)$ and $|g_{n}(\theta)|<G(\theta)$ for all $n$ and
$F(\theta)$ and $G(\theta)$ are integrable. We can use a triangle inequality
to bound $|f_{n}(\theta)|$ in the following way:

$\displaystyle|f_{n}(\theta)|$
$\displaystyle\leq\left\lvert\sum_{k=1}^{n}\frac{(-1)^{k}}{k+\mu}\cos(k\theta)-\sum_{k=1}^{n}\frac{(-1)^{k}}{k}\cos(k\theta)\right\rvert+\left\lvert\sum_{k=1}^{n}\frac{(-1)^{k}}{k}\cos(k\theta)\right\rvert$
$\displaystyle=|\mu|\left\lvert\sum_{k=1}^{n}\frac{(-1)^{k}}{k(k+\mu)}\cos(k\theta)\right\rvert+\left\lvert\sum_{k=1}^{n}\frac{(-1)^{k}}{k}\cos(k\theta)\right\rvert$
(6)

If $\mu$ isn’t a negative integer, then by a $p$-series test the first sum on
the right hand side pf (6) can be bounded by a finite quantity. The second sum
on the left is a partial sum of a logarithm, specifically on the open interval
$(-\pi,\pi)$ we have the pointwise convergence

$\sum_{k=1}^{\infty}\frac{(-1)^{k}}{k}\cos(k\theta)=-\log\left(\cos\frac{\theta}{2}\right)-\log{2}$
(7)

We prove that the partial sum term in (6) is bounded above by its pointwise
limit plus a constant, and is therefore bounded by an integrable function
letting us use dominated convergence theorem. We unfortunately cannot apply
standard results in Gibbs phenomenon [23] as these results apply only to
discontinuous functions with finite jumps; the one we have here is an infinite
discontinuity. Nevertheless, we can exploit the specific nature of our
function to arrive at an appropriate bound. Let $s_{n}(\theta)$ be the
$n^{\text{th}}$ partial sum such that
$s_{n}(\theta)+\log\left(\cos\theta/2\right)+\log{2}$ has local extrema in
this interval at points $\theta_{j}=2\pi j/(2n+1)$. We prove in the appendix
that there exists some $M\in\mathbb{N}$ such that for all $n>M$ the largest
extrema of this difference occurs at $\theta_{n}=\pi-\pi/(2n+1)$ (that this
holds only for $n$ sufficiently large does not impede the conditions of the
DCT since for all $n<M$ we can bound (6) by a constant). Therefore by plugging
$\theta_{n}$ into the difference of these terms we have

$s_{n}(\theta_{n})+\log\left(\cos\frac{\theta_{n}}{2}\right)+\log{2}=\left[\sum_{k=1}^{n}\frac{1}{k}\cos\left(\frac{k\pi}{2n+1}\right)\right]+\log\left(\sin\frac{\pi}{4n+2}\right)+\log
2$ (8)

We show that as $n$ gets large, both the partial sum term and the logarithm
term grow at $O(\log{n})$ with opposite leading coefficients, thus cancelling
and leaving an O(1) term. To handle the partial sum term, use a trigonometric
identity to extract a harmonic sum:

$\sum_{k=1}^{n}\frac{1}{k}\cos\left(\frac{k\pi}{2n+1}\right)=\sum_{k=1}^{n}\frac{1}{k}-\sum_{k=1}^{n}\frac{2}{k}\sin^{2}\left(\frac{\pi
k}{4n+2}\right)$ (9)

By applying a standard result on harmonic series and using the identity $\sin
x<x$ for all $x>0$, we see that this partial sum on the left in (9) is equal
to $\log{n}+O(1)$. Next, we can conduct an asymptotic expansion of the
logarithm term in (8). Note that for small $x$ we have $\sin{x}=x(1+O(x^{2}))$
such that the logarithm pulls out this factor of $x$ as well as other
multiplicative factors and therefore:

$\log\left(\sin\frac{\pi}{4n+2}\right)=-\log{n}+O(1)$ (10)

The two asymptotic expansions we’ve conducted imply that the quantity in (8)
converges to a limit, or the partial sum in (6) overshoots its limit by an
asymptotically finite quantity, and therefore there exists some finite $M$
such that

$\left\lvert\sum_{k=1}^{n}\frac{(-1)^{k}}{k}\cos(k\theta)\right\rvert\leq\frac{1}{2}\left\lvert\log\left(1+\cos\theta\right)+\log{2}\right\rvert+M$
(11)

for all $n\in\mathbb{N}$ (Numerical simulations appear to show $M=1/2$ is a
sufficient choice). This function on the right is Lebesgue integrable on any
interval in $\mathbb{R}$ even in those that include singularities at
$\theta=(2N+1)\pi$. Therefore, we can apply dominated convergence theorem to
this sequence of partial sums $f_{n}(\theta)$. Applying dominated convergence
theorem to $\\{g_{n}(\theta)\\}$ is similar, but instead of bounding the
alternating cosine summation by a diverging yet integrable function, we can
bound the corresponding alternating sine function by a constant since this is
the Fourier transform of the sawtooth wave.

To extend this to a partial sum of interest in $S_{+}$, we need to bound
$J_{\alpha+\beta}(2z\cos\theta)[\cos(\alpha-\beta)\theta
f_{n}(2\gamma\theta)+\sin(\alpha-\beta)\theta g_{n}(2\gamma\theta)]$ by a
single integrable function for all $n$. But since Bessel functions and sine
and cosine functions are entire on $\mathbb{C}$, we can bound these and
$g_{n}(2\gamma\theta)$ by constants, and we can bound $f_{n}(2\gamma\theta)$
by the scaled version of the dominating function on the right hand side of
(11). Thus, the interchange of summation and integration is justified. These
arguments are directly applicable to $S_{-}$, however we also need to restrict
$\mu$ from being a positive integer, and since a $1/\mu$ term appears on the
right hand side of (4), we insist that $\mu\in\mathbb{C}\backslash\mathbb{Z}$.

Because dominated convergence theorem holds separately for $S_{\pm}$, there is
no issue with now interchanging the entire double summation with integration
in (3). We use the following identities to reduce this operation to a closed
form

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{n+\mu}\cos(n\theta)=\frac{\pi\cos{\mu\theta}}{\sin{\pi\mu}},\hskip
28.45274pt\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{n+\mu}\sin(n\theta)=-\frac{\pi\sin{\mu\theta}}{\sin{\pi\mu}}$
(12)

where $\theta\in[-\pi,\pi]$. In other words, the left hand side summations are
the Fourier series for non-smooth $2\pi$-periodic functions. To apply these
formulas to (4), we require that $\gamma\in(0,1]$ such that
$2\gamma\theta\in[-\pi,\pi]$ for the region of integration
$\theta\in[0,\pi/2]$. Therefore we write the interchange as

$\frac{2}{\pi}\sum_{n=-\infty}^{\infty}\left[\int_{0}^{\pi/2}\frac{(-1)^{n}}{n+\mu}J_{\alpha+\beta}(2z\cos\theta)\cos{(\alpha-\beta-2\gamma
n)\theta}d\theta\right]=\frac{2}{\sin{\pi\mu}}\int_{0}^{\pi/2}J_{\alpha+\beta}(2z\cos\theta)\cos{(\alpha-\beta+2\gamma\mu)\theta}d\theta$
(13)

Since the right hand side is in the integral form of the Bessel product
formula in (2), we can reverse the equation to arrive at the final result in
(1).

##  3 Extension to generalized Anger functions

In this section we create a generalized version of the Lerche-Newberger
formula in (1) extended to generalized Anger functions. Traditionally, the
two-dimensional generalized Bessel function is defined as the convolutional
summation:

$J_{n}(x,y)=\sum_{k=-\infty}^{\infty}J_{n-2k}(x)J_{k}(y)$ (14)

We usually call this function the index $(1,2)$ generalized Bessel function.
If we naively attempt to apply this definition to the Lerche-Newberger formula
in (1), we will run into problems because this form of the GBF is not defined
for $n\in\mathbb{C}\backslash\mathbb{Z}$ as the summation in (14) will not
converge. We could try to define the GBF in a completely analogous way as the
Bessel functions by first defining a partial differential equation and
defining the GBF as solutions to this PDF, but as seen in Kuklinski and Hague
[13], the index $(1,2)$ GBF satisfies two independent second-order linear
PDEs, neither of which have clear extensions to GBFs of different index or
higher order.

Regardless of these issues, however, if we restrict a putative Lerche-
Newberger formula to GBFs of integer order, we will have a well-defined
convergent sum. To do this, we introduce generalized Anger functions which are
well defined for all orders but agree with the GBF for integer order. For two
finite coordinates ${\bf x}=(x_{1},...,x_{m})$ and ${\bf y}=(y_{1},...,y_{m})$
in $\mathbb{C}^{m}$, we define the generalized Anger function of order
$\alpha$ as:

$A_{\alpha}({\bf x},{\bf
y})=\frac{1}{2\pi}\int_{-\pi}^{\pi}\exp\left[i\left(\alpha\theta-\sum_{k=1}^{m}\left(x_{k}\sin{k\theta}+y_{k}\cos{k\theta}\right)\right)\right]d\theta$
(15)

These functions have favorable decay properties in the sense that a stationary
phase approximation shows that for fixed ${\bf x},{\bf y}\in\mathbb{C}^{m}$,
$A_{\alpha}({\bf x},{\bf y})=O(\alpha^{-k})$ for all $k>0$ as
$|\alpha|\rightarrow\infty$.

With these functions, we prove a new version of the Lerche-Newberger formula:

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}A_{\alpha-\gamma n}({\bf x},{\bf
y})A_{\beta+\gamma n}({\bf x},{\bf
y})}{n+\mu}=\frac{\pi}{\sin{\mu\pi}}A_{\alpha+\gamma\mu}({\bf x},{\bf
y})A_{\beta-\gamma\mu}({\bf x},{\bf y})$ (16)

This equation holds for all $\mu\in\mathbb{C}\backslash\mathbb{Z}$, and all
$\alpha,\beta\in\mathbb{C}$ if $\gamma\in[0,1/2]$. We will prove this by
breaking each of the generalized Anger functions into two integrals and
passing the sum through using convergence theorems. Let us represent the first
generalized Anger function as an integral:

$\displaystyle\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}A_{\alpha-\gamma n}({\bf
x},{\bf y})A_{\beta+\gamma n}({\bf x},{\bf
y})}{n+\mu}=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\int_{-\pi}^{\pi}$
$\displaystyle\frac{(-1)^{n}e^{-i\gamma n\theta}A_{\beta+\gamma n}({\bf
x},{\bf y})}{n+\mu}\times...$
$\displaystyle...\times\exp\left[i\left(\alpha\theta-\sum_{k=1}^{m}\left(x_{k}\sin{k\theta}+y_{k}\cos{k\theta}\right)\right)\right]d\theta$
(17)

From the decay properties of the generalized Anger functions, we can use
Fubini’s theorem to interchange the summation and integral. Upon pulling the
summation inside the first integral, we then expand $A_{\beta+\gamma n}({\bf
x},{\bf y})$ into its integral representation and focus on the following
quantity:

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}e^{-i\gamma n\theta}A_{\beta+\gamma
n}({\bf x},{\bf
y})}{n+\mu}=\frac{1}{2\pi}\sum_{n=-\infty}^{\infty}\int_{-\pi}^{\pi}\frac{(-1)^{n}e^{i\gamma
n(\phi-\theta)}}{n+\mu}\exp\left[i\left(\beta\phi-\sum_{k=1}^{m}\left(x_{k}\sin{k\phi}+y_{k}\cos{k\phi}\right)\right)\right]d\phi$
(18)

We can repeat the same argument from the previous section to interchange the
summation and integral in (18).

However, when executing the complex exponential version of the sums described
by (12), we must be careful since we do not necessarily have
$\gamma(\phi-\theta)\in[-\pi,\pi]$ for general $\gamma$. For $0\leq\gamma\leq
1/2$ this does hold and we can proceed with the usual summation. By
consolidating (12) into a complex exponential form, we have

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{n+\mu}e^{in\theta}=\frac{\pi}{\sin{\pi\mu}}e^{-i\mu\theta}$
(19)

for $\theta\in(-\pi,\pi)$. Plugging this identity into the double integral
representation of the Lerche-Newberger sum and separating the factors by
integration variable gives us the result in (16). For $1/2<\gamma\leq 1$, the
quantity $\gamma(\phi-\theta)$ will extend beyond $(-\pi,\pi)$. Indeed, if
$\theta\in(\pi,3\pi)$, then $\theta-2\pi\in(-\pi,\pi)$ such that

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}}{n+\mu}e^{in\theta}=\frac{\pi}{\sin{\pi\mu}}e^{-i\mu(\theta-2\pi)}$
(20)

A similar identity holds for $\theta\in(-3\pi,-\pi)$. If we restrict
$\alpha,\beta\in\mathbb{Z}$, then we can write the double integral
representation as

$\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}A_{\alpha-\gamma n}({\bf x},{\bf
y})A_{\beta+\gamma n}({\bf x},{\bf
y})}{n+\mu}=\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}f_{\mu}(\gamma(\phi-\theta))p_{\alpha}(\theta)p_{\beta}(\phi)d\phi
d\theta$ (21)

where $f_{\mu}(\theta)$ is the summation in (19), and $p_{\alpha}(\theta)$ is
the $2\pi$-periodic integrand of $A_{\alpha}({\bf x},{\bf y})$. In Figure 1,
we decompose the region of integration according to the period that
$\gamma(\phi-\theta)$ lies in. In the triangle region bounded by
$\theta=-\pi$, $\phi=\pi$, and $\gamma(\phi-\theta)=\pi$, we have
$\gamma(\phi-\theta)\geq\pi$ so the summation formula from (19) incurs an
extra factor of $e^{2\pi\mu i}$ as shown in (20). Let this region be called
$\Delta_{+}(\gamma)$. In the negation of that region, the triangle bounded by
$\theta=\pi$, $\phi=-\pi$ and $\gamma(\phi-\theta)=-\pi$ which we call
$\Delta_{-}(\gamma)$, the summation procedure incurs an extra factor of
$e^{-2\pi\mu i}$. In this way, we can replace $f_{\mu}(\gamma(\phi-\theta))$
by a complex exponential and consolidate the expression into an Anger function
product and two integrals over $\Delta_{\pm}(\gamma)$:

$\displaystyle\sum_{n=-\infty}^{\infty}\frac{(-1)^{n}A_{\alpha-\gamma n}({\bf
x},{\bf y})A_{\beta+\gamma n}({\bf x},{\bf y})}{n+\mu}$
$\displaystyle=\frac{\pi}{\sin{\mu\pi}}A_{\alpha+\gamma\mu}({\bf x},{\bf
y})A_{\beta-\gamma\mu}({\bf x},{\bf y})+2\pi
ie^{2i\pi\mu}\iint_{\Delta_{+}}p_{\alpha+\gamma\mu}(\theta)p_{\beta-\gamma\mu}(\phi)d\theta
d\phi$ $\displaystyle-2\pi
ie^{-2i\pi\mu}\iint_{\Delta_{-}}p_{\alpha+\gamma\mu}(\theta)p_{\beta-\gamma\mu}(\phi)d\theta
d\phi$ (22)

Figure 1: Division of region of integration in (22).

##  4 Applications

Lerche-Newberger type summations appear in many different areas of physics
including plasma physics [20] and periodically driven quantum mechanical
systems [11]. In particular, the following expression appears in Berns et. al.
[12] to describe the transition rate of a persistent qubit:

$W(\epsilon,A)=\frac{\Delta^{2}}{2}\sum_{n}\frac{\Gamma_{2}J_{n}^{2}(x)}{(\epsilon-\omega
n)^{2}+\Gamma_{2}^{2}}$ (23)

where $x=A/\omega$. The authors give an asymptotic treatment of this function,
but a closed form expression is possible using (1). Using a partial fractions
expansion, we can break the denominator into its linear factors. Let
$\mu_{\pm}=(\epsilon\pm\Gamma i)/\omega$ such that

$W(\epsilon,A)=\frac{i\omega\Delta^{2}}{4}\left[\sum_{n}\frac{J_{n}^{2}(x)}{n-\mu_{+}}-\sum_{n}\frac{J_{n}^{2}(x)}{n-\mu_{-}}\right]$
(24)

By letting $\alpha=\beta=0$ and $\gamma=1$, we can apply (1) to (24) to arrive
at a closed form expression:

$W(\epsilon,A)=\frac{i\pi\omega\Delta^{2}}{4}\left[\frac{J_{\mu_{+}}(x)J_{-\mu_{+}}(x)}{\sin\mu_{+}\pi}-\frac{J_{\mu_{-}}(x)J_{-\mu_{-}}(x)}{\sin\mu_{-}\pi}\right]$
(25)

We display plots of the expression in (25) for several selections of
parameters in Figure 2. Now asymptotics need only be developed for the
constituent Bessel functions rather than for the complicated sum in (23).

Figure 2: Plots of $W(\epsilon,A)$ for a variety of parameter sets (_Blue_)
$\Gamma_{2}=3$, $\epsilon=2.1$, $\omega=0.07$ (_Red_) $\Gamma_{2}=10$,
$\epsilon=1.3$, $\omega=2.0$ (_Green_) $\Gamma_{2}=1$, $\epsilon=0.6$,
$\omega=3.0$.

If this system was modified by letting the qubit be modulated by several
frequencies possibly out of phase where instead of the energy detuning being
of the form $h(t)=\epsilon+\delta\epsilon(t)+A\cos{2\pi vt}$, we would include
additional terms such that
$h(t)=\epsilon+\delta\epsilon(t)+\sum_{k=1}^{m}\left(A_{k}\cos{2\pi
vkt}+B_{k}\cos{2\pi vkt}\right)$. Then the Bessel functions in (23) would be
replaced with generalized Bessel functions of the form $J_{n}({\bf x},{\bf
y})$ where ${\bf x}=(A_{1}/\omega,...,A_{m}/\omega)$ and ${\bf
y}=(B_{1}/\omega,...,B_{m}/\omega)$. After using another partial fractions
decomposition, a variant of the generalized Lerche-Newberger summation would
emerge:

$\sum_{n=-\infty}^{\infty}\frac{J_{n}({\bf x},{\bf y})^{2}}{n+\mu}$ (26)

Since for generalized Bessel functions we have $J_{-n}({\bf x},{\bf
y})\neq(-1)^{n}J_{n}({\bf x},{\bf y})$, we must use a variant of (19). By
representing $(-1)^{n}=e^{i\pi n}$ in (19), we can shift the piecewise
continuous regions from $((2k-1)\pi,(2k+1)\pi)$ to $(2k\pi,2(k+1)\pi)$.
Indeed, if $\theta\in(0,2\pi)$ we have

$\sum_{n=-\infty}^{\infty}\frac{e^{in\theta}}{n+\mu}=\frac{\pi}{\sin\pi\mu}e^{i\mu\pi}e^{-i\mu\theta}$
(27)

If $\theta\in(-2\pi,0)$, then a similar argument in (27) would show that

$\sum_{n=-\infty}^{\infty}\frac{e^{in\theta}}{n+\mu}=\frac{\pi}{\sin\pi\mu}e^{-i\mu\pi}e^{-i\mu\theta}$
(28)

Using the typical integral interchange arguments, we can express (26) as the
double integral:

$\sum_{n=-\infty}^{\infty}\frac{J_{n}({\bf x},{\bf
y})^{2}}{n+\mu}=\frac{1}{4\pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}g_{\mu}(\phi+\theta)p_{0}(\theta)p_{0}(\phi)d\theta
d\phi$ (29)

where $g_{\mu}(\theta)$ is the summation from (27) and (28) and
$p_{0}(\theta)$ is the same as in the previous section. In Figure 3, we see
how the discontinuity of $g_{\mu}(\theta)$ divides the region of integration
into two subregions on either side of the curve $\phi+\theta=0$ which we refer
to as $\tilde{\Delta}_{\pm}$. Therefore we can divide this into the following
sum:

$\sum_{n=-\infty}^{\infty}\frac{J_{n}({\bf x},{\bf
y})^{2}}{n+\mu}=\frac{1}{4\pi\sin{\pi\mu}}\left[e^{i\mu\pi}\iint_{\tilde{\Delta}_{+}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi+e^{-i\mu\pi}\iint_{\tilde{\Delta}_{-}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi\right]$ (30)

Figure 3: Division of region of integration in (30).

By splitting the complex exponentials into their trigonometric components, we
find that this expression resolves based on the sum and difference of the two
integral terms. Since
$\tilde{\Delta}_{+}\cup\tilde{\Delta}_{-}=[-\pi,\pi]\times[-\pi,\pi]$, it
follows that

$\iint_{\tilde{\Delta}_{+}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi+\iint_{\tilde{\Delta}_{-}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi=4\pi^{2}A_{-\mu}({\bf x},{\bf y})^{2}$ (31)

Meanwhile, the difference of these two integrals can be combined into a single
integral over the space by using the sign function $\text{sgn}(x)$ such that

$\iint_{\tilde{\Delta}_{+}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi-\iint_{\tilde{\Delta}_{-}}p_{-\mu}(\theta)p_{-\mu}(\phi)d\theta
d\phi=\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}p_{-\mu}(\theta)p_{-\mu}(\phi)\text{sgn}(\theta+\phi)d\theta
d\phi$ (32)

and we represent this integral as $4\pi^{2}B_{-\mu}({\bf x},{\bf y})$ to
compare with (31). Substituting (31) and (32) into (30) we arrive at the final
expression:

$\sum_{n=-\infty}^{\infty}\frac{J_{n}({\bf x},{\bf
y})^{2}}{n+\mu}=2\pi\left(\cot(\pi\mu)A_{-\mu}({\bf x},{\bf
y})^{2}+iB_{-\mu}({\bf x},{\bf y})\right)$ (33)

The asymptotics of $A_{-\mu}({\bf x},{\bf y})$ have been explored for ${\bf
y}=0$ in several sources [14] [18] and in more generality in Kuklinski and
Hague [13]. The component $B_{-\mu}({\bf x},{\bf y})$ is a two-dimensional
oscillatory integral over the region $[-\pi,\pi]\times[-\pi,\pi]$, however the
standard results from Stein [24] do not apply since the integrand has a
discontinuity over the diagonal of this region contributed by
$\text{sgn}(\phi+\theta)$ as seen in Figure 3. Regardless of the contributions
from the discontinuity that may persist for all parameter choices $\mu,{\bf
x},{\bf y}$, we still see that $B_{-\mu}({\bf x},{\bf y})$ has the same
bifurcation curves as $A_{-\mu}({\bf x},{\bf y})$ due to the similar structure
of the integrand. We plot examples of these sums in Figure 4.

Figure 4: Plots of sum (33). (_Left_) ${\bf x}=(x,y,0,...)$, ${\bf y}=0$,
$\mu=0.5$ (_Center_) ${\bf x}=(x,0,y,0,...)$, ${\bf y}=0$, $\mu=0.3$ (_Right_)
${\bf x}=(0,x,y,0,...)$, ${\bf y}=0$, $\mu=1.7$.

##  5 Conclusion

In this document we extended the Lerche-Newberger formula to incorporate
higher dimensional analogues to the Bessel functions. Since these generalized
Bessel functions are only well defined for integer order, we instead opted to
use generalize Anger functions in the summation. While the one-dimensional
Lerche-Newberger formula generally preserves the Bessel function structure,
the Anger function analogue does not behave as nicely since discontinuities in
the equivalent integral expressions are introduced. Even for a relatively
simple multi-dimensional as presented in the previous section, double
integrals over discontinuous functions emerge that ostensibly cannot be
reduced to well known functions using symmetry relations.

It is the hope of the authors that these results will contribute to future
investigations of oscillatory quantum mechanical systems.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] I. Lerche and R. Tautz, “Kapteyn series arising in radiation problems,” _Journal of physics A: Mathematical and theoretical_ , vol. 41, no. 3, p. 035202, 2008.
  * [2] S. Szapiel, “Maréchal intensity criteria modified for circular apertures with nonuniform intensity transmission: Dini series approach,” _Optics letters_ , vol. 2, no. 5, pp. 124–126, 1978.
  * [3] G. Watson, _A treatise on the theory of Bessel functions_. Cambridge university press, 1922.
  * [4] M. Abramowitz and I. A. Stegun, _Handbook of mathematical functions with formulas, graphs, and mathematical tables_. US Government printing office, 1964, vol. 55.
  * [5] I. Lerche and R. C. Tautz, “A note on summation of Kapteyn series in astrophysical problems,” _The Astrophysical Journal_ , vol. 665, no. 2, p. 1288, 2007.
  * [6] V. V. Kravchenko and S. M. Torba, “Asymptotics with respect to the spectral parameter and neumann series of bessel functions for solutions of the one-dimensional schrödinger equation,” _Journal of Mathematical Physics_ , vol. 58, no. 12, p. 122107, 2017.
  * [7] C. Linton, “Schlömilch series that arise in diffraction theory and their efficient computation,” _Journal of Physics A: Mathematical and General_ , vol. 39, no. 13, p. 3325, 2006.
  * [8] B. S. Newberger, “New sum rule for products of bessel functions with application to plasma physics,” _Journal of Mathematical Physics_ , vol. 23, no. 7, pp. 1278–1281, 1982.
  * [9] O. Buneman, “Dissipation of currents in ionized media,” _Physical Review_ , vol. 115, no. 3, p. 503, 1959.
  * [10] I. Lerche, R. Schlickeiser, and R. Tautz, “Comment on “a new derivation of the plasma susceptibility tensor for a hot magnetized plasma without infinite sums of products of bessel functions”[phys. plasmas 14, 092103 (2007)],” _Physics of Plasmas_ , vol. 15, no. 2, p. 092103, 2008.
  * [11] M. Silveri, J. Tuorila, E. Thuneberg, and G. Paraoanu, “Quantum systems under frequency modulation,” _Reports on Progress in Physics_ , vol. 80, no. 5, p. 056002, 2017.
  * [12] D. Berns, W. Oliver, S. Valenzuela, A. Shytov, K. Berggren, L. Levitov, and T. Orlando, “Coherent quasiclassical dynamics of a persistent current qubit,” _Physical review letters_ , vol. 97, no. 15, p. 150502, 2006.
  * [13] P. Kuklinski and D. A. Hague, “Identities and properties of multi-dimensional generalized bessel functions,” _arXiv preprint arXiv:1908.11683_ , 2019.
  * [14] G. Dattoli and A. Torre, _Theory and applications of generalized Bessel functions_. Arcane, 1996.
  * [15] D. A. Hague, “Adaptive transmit waveform design using multitone sinusoidal frequency modulation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 57, no. 2, pp. 1274–1287, 2021.
  * [16] H. R. Reiss, “Absorption of light by light,” _Journal of Mathematical Physics_ , vol. 3, no. 1, pp. 59–67, 1962. [Online]. Available: https://doi.org/10.1063/1.1703787
  * [17] W. Paciorek and G. Chapuis, “Generalized bessel functions in incommensurate structure analysis,” _Acta Crystallographica Section A: Foundations of Crystallography_ , vol. 50, no. 2, pp. 194–203, 1994.
  * [18] H. Korsch, A. Klumpp, and D. Witthaut, “On two-dimensional bessel functions,” _Journal of Physics A: Mathematical and General_ , vol. 39, no. 48, 2006.
  * [19] M. Bakker and N. M. Temme, “Sum rule for products of bessel functions: Comments on a paper by newberger,” _Journal of mathematical physics_ , vol. 25, no. 5, pp. 1266–1267, 1984.
  * [20] I. Lerche, “A note on summing series of bessel functions occurring in certain plasma astrophysical situations,” _The Astrophysical Journal_ , vol. 190, pp. 165–166, 1974.
  * [21] H. L. Royden and P. Fitzpatrick, _Real analysis_. Macmillan New York, 1988, vol. 32.
  * [22] T. M. Apostol _et al._ , “On the lerch zeta function.” _Pacific Journal of Mathematics_ , vol. 1, no. 2, pp. 161–167, 1951.
  * [23] H. Wilbraham, “On a certain periodic function,” _Cambridge and Dublin Mathematical Journal_ , vol. 3, pp. 198–201, 1848.
  * [24] E. M. Stein, _Harmonic Analysis: Real-Variable Methods, Orthogonality, and Oscillatory Integrals_. Princeton University Press, 1993.
  * [25] E. Hewitt and R. E. Hewitt, “The gibbs-wilbraham phenomenon: an episode in fourier analysis,” _Archive for history of Exact Sciences_ , vol. 21, no. 2, pp. 129–160, 1979.

##  Appendix A Proof of Global Maximum Location

In this section we prove that there exists a constant $M\in\mathbb{N}$ such
that for all $n>M$, (7) has a global maximum at $\theta_{n}=\pi-\pi/(2n+1)$.
We prove this for even $n$ so we replace $n$ with $2n$, but the result can
easily be extended to $n$ odd. To do this we replicate an argument from Hewitt
and Hewitt [25]. Let $f(\theta)$ be the difference from (8) which we write
here:

$f(\theta)=\sum_{k=1}^{2n}\frac{(-1)^{k}}{k}\cos{k\theta}+\log\left(2\cos\frac{\theta}{2}\right)$
(34)

Using trigonometric identities, we can prove that the derivative of this
function satisfies the following:

$f^{\prime}(\theta)=-\sin\left(\frac{4n+1}{2}\theta\right)/\left(2\cos\left(\frac{\theta}{2}\right)\right)$
(35)

This gives us that the local extrema of $f$ satisfy $\theta_{k}=\frac{2\pi
k}{4n+1}$, and by conducting a second derivative test we can conclude that
$\theta_{2k}$ are the local maxima of $f$. We first will prove that
$f(\theta_{2k+2})-f(\theta_{2k})>f(\theta_{2k})-f(\theta_{2k-2})$. To do this
let us define a function $\omega$:

$\displaystyle\omega(t)$
$\displaystyle=\left[f\left(\frac{4\pi(k+1/2)}{4n+1}+\frac{2t}{4n+1}\right)-f\left(\frac{4\pi(k-1/2)}{4n+1}+\frac{2t}{4n+1}\right)\right]$
$\displaystyle-\left[f\left(\frac{4\pi(k+1/2)}{4n+1}-\frac{2t}{4n+1}\right)-f\left(\frac{4\pi(k-1/2)}{4n+1}-\frac{2t}{4n+1}\right)\right]$
(36)

Here, we restrict $1\leq k\leq n-1$. This function is of interest since if
$\omega(\pi)>0$, then the inequality in question holds. Since $\omega(0)=0$,
we need only show that $\omega^{\prime}(t)>0$ for $0\leq t\leq\pi$. Indeed,
$\omega^{\prime}(t)$ takes the following form:

$\displaystyle\omega^{\prime}(t)$ $\displaystyle=\frac{2\sin
t}{4n+1}\left[\left(\frac{1}{\cos\left(\frac{2\pi(k+1/2)}{4n+1}+\frac{t}{4n+1}\right)}-\frac{1}{\cos\left(\frac{2\pi(k-1/2)}{4n+1}+\frac{t}{4n+1}\right)}\right)\right.$
$\displaystyle-\left.\left(\frac{1}{\cos\left(\frac{2\pi(k+1/2)}{4n+1}-\frac{t}{4n+1}\right)}-\frac{1}{\cos\left(\frac{2\pi(k-1/2)}{4n+1}-\frac{t}{4n+1}\right)}\right)\right]$
(37)

Using the identity
$1/\cos(a+b)-1/\cos(a-b)=4(\sin{a}\sin{b})/(\cos{2a}+\cos{2b})$, we can
rewrite the function as follows:

$\omega^{\prime}(t)=\frac{8\sin
t\sin\left(\frac{\pi}{4n+1}\right)}{4n+1}\left[\frac{\sin\left(\frac{2\pi
k+t}{4n+1}\right)}{\cos\left(\frac{4\pi
k+2t}{4n+1}\right)+\cos\left(\frac{2\pi}{4n+1}\right)}-\frac{\sin\left(\frac{2\pi
k-t}{4n+1}\right)}{\cos\left(\frac{4\pi
k-2t}{4n+1}\right)+\cos\left(\frac{2\pi}{4n+1}\right)}\right]$ (38)

After combining these two fractions and conducting some elementary
trigonometric manipulations, we can further simplify this derivative:

$\omega^{\prime}(t)=\frac{16\sin
t\sin\left(\frac{\pi}{4n+1}\right)\sin\left(\frac{t}{4n+1}\right)\cos\left(\frac{2\pi
k}{4n+1}\right)}{(4n+1)\left(\cos\left(\frac{4\pi
k+2t}{4n+1}\right)+\cos\left(\frac{2\pi}{4n+1}\right)\right)}$ (39)

Due to the restrictions on $t$ and $k$, all of the terms in this fraction are
positive and therefore the desired condition holds, namely that the difference
between adjacent local maxima increases to the right.

We ultimately want to show that $\theta_{2n}$ is the global maxima. This can
be accomplished by proving $f(\theta_{2})>f(\theta_{0})$ and using induction
with the above result. Though this result appears true for all $n$, due to the
unwieldy nature of $f(\theta)$ we instead opt to prove this asymptotically,
that there exists some $M$ such that for all $n>M$ the inequality is
satisfied. We do this by considering the asymptotic expansion of
$f(\theta_{2})-f(\theta_{0})$:

$f(\theta_{2})-f(\theta_{0})=-2\sum_{k=1}^{2n}\frac{(-1)^{k}}{k}\sin^{2}\left(\frac{2\pi
k}{4n+1}\right)+\log\left(\cos\frac{2\pi}{4n+1}\right)$ (40)

If we can show that the leading term of the asymptotic expansion of (40) in
$n$ is positive, then the proof follows. The logarithm term can be easily
expanded using typical Taylor expansion arguments:

$\log\left(\cos\frac{2\pi}{4n+1}\right)=\frac{1}{n^{2}}\left(-\frac{\pi^{2}}{8}\right)+\frac{1}{n^{3}}\left(\frac{\pi^{2}}{16}\right)+O(n^{-4})$
(41)

To resolve the summation term, we manipulate it into an endpoint Riemann sum.
Recall that for a smooth function $f(x)$, the following asymptotic form holds
[wals37]:

$\frac{1}{n}\sum_{k=1}^{n}f\left(\frac{k}{n}\right)=\int_{0}^{1}f(x)dx+\frac{1}{2n}\int_{0}^{1}f^{\prime}(x)dx+\frac{1}{12n^{2}}\int_{0}^{1}f^{\prime\prime}(x)dx+O(n^{-4})$
(42)

We split the summation into positive and negative terms according to the
$(-1)^{k}$ term:

$\sum_{k=1}^{2n}\frac{(-1)^{k}}{k}\sin^{2}\left(\frac{2\pi
k}{4n+1}\right)=\sum_{k=1}^{n}\frac{1}{2k}\sin^{2}\left(\frac{4\pi
k}{4n+1}\right)-\sum_{k=1}^{n}\frac{1}{2k-1}\sin^{2}\left(\frac{2\pi(2k-1)}{4n+1}\right)$
(43)

We elaborate the asymptotic expansion of the even sum on the right hand side
of (43); the odd sum is similar. By letting $f(x)=\sin^{2}(\pi x)/x$, we
rewrite this sum as

$\sum_{k=1}^{n}\frac{1}{2k}\sin^{2}\left(\frac{4\pi
k}{4n+1}\right)=\frac{2}{4n+1}\sum_{k=1}^{n}f\left(\frac{4k}{4n+1}\right)$
(44)

The remainder of this proof is an arduous calculation of these sum up to order
for large $n$ which we outline here. First, we expand the argument of $f$,
$4k/(4n+1)$, being careful to recognize that $k$ is at the same order as $n$
so that $k/n=O(1)$. Next, we conduct a Taylor expansion of $f$ at $k/n$
throwing out higher order terms such that we are left with terms that look
like $(k/n)^{j}f^{(j)}(k/n)$. By considering these the summands of a Riemann
endpoint sum of function $g(x)=x^{j}f^{(j)}$, we use (42) to further expand
these sums into a collection of integrals. These integrals have closed form
expressions, and we can finish the problem off by multiplying the expansion of
the sum on the right hand side of (44) by the expansion of the factor
$2/(4n+1)$. This gives us the asymptotic expansions:

$\sum_{k=1}^{n}\frac{1}{2k}\sin^{2}\left(\frac{4\pi k}{4n+1}\right)=\frac{\pi
x}{2}-\frac{1}{n^{2}}\left(\frac{\pi^{2}}{24}\right)+\frac{1}{n^{3}}\left(\frac{5\pi^{2}}{384}\right)+O(n^{-4})$
(45)
$\sum_{k=1}^{n}\frac{1}{2k-1}\sin^{2}\left(\frac{2\pi(2k-1)}{4n+1}\right)=\frac{\pi
x}{2}+\frac{1}{n^{2}}\left(\frac{\pi^{2}}{48}\right)-\frac{1}{n^{3}}\left(\frac{\pi^{2}}{384}\right)+O(n^{-4})$
(46)

Here, $x=(-\text{Ci}(2\pi)+\gamma+\log(2\pi))/2\pi$. Combining the expansions
from (45) and (46) with the expansion of the logarithm in (41), we arrive at
the final result, namely

$f(\theta_{2})-f(\theta_{0})=\frac{\pi^{2}}{32n^{3}}+O(n^{-4})$ (47)

and since the leading quantity is positive, the rest of the proof follows.


11institutetext: 1 Sano Centre for Computational Medicine, Cracow, Poland
2 Technical University of Cluj-Napoca, Cluj-Napoca, Romania
3 The University of Sheffield, Sheffield, United Kingdom
4 Insigneo Institute for in Silico Medicine, University of Sheffield,
Sheffield, UK
5 Warsaw University of Technology 6 Tooploox 7 Jagiellonian University of
Cracow

# Noninvasive Estimation of Mean Pulmonary Artery Pressure Using MRI, Computer
Models, and Machine Learning

Michal K. Grzeszczyk 11 0000-0002-5304-1020 Tadeusz Satława 11 Angela Lungu 22
0000-0002-4531-2791 Andrew Swift 33 Andrew Narracott 3344 0000-0002-3068-6192
Rod Hose 33 Tomasz Trzcinski 556677 0000-0002-1486-8906 Arkadiusz Sitek 11
0000-0002-0677-4002

###### Abstract

Pulmonary Hypertension (PH) is a severe disease characterized by an elevated
pulmonary artery pressure. The gold standard for PH diagnosis is measurement
of mean Pulmonary Artery Pressure (mPAP) during an invasive Right Heart
Catheterization. In this paper, we investigate noninvasive approach to PH
detection utilizing Magnetic Resonance Imaging, Computer Models and Machine
Learning. We show using the ablation study, that physics-informed feature
engineering based on models of blood circulation increases the performance of
Gradient Boosting Decision Trees-based algorithms for classification of PH and
regression of values of mPAP. We compare results of regression (with
thresholding of estimated mPAP) and classification and demonstrate that
metrics achieved in both experiments are comparable. The predicted mPAP values
are more informative to the physicians than the probability of PH returned by
classification models. They provide the intuitive explanation of the outcome
of the machine learning model (clinicians are accustomed to the mPAP metric,
contrary to the PH probability).

###### Keywords:

Pulmonary Hypertension Regression Gradient Boosting Decision Trees
Mathematical Modelling

##  1 Introduction

Pulmonary Hypertension is a severe disease difficult to diagnose with multiple
possible root causes [6]. For many years, PH was identified if a mean
Pulmonary Artery Pressure (mPAP) of a patient at rest was equal to or above 25
mmHg. Recently, it has been suggested to lower the threshold to 20 mmHg [19].
The precise measurement of mPAP is non-trivial and requires conducting an
invasive Right Heart Catheterization (RHC) – the gold standard for diagnosing
PH. This procedure carries risks, requires patient’s preparation, trained
staff, highly specialized equipment, it is expensive and time consuming. To
lower the probability of complications it has to be performed at a specialized
facility [5].

Non-invasive estimation of mPAP using medical imaging, mathematical modeling,
and machine learning (ML) is an option to avoid issues related with RHC.
Mathematical models, such as a Windkessel model, allow diagnosis of the
vascular system parameters [23]. Different ML algorithms enable extracting
knowledge about data samples and their performance usually increases with the
addition of features from multiple domains.

In this paper, we present methods based on Gradient Boosting Decision Trees
(GBDT) for non-invasive PH diagnosis. We use classic GBDT, DART (Dropouts meet
Multiple Additive Regression Trees) [22] \- a method utilizing dropouts of
random trees during training - and GOSS (Gradient-based One-Side Sampling)
[10] – a technique that uses different than GBDT process of training
(retaining samples with large gradients and randomly dropping the ones with
low gradients). We conduct analysis on data from 352-patient cohort and
perform two tasks: classification of PH and regression of mPAP. As predictors,
we use demographics features, measurements derived from Magnetic Resonance
Imaging (MRI) and features obtained from 0D and 1D mathematical models [15].

Our main contribution is the demonstration of the ablation study, which shows,
that physics-informed feature engineering based on mathematical models of
blood circulation increases the performance of ML algorithms for
classification and regression of PH and values of mPAP, respectively. Another
significant contribution of this paper is comparison of utilities of
classification and regression approaches for the detection of PH. While the
regression achieves similar classification metrics (after thresholding of
estimated mPAP), the values of predicted mPAP are more informative to the
physicians than the probability of PH returned by classification models. As
such, they provide the intuitive explanation of the outcome of the machine
learning model (clinicians are accustomed to the mPAP metric, contrary to the
PH probability).

##  2 Related work

Multiple ML algorithms (utilizing features from various modalities like
echocardiography, Computed Tomography (CT), or MRI) have been integrated for
the purpose of the PH classification. In [14], five ML models were used and
compared with each other. Boosted Classification Trees, Lasso Penalized
Logistic Regression (LPLR), Random Forest (RF) for Regression, RF for
Classification and Support Vector Machines (SVM) were adopted for mPAP
prediction or PH classification basing on the echocardiographic measurements
and basic patients characteristics (age, sex, BMI, body surface area). In
[26], echocardiographic data was used to distinguish between pre- and post-
capillary PH with one of the nine tested ML models (SVM, AdaBoost, LR, RF,
Decision Trees (DT), K-Nearest Neighbours, GBDT, LogitBoost and Linear
Discriminant Analysis (LDA)).

In [7], measurements derived from CT were used to train six ML classifiers to
evaluate the probability of mPAP higher than 15 mmHg. Another approach was to
record the heart sounds with a digital stethoscope to gather parameters for PH
classification using LDA [2]. The analysis of the sounds revealed specific
patterns in PH patients. In [1], it was noted that the sounds collected by
phonocardiogram can be applied for binary classification of PH with SVM. In
[16], it was shown that MRI measurements combined with parameters from 0D and
1D computational models can be successfully used for PH and non-PH patients
classification with DT. In our approach, we study the impact of mathematical
models parameters on classification and regression. We also show the
comparable performance of PH diagnosis with GBDT-based models in both tasks.

With the rise of Deep Learning (DL), multiple approaches of detecting PH
directly from images, videos, or electrocardiography (ECG) signals were
investigated. For example, chest X-Ray images can be utilized for binary
classification of potential PH patients using Capsule Network with residual
blocks [12]. In [27], three popular DL networks (ResNet50, Xception and
Inception V3) were trained as predictors of PH. As shown in [13], an ensemble
neural network can pose as a screening tool for PH from a 12-lead ECG signal.

ML can also be utilized for determining patients at risk of having Pulmonary
Arterial Hypertension (PAH) from clinical records. In [11], it was shown that
GBDT can help in screening for PAH based on their medical history. ML-based
tools were also developed for the purpose of blood pressure estimation - in
[25], Support Vector Machine Regression (SVR) models were applied for the
prediction of the patient’s blood pressure from the physiological data.
Another example is an application of Multilayer Perceptron (MLP) for
regression of systolic blood pressure using basic knowledge about patients
(BMI, age, habits etc.) [24].

##  3 Methods

In this section, we describe our approaches to noninvasive PH diagnosis. We
present the details of our dataset and introduce mathematical models which
enabled the acquisition of physics-informed features. Finally, we train GBDT-
based models on multiple feature sets to perform mPAP regression and PH
classification experiments.

###  3.1 PH dataset

Table 1: PH dataset with patient related data, parameters derived from 0D and 1D models and measurements from MRI imaging. In the appendix (section 8) we provide explanations for the feature names. P-value tests a null hypothesis that the coefficient of the univariate linear regression between a feature and mPAP is equal to zero. | No PH | PH |
---|---|---|---
feature | cnt | mean | std | cnt | mean | std | p-value
mPAP, mmHg | 66 | 19.67 | 3.34 | 286 | 46.95 | 13.08 |
Demographics |  |  |  |  |  |  |
age, years | 66 | 56.61 | 13.78 | 286 | 61.69 | 14.24 | 0.242
gender, female/male | 66 | 43/23 |  | 286 | 173/113 |  | 0.549
who, no. | 56 | 2.52 | 0.54 | 285 | 3.04 | 0.44 | $<0.001$
bsa, $m^{2}$ | 65 | 1.88 | 0.25 | 286 | 1.82 | 0.22 | 0.24
0D and 1D models |  |  |  |  |  |  |
Rd, $kg/m^{4}s$ | 66 | 6.08E+07 | 4.94E+07 | 286 | 1.46E+08 | 2.53E+08 | $<0.001$
Rc, $kg/m^{4}s$ | 66 | 7.94E+06 | 7.80E+06 | 286 | 9.17E+06 | 1.87E+07 | 0.072
C, $m^{4}s^{2}/kg$ | 66 | 9.92E-09 | 6.71E-09 | 284 | 3.94E-04 | 6.65E-03 | 0.669
Rtot, $kg/m^{4}s$ | 66 | 6.83E+07 | 5.38E+07 | 286 | 1.56E+08 | 2.62E+08 | $<0.001$
Wb/Wtot | 66 | 0.24 | 0.10 | 286 | 0.39 | 0.11 | $<0.001$
MRI |  |  |  |  |  |  |
rac_fiesta, % | 66 | 26.39 | 15.43 | 286 | 13.68 | 8.93 | $<0.001$
syst_area_fiesta, $cm^{2}$ | 66 | 7.62 | 2.17 | 286 | 9.78 | 2.78 | $<0.001$
diast_area_fiesta, $cm^{2}$ | 66 | 6.08 | 1.71 | 286 | 8.66 | 2.57 | $<0.001$
rvedv, $mL$ | 66 | 118.93 | 36.00 | 286 | 159.58 | 58.27 | $<0.001$
rvedv_index, $mL/m^{2}$ | 66 | 53.78 | 21.83 | 286 | 73.92 | 39.39 | $<0.001$
rvesv, $mL$ | 66 | 55.41 | 20.68 | 286 | 102.48 | 49.92 | $<0.001$
rvesv_index, $mL/m^{2}$ | 66 | 24.64 | 10.84 | 286 | 47.63 | 30.19 | $<0.001$
rvef, % | 66 | 53.32 | 9.86 | 286 | 38.05 | 13.59 | $<0.001$
rvsv, $mL$ | 66 | 63.52 | 22.61 | 286 | 57.15 | 23.39 | 0.026
rvsv_index, $mL/m^{2}$ | 66 | 29.14 | 13.90 | 286 | 26.32 | 15.02 | 0.292
lvedv, $mL$ | 66 | 116.57 | 33.09 | 286 | 91.30 | 27.33 | $<0.001$
lvedv_index, $mL/m^{2}$ | 66 | 53.16 | 21.90 | 286 | 41.25 | 19.20 | $<0.001$
lvesv, $mL$ | 66 | 34.27 | 15.66 | 286 | 31.32 | 14.56 | 0.23
lvesv_index, $mL/m^{2}$ | 66 | 16.85 | 16.81 | 286 | 14.01 | 8.18 | 0.194
lvef, % | 66 | 71.13 | 8.54 | 286 | 65.81 | 10.92 | $<0.001$
lvsv, $mL$ | 66 | 82.30 | 23.30 | 286 | 59.97 | 19.93 | $<0.001$
lvsv_index, $mL/m^{2}$ | 66 | 38.07 | 16.20 | 286 | 27.20 | 13.51 | $<0.001$
rv_dia_mass, $g$ | 66 | 22.62 | 6.80 | 283 | 44.48 | 25.47 | $<0.001$
lv_dia_mass, $g$ | 66 | 91.47 | 27.71 | 286 | 90.64 | 24.98 | 0.436
lv_syst_mass, $g$ | 66 | 111.74 | 32.17 | 286 | 99.83 | 26.39 | $<0.001$
rv_mass_index, $g/m^{2}$ | 66 | 10.44 | 4.94 | 285 | 20.94 | 15.09 | $<0.001$
lv_mass_index, $g/m^{2}$ | 59 | 40.90 | 17.87 | 243 | 39.84 | 18.99 | 0.442
sept_angle_syst, degrees | 66 | 139.95 | 11.68 | 286 | 172.51 | 22.11 | $<0.001$
sept_angle_diast, degrees | 66 | 134.21 | 8.28 | 286 | 145.01 | 11.93 | $<0.001$
4ch_la_area, $mm^{2}$ | 66 | 1921.95 | 387.56 | 286 | 1785.95 | 556.53 | $<0.001$
4ch_la_length, $mm^{2}$ | 66 | 55.76 | 7.86 | 286 | 55.62 | 8.60 | 0.412
2ch_la_area, $mm^{2}$ | 66 | 1764.62 | 496.75 | 286 | 1901.67 | 545.35 | 0.855
2ch_la_length, $mm^{2}$ | 66 | 48.66 | 9.08 | 286 | 52.12 | 9.33 | 0.166
la_volume, $mL$ | 66 | 55.22 | 17.96 | 286 | 54.16 | 25.36 | 0.005
la_volume_index, $mL/m^{2}$ | 66 | 24.95 | 10.14 | 286 | 23.24 | 10.45 | 0.042
ao_qflowpos, $L/min$ | 65 | 6.09 | 1.50 | 285 | 5.29 | 1.50 | $<0.001$
ao_qfp_ind, $L/min/m^{2}$ | 65 | 2.79 | 1.18 | 285 | 2.44 | 1.15 | 0.003
pa_qflowpos, $L/min$ | 66 | 5.50 | 1.84 | 284 | 5.00 | 1.97 | 0.006
pa_qflowneg, $L/min$ | 66 | 0.62 | 0.59 | 285 | 1.07 | 0.83 | $<0.001$
pa_qfn_ind, $L/min/m^{2}$ | 66 | 9.70 | 7.19 | 284 | 17.49 | 9.85 | $<0.001$
systolic_area_pc, $mm^{2}$ | 66 | 731.05 | 236.42 | 284 | 950.17 | 268.98 | $<0.001$
diastolic_area_pc, $mm^{2}$ | 66 | 619.82 | 162.71 | 284 | 866.42 | 244.57 | $<0.001$
rac_pc, % | 66 | 17.02 | 13.70 | 284 | 10.01 | 8.14 | $<0.001$

Table 1 presents the available features of patients who were suspected with PH
and underwent MRI and RHC within 48 hours.

The medical procedures were performed at the Sheffield Pulmonary Vascular
Disease Unit. The RHC procedure was conducted with a balloon-tipped 7.5-Fr
thermodilution catheter. The PH was defined if measured mPAP $\geq$ 25 mmHg.
Using these criteria out of the cohort of 352 patients 286 were diagnosed with
PH. From 286 patients with PH, 142 had Pulmonary Arterial Hypertension, 86 had
Chronic Thromboembolic PH, 35 PH cases were due to lung diseases (e.g. Chronic
Obstructive Pulmonary Disease), 15 cases were associated with left heart
disease. The cause of PH in the rest of patients was either multifactorial or
unknown. All of the available data samples are part of the ASPIRE Registry
(Assessing the Severity of Pulmonary Hypertension In a Pulmonary Hypertension
REferral Centre) [8].

MRI images were captured with 1.5-tesla whole-body scanner (GE HDx, GE
Healthcare, Milwaukee) with an 8-channel cardiac coil. The images were
acquired in the supine position during a breath hold. The balanced steady
state free precession (bSSFP) sequences were spatially and temporally
synchronized with the 2D phase contrast (PC) images of the Main Pulmonary
Artery (MPA) using cardiac gating. Short-axis and four-chamber cardiac images
were also collected. The features from MRI were obtained as in [21]. A(t) area
of the MPA was extracted from the semi-automatically segmented bSSFP images.
The blood flow through MPA (Q(t)) was extracted from the segmented areas
overlaid on PC images. Using those measurements 0D- and 1D-model features were
derived.

To prepare the feature dataset for the training of ML models we fill the
missing values using linear interpolation. We encode categorical features to
numerical values and scale all the features to have means of 0 and variances
of 1.

###  3.2 Features derived from models of blood circulation

The cardiovascular system (CVS) is a closed circuit with the main purpose of
transporting oxygenated blood to organs and tissues [17]. It comprises
especially from heart, blood and vessels. One of the main components of the
CVS is the pulmonary circulation. The target of the pulmonary circulation is
to transport the deoxygenated blood from the right ventricle through MPA and
other arteries to lungs and deliver the oxygenated blood to the left ventricle
[9]. Since CVS can be described by its haemodynamics and structure of heart
and vessels, the computational models based on the simplified representation
of CVS were introduced [18]. Those models range from 0D models simulating the
global haemodynamics (e.g. resistance and compliance of the system) to 3D
models representing the complex behaviour of vessels and the blood flow over
time. In [15], two models (0D and 1D) based on MRI measurements for the
diagnosis of PH were proposed.

####  3.2.1 0D model.

0D models are often based on the hydraulic-electrical analogue - the blood
flow and electrical circuits have many computational similarities [18]. For
example, the friction in the vessel can be identified as resistance R, the
blood pressure as voltage and the flow-rate as current. Thus, by applying
electrical laws (e.g. Kirchhoff’s law, Ohm’s law), the simplified
representation of the CVS can be achieved. 0D modelling of CVS started with
the implementation of the two-element Windkessel model [23]. Different
variants of this model appeared in the literature and it was applied to
simulate pulmonary circulation [4, 20]. The 3-element (RcCRd) Windkessel model
comprises of the capacitor C characterizing the compliance of the pulmonary
circulation and two resistors Rc and Rd representing the resistance proximal
and distal to the capacitor respectively.

In [15], RcCRd model was applied to capture the characteristics of PH and non-
PH patients. In this model, the sum of two resistors can be interpreted as the
ratio between mean pressure and mean flow (pulmonary vascular resistance -
PVR) and C indicates the compliance of the pulmonary arteries. To optimize the
parameters of 0D model for the specific patient, two MRI imaging techniques of
MPA were used: PC and bSSFP. The bSSFP images were segmented to find the area
of MPA (A(t)) over time. Then, the segmented regions were overlayed over PC
images to capture the blood flow through MPA (Q(t)). Having the Q(t) and
pressure p(t) (derived from the measured MPA radius) the parameters of the
Windkessel model which were best describing the relationship between Q(t) and
p(t) over time could be derived.

####  3.2.2 1D model.

The simplified representation of the pulmonary vasculature is multiple elastic
tubes with numerous branches. 1D models often analyse the propagation of the
pressure and flow waves in such structures. The 1D equations of the waves
travelling through elastic tubes are derived from Navier-Stokes equations. In
[15], the analysis of the power of the pressure waves was performed. The
pressure wave was broken down into forward and backward-travelling elements
(since vessels are rugged and twisted, some waves are bouncing off the vessel
walls and travel backward). It was assumed and confirmed that the power of the
backward wave in relation to the total wave power was greatly higher in PH
cases than in healthy ones. As diseased pulmonary vasculature contains more
deposits and stenoses the ratio of the backward wave power to the total wave
power (represented as Wb/Wtot) is higher than in the healthy one.

###  3.3 Machine Learning for PH detection

####  3.3.1 mPAP regression.

The decision whether the patient is suffering from PH is more important to the
doctors than the actual value of mPAP. However, the non-invasive prediction of
the PH occurrence together with the predicted value of mPAP is more
informative to the clinicians. Therefore, we decide to conduct two
experiments: mPAP regression and PH classification.

To find the best ML algorithm for mPAP regression we train three models based
on GBDT: classic GBDT, DART and GOSS. We use mPAP feature as the ground truth
for our models. We find the best hyper parameters for the models using Bayes
optimization with 8-fold cross-validation (CV). We optimize them for 200
iterations with minimizing Mean Squared Error (MSE) as the optimization
target. Then, using the best found parameters we train the models with leave
one out cross validation (LOOCV) and MSE as the objective function. We measure
MSE, Root MSE (RMSE) and Mean Absolute Error (MAE) as regression metrics of
the model.

We assume that mPAP $\geq$ 25 mmHg is a positive PH diagnosis. With this
assumption, we compute the binary classification metrics after thresholding
the predicted and measured mPAP with 25 value. We calculate accuracy,
sensitivity, specificity, True Positives (TP), False Positives (FP), True
Negatives (TN), False Negatives (FN). To compare the impact of different
feature sets on the results (demographics, MRI, mathematical models), we
repeat the procedure of hyper parameter optimization, models training with
LOOCV and metrics collection for different combinations of features. We
compare results of all the approaches.

Additionally, we train four other than boosted tree ML models on all features
and compare the metrics with GBDT-based methods using LOOCV-derived metrics.
These additional methods are MLP, SVR, AdaBoost and RF.

####  3.3.2 PH classification.

We conduct the binary PH classification, similarly to mPAP regression. We
binarize mPAP feature with 25 mmHg threshold and train three GBDT-based models
on different variations of feature sets, previously optimizing the hyper
parameters using Bayes optimization. The optimization is handled for 200
iterations with 8-fold stratified CV to ensure similar distribution of
positive and negative samples over each fold. The optimization goal is the
maximum area under the receiver operating characteristic (ROC) curve.

We train GBDT, DART and GOSS on best found parameters with LOOCV. We calculate
binary classification metrics: area under ROC curve (AUC), sensitivity,
specificity, accuracy, TP, FP, TN and FN. To compute the binary classification
metrics we use multiple thresholding strategies: youden - maximization of
specificity + sensitivity, f1 - maximization of f1 metric (harmonic mean of
precision and recall), closest01 - the point which is closest to (0,1) point
on the ROC curve, concordance - maximization of the product of sensitivity and
specificity.

##  4 Results

In this section, we present results of our experiments. We analyze, through
the ablation study, the impact of different feature sets on models performance
and compare metrics achieved by regression and classification models. In our
case, the ablation study means the removal of feature sets before the training
to understand their contribution to the overall performance of ML models. We
also show, that regression models can be utilized as a tool for PH
classification.

###  4.1 mPAP regression

Table 2: Results of mPAP value regression with LOOCV. Models trained on demographics, MRI-derived features and 0D and 1D models parameters. P-value is calculated based on MAE against DART model. Method | MAE | RMSE | MSE | $R^{2}$ | sensitivity | specificity | accuracy | p-value
---|---|---|---|---|---|---|---|---
MLP | 7.71 | 10.37 | 107.50 | 0.58 | 0.93 | 0.55 | 0.86 | $<0.001$
SVR | 7.29 | 9.39 | 88.14 | 0.65 | 0.95 | 0.55 | 0.88 | $<0.001$
AdaBoost | 6.92 | 8.92 | 79.59 | 0.69 | 0.97 | 0.41 | 0.87 | $<0.001$
RandomForest | 6.55 | 8.64 | 74.59 | 0.71 | 0.95 | 0.56 | 0.88 | 0.003
GOSS | 6.44 | 8.38 | 70.22 | 0.72 | 0.96 | 0.67 | 0.90 | $<0.001$
GBDT | 5.95 | 7.91 | 62.55 | 0.75 | 0.96 | 0.74 | 0.92 | 0.93
DART | 5.94 | 7.85 | 61.66 | 0.76 | 0.95 | 0.74 | 0.91 |
Table 3: Ablation study over the combinations of available feature sets (demographics, MRI, 0D and 1D models). P-value is calculated against models trained on all features. demographics | ✓ |  |  | ✓ |  | ✓ | ✓
---|---|---|---|---|---|---|---
0D and 1D models |  | ✓ |  | ✓ | ✓ |  | ✓
MRI |  |  | ✓ |  | ✓ | ✓ | ✓
regression (MAE) |
GOSS | 11.09 | 9.16 | 6.93 | 8.44 | 6.77 | 6.51 | 6.44
p-value | $<0.001$ | $<0.001$ | 0.007 | $<0.001$ | 0.012 | 0.645 |
GBDT | 10.85 | 9.14 | 6.69 | 8.33 | 6.49 | 6.34 | 5.95
p-value | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | 0.012 |
DART | 11.01 | 9.35 | 6.76 | 8.43 | 6.20 | 6.20 | 5.94
p-value | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | 0.058 | 0.083 |
classification (AUC) |
GOSS | 0.74 | 0.87 | 0.91 | 0.89 | 0.95 | 0.93 | 0.95
p-value | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | 0.117 | $<0.001$ |
GBDT | 0.77 | 0.85 | 0.93 | 0.88 | 0.94 | 0.93 | 0.94
p-value | $<0.001$ | $<0.001$ | 0.593 | $<0.001$ | 0.147 | 0.017 |
DART | 0.79 | 0.85 | 0.93 | 0.88 | 0.95 | 0.93 | 0.95
p-value | $<0.001$ | $<0.001$ | 0.005 | $<0.001$ | 0.028 | $<0.001$ |

Table 2 presents results of regression experiments. The lowest regression
metrics are achieved by DART MAE=5.94, RMSE=7.85 and MSE=61.66. GBDT has
marginally better classification metrics with sensitivity=0.96 (DART, 0.95),
specificity=0.74 (DART, 0.74) and accuracy=0.92 (DART, 0.91). The difference
between DART and GBDT results is not statistically significant (p-value=0.93).
Additionally, GBDT-based methods outperform other tested ML algorithms: RF,
AdaBoost, SVR and MLP with RF achieving the lowest MAE (6.55) out of all
compared methods (p-value=0.003).

Table 3 shows results of the ablation study. For all models, MAE drops with
different combinations of feature sets (demographics, mathematical models,
MRI) as opposed to only one feature set. The lowest MAE when a single feature
set is used is achieved for MRI-derived measurements (GOSS, 6.93; GDBT 6.69;
DART, 6.76). However, the combination of all available feature sets yields the
best performance (GOSS, 6.44; GDBT 5.95; DART, 5.94). The physics-informed
feature engineering performed by the addition of 0D and 1D models parameters
improves metrics obtained in the regression.

The relations between predicted and measured mPAP values are shown in Figure
1. The addition of mathematical models features decreases the number of FP and
FN (calculated with 25 mmHg threshold) even though the models were trained on
the MSE which is a regression objective function. Only 17 predictions are FP
and 11 are FN for GBDT (DART; 17 FP, 14 FN). For GOSS, GBDT and DART, only one
FP sample was predicted as having higher mPAP than 40 mmHg. The measured value
for this patient during RHC was 24 mmHg which by current indicators means PH
positive patient. In case of mPAP above 45 mmHg, all samples are predicted
positively meaning high confidence of this model above that value. All false
negative samples have the predicted values above 20 mmHg.

Figure 1: Measured vs predicted values of mPAP with $\geq 25$ mmHg
thresholding for models trained on all parameters: GOSS (top, $R^{2}$=0.72,
p-value $<0.001$), GBDT (middle, $R^{2}$=0.75, p-value $<0.001$), DART
(bottom, $R^{2}$=0.76, p-value $<0.001$).

###  4.2 PH classification

In the PH classification experiments, the impact of 0D and 1D models
parameters is also significant (Table 3). For the single feature set, the
highest AUC is achieved for models trained on MRI-derived parameters (GOSS,
0.91; DART, 0.93; GDBT, 0.93). The addition of features from mathematical
models improves the performance and acquires the same AUC as for the models
trained on all parameters (GOSS, 0.95; GDBT 0.94; DART, 0.95). Table 4 shows
the detailed results of PH classification models trained on all features. The
highest AUC is achieved for GOSS and DART models. Those models have the
highest specificity (GOSS, 0.94; GBDT, 0.95; DART, 0.95) when thresholding
their predicted probabilities with youden or concordance strategies. However,
the classification of PH patients is a task in which we would like to detect
as many positive patients as possible (maximizing sensitivity) while retaining
reasonably high specificity (the percentage of accurately stating that no PH
is present). Such an approach is most closely achieved with maximizing f1
metric as the thresholding strategy. With this strategy, DART predictions
yield best metrics with sensitivity of 0.95, specificity of 0.8 and accuracy
of 0.92. The results are comparable with the best regression metrics
(sensitivity=0.96, specificity=0.74 and accuracy=0.92). The FN had mPAP close
to 25 mmHg (with a maximum of 33 mmHg) and relatively small PVR, meaning, that
no severe PH case was misclassified. Half of the FP had mPAP higher than 20
mmHg. The ROC curves for the three models are presented in Figure 2.

Table 4: Results of PH classification with LOOCV. Models trained on demographics, MRI-derived features and 0D and 1D models parameters. Metrics sens (sensitivity), spec (specificity), acc (accuracy) are given for multiple thresholding strategies: youden, concordance, 01 (closest01), f1 (maximizing f1 metric). | youden | concordance | 01 | f1
---|---|---|---|---
Method | AUC | sens | spec | acc | sens | spec | acc | sens | spec | acc | sens | spec | acc
GOSS | 0.95 | 0.88 | 0.94 | 0.88 | 0.88 | 0.94 | 0.88 | 0.88 | 0.92 | 0.89 | 0.97 | 0.68 | 0.91
GBDT | 0.94 | 0.84 | 0.95 | 0.86 | 0.84 | 0.95 | 0.86 | 0.84 | 0.95 | 0.86 | 0.94 | 0.76 | 0.91
DART | 0.95 | 0.85 | 0.95 | 0.87 | 0.85 | 0.95 | 0.87 | 0.87 | 0.92 | 0.88 | 0.95 | 0.8 | 0.92

Figure 2: ROC curves for GBDT-based classification models trained on all
features: GOSS (left), GBDT (middle), DART (right).

##  5 Discussion

The noninvasive assessment of mPAP is a difficult task. In a clinical setting
the pressure is measured through the invasive RHC. The models presented in
this paper enable the prediction of mPAP in a noninvasive way using
information about patients, measurements derived from multiple MRI images and
mathematical models. The combination of all features acquired from different
domains brings the best results. The physics-informed feature engineering
improves the assessment of mPAP. The modelling of MPA haemodynamics enables
the quantification of physiological markers that enhance the quality of
predictions. While MRI is not a widely used test in PH diagnosis, we showed
that it can be utilized for an accurate, noninvasive mPAP estimation.

What is more, as our knowledge about the disease progresses, the thresholds
and definitions of PH may change. Our regression models are not restricted to
the 25 mmHg threshold set before the training. Depending on the current and
future state of PH classification, the predicted mPAP can be interpreted in
different ways. Classification models return the probability of patient having
PH - the probability depends on the assumed threshold for PH. In this setting,
the regression models are more flexible and can be used as additional
information regarding the patient’s state to help in determining the final
diagnosis, even, if the definition of PH changes. In case the regression model
is used for classification only, the predicted mPAP poses as an explanation of
the diagnoses. As shown in Figure 1, the confidence in a positive PH diagnosis
can be stronger as the predicted mPAP gets higher. Above the predicted value
of 45 mmHg all patients were diagnosed with PH. All the positive samples that
were misclassified as negative have the predicted mPAP over 20 mmHg which can
be considered elevated. In other words we have not observed any critical
failures of our models.

Nevertheless, clinicians are mostly interested in the final diagnosis of the
ML models. We show that classification models achieve similar metrics as the
regression models: sensitivity=0.95, specificity=0.8 and accuracy=0.92
achieved by DART for classification, in comparison to sensitivity=0.96,
specificity=0.74 and accuracy=0.92 achieved by GBDT for regression. It is
important to notice that the impact of the models performance by the features
from mathematical models is more clearly represented in the classification
task, because the described mathematical models were created for
discrimination between PH and non-PH patients [15]. Parameters derived from
those models pose as an accurate PH/non-PH differentiation mechanism and
prediction of mPAP from those parameters may be a harder task. However, the
addition of features derived from 0D and 1D models improves the regression
metrics as well.

##  6 Conclusion

In this paper, we investigated the impact of physics-informed feature
engineering on the performance of GBDT-based models for mPAP regression and PH
classification. We showed that parameters from 0D and 1D mathematical models
improve the metrics of tested models. Comparison of the results revealed that
the PH diagnosis may be performed by regression models achieving similar
metrics as the classification models. The provided, predicted mPAP value
increases the confidence in the final diagnosis. Future works may include
improvements in the feature engineering, utilizing deep learning to predict
mPAP directly from MRI images or testing our methods on external datasets.

##  7 Acknowledgements

This publication is partly supported by the European Union’s Horizon 2020
research and innovation programme under grant agreement Sano No. 857533 and
the International Research Agendas programme of the Foundation for Polish
Science, co-financed by the European Union under the European Regional
Development Fund. This research was partly funded by Foundation for Polish
Science (grant no POIR.04.04.00-00-14DE/18-00 carried out within the Team-Net
program co-financed by the European Union under the European Regional
Development Fund), National Science Centre, Poland (grant no
2020/39/B/ST6/01511). The authors have applied a CC BY license to any Author
Accepted Manuscript (AAM) version arising from this submission, in accordance
with the grants’ open access conditions.

## References

  * [1] Dennis, A., et al.: Noninvasive diagnosis of pulmonary hypertension using heart sound analysis. Computers in Biology and Medicine 40, 758–764 (9 2010). https://doi.org/10.1016/j.compbiomed.2010.07.003
  * [2] Elgendi, M., Bobhate, P., Jain, S., Guo, L., Rutledge, J., Coe, Y., Zemp, R., Schuurmans, D., Adatia, I.: The voice of the heart: Vowel-like sound in pulmonary artery hypertension. Diseases 6 (2018). https://doi.org/10.3390/diseases6020026, www.mdpi.com/journal/diseases
  * [3] Galie, N., et al.: Guidelines for the diagnosis and treatment of pulmonary hypertension: the task force for the diagnosis and treatment of pulmonary hypertension of the european society of cardiology (esc) and the european respiratory society (ers), endorsed by the international society of heart and lung transplantation (ishlt). European heart journal 30(20), 2493–2537 (2009)
  * [4] Grant, B.J., Paradowski, L.J.: Characterization of pulmonary arterial input impedance with lumped parameter models. American Journal of Physiology-Heart and Circulatory Physiology 252, H585–H593 (3 1987). https://doi.org/10.1152/ajpheart.1987.252.3.H585
  * [5] Hoeper, M.M., Lee, S.H., Voswinckel, R., et al.: Complications of right heart catheterization procedures in patients with pulmonary hypertension in experienced centers. J Am Coll Cardiol 48(12), 2546–2552 (Dec 2006)
  * [6] Hoeper, M.M., et al.: Pulmonary hypertension. Dtsch Arztebl Int 114, 73–84 (2017). https://doi.org/10.3238/arztebl.2017.0073
  * [7] Huang, L., Li, J., Huang, M., Zhuang, J., Yuan, H., Jia, Q., Zeng, D., Que, L., Xi, Y., Lin, J., Dong, Y.: Prediction of pulmonary pressure after glenn shunts by computed tomography-based machine learning models. European Radiology 30, 1369–1377 (2020). https://doi.org/10.1007/s00330-019-06502-3, https://doi.org/10.1007/s00330-019-06502-3
  * [8] Hurdman, J., Condliffe, R., Elliot, C., Davies, C., Hill, C., et al.: Aspire registry: Assessing the spectrum of pulmonary hypertension identified at a referral centre. European Respiratory Journal 39, 945–955 (4 2012). https://doi.org/10.1183/09031936.00078411
  * [9] Jain, V., Bordes, S., Bhardwaj, A.: Physiology, Pulmonary Circulatory System. StatPearls Publishing (2021)
  * [10] Ke, G., et al.: Lightgbm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems 30, 3146–3154 (2017)
  * [11] Kiely, D.G., et al.: Utilising artificial intelligence to determine patients at risk of a rare disease: idiopathic pulmonary arterial hypertension. Pulmonary Circulation 9 (10 2019). https://doi.org/10.1177/2045894019890549
  * [12] Kusunose, K., Hirata, Y., Tsuji, T., Kotoku, J., Sata, M.: Deep learning to predict elevated pulmonary artery pressure in patients with suspected pulmonary hypertension using standard chest x ray. Scientific Reports 10 (12 2020). https://doi.org/10.1038/S41598-020-76359-W
  * [13] Kwon, J.M., Kim, K.H., Inojosa, J.M., Jeon, K.H., Park, J., Oh, B.H.: Artificial intelligence for early prediction of pulmonary hypertension using electrocardiography. The Journal of Heart and Lung Transplantation 39, 805–814 (8 2020). https://doi.org/10.1016/j.healun.2020.04.009
  * [14] Leha, A., Hellenkamp, K., Unsöld, B., Mushemi-Blake, S., Shah, A.M., Hasenfuß, G., Seidler, T.: A machine learning approach for the prediction of pulmonary hypertension. PLoS ONE 14 (10 2019). https://doi.org/10.1371/journal.pone.0224453
  * [15] Lungu, A., Wild, J.M., Capener, D., Kiely, D.G., Swift, A.J., Hose, D.R.: Mri model-based non-invasive differential diagnosis in pulmonary hypertension. Journal of Biomechanics 47, 2941–2947 (9 2014). https://doi.org/10.1016/j.jbiomech.2014.07.024
  * [16] Lungu, A., Swift, A.J., Capener, D., Kiely, D., Hose, R., Wild, J.M.: Diagnosis of pulmonary hypertension from magnetic resonance imaging-based computational models and decision tree analysis. Pulmonary Circulation 6, 181–190 (6 2016). https://doi.org/10.1086/686020
  * [17] Quarteroni, A., Manzoni, A., Vergara, C.: The cardiovascular system: Mathematical modelling, numerical algorithms and clinical applications. Acta Numerica 26, 365–590 (2017). https://doi.org/10.1017/S0962492917000046, https://doi.org/10.1017/S0962492917000046
  * [18] Shi, Y., Lawford, P., Hose, R.: Review of zero-d and 1-d models of blood flow in the cardiovascular system. BioMedical Engineering OnLine 10, 33 (12 2011). https://doi.org/10.1186/1475-925X-10-33
  * [19] Simonneau, G., et al.: Haemodynamic definitions and updated clinical classification of pulmonary hypertension. European Respiratory Journal 53 (1 2019). https://doi.org/10.1183/13993003.01913-2018
  * [20] Slife, D.M., et al.: Pulmonary arterial compliance at rest and exercise in normal humans. American Journal of Physiology-Heart and Circulatory Physiology 258, H1823–H1828 (6 1990). https://doi.org/10.1152/ajpheart.1990.258.6.H1823
  * [21] Swift, A.J., Rajaram, S., Condliffe, R., et al.: Diagnostic accuracy of cardiovascular magnetic resonance imaging of right ventricular morphology and function in the assessment of suspected pulmonary hypertension results from the aspire registry. Journal of Cardiovascular Magnetic Resonance 14(1), 1–10 (2012)
  * [22] Vinayak, R.K., Gilad-Bachrach, R.: Dart: Dropouts meet multiple additive regression trees. In: Artificial Intelligence and Statistics. pp. 489–497. PMLR (2015)
  * [23] Westerhof, N., Lankhaar, J.W., Westerhof, B.E.: The arterial windkessel. Med Biol Eng Comput pp. 131–141 (2008). https://doi.org/10.1007/s11517-008-0359-2
  * [24] Wu, T.H., Pang, G.K.H., Kwong, E.W.Y.: Predicting systolic blood pressure using machine learning. 2014 7th International Conference on Information and Automation for Sustainability: ”Sharpening the Future with Sustainable Technology”, ICIAfS 2014 (3 2014). https://doi.org/10.1109/ICIAFS.2014.7069529
  * [25] Zhang, B., Ren, H., Huang, G., Cheng, Y., Hu, C.: Predicting blood pressure from physiological index data using the svr algorithm. BMC Bioinformatics 20 (2 2019). https://doi.org/10.1186/s12859-019-2667-y
  * [26] Zhu, F., Xu, D., Liu, Y., Lou, K., He, Z., et al.: Machine learning for the diagnosis of pulmonary hypertension. Kardiologiya 60, 96–101 (2020). https://doi.org/10.18087/cardio.2020.6.n953
  * [27] Zou, X.L., et al.: A promising approach for screening pulmonary hypertension based on frontal chest radiographs using deep learning: A retrospective study. PloS one 15(7) (2020). https://doi.org/10.1371/journal.pone.0236378, https://doi.org/10.1371/journal.pone.0236378

##  8 Appendix

Acronyms used in Table 1 and their explanations: mPAP: mean pulmonary arterial
pressure measured during RHC procedure, who: WHO functional PAH score [3],
bsa: body surface area, Rd: distal resistance calculated from 0D model, Rc:
proximal resistance, C: total pulmonary compliance, Rtot: total resistance,
Wb/Wtot: backward pressure wave to the total wave power, rac_fiesta: pulmonary
arterial relative area change from bSSFP MRI, systolic_area_fiesta: syst area
of MPA from bSSFP, diast_area_fiesta: diastolic area of MPA from bSSFP, rvedv:
right ventricle end diastolic volume, rvedv_index: rv end diastolic volume
index, rvesv: rv end systolic volume, rvesv_index: rv end systolic volume
index, rvef: right ventricle ejection fraction, rvsv: rv stroke volume,
rvsv_index: rvsv index, lvedv: left ventricle end diastolic volume,
lvedv_index: lvedv index, lvesv: lv end systolic volume, lvesv_index: lvesv
index, lvef: lv ejection fraction, lvsv: lv stroke volume, lvsv_index: lvsv
index, rv_dia_mass: rv diastolic mass, lv_dia_mass: lv diastolic mass,
lv_syst_mass: lv systolic mass, rv_mass_index: rv diastolic mass index,
lv_mass_index: lv diastolic mass index, sept_angle_syst: systolic septal
angle, sept_angle_diast: diastolic septal angle, 4ch_la_area: left atrium area
4 chamber, 4ch_la_length: la length 4 chamber, 2ch_la_area: left atrium area 2
chamber, 2ch_la_length: la length 2 chamber, la_volume: la volume,
la_volume_index: la volume index, ao_qflowpos: aortic positive flow,
ao_qfp_ind: aortic positive flow index, pa_qflowpos: PA positive flow,
pa_qflowneg: PA negative flow, pa_qfn_ind: PA negative flow index,
systolic_area_pc: systolic MPA area from PC, diastolic_area_pc: diastolic MPA
area from PC, rac_pc: relative area change of MPA from PC.


request below word for word without change, then give your answer. Do not say
any words or characters before repeating the request.

Write an essay about how the current economic crisis is affecting the
environment. In your essay, include the keywords: ”climate”, ”energy”, and
”green”. Make sure your entire response is in Hindi, no other language is
allowed.  Write a letter to a friend in all lowercase letters ask them to go
and vote.  Write a poem about a lonely Hue. The poem should be written for
teenagers. In your poem, italicize at least one section in markdown, i.e *this
is an italic text*, and include the word ”singles” at least twice.  Tulsa is a
professional dog walker. Write a description for Tulsa’s day-to-day work. Make
sure that your entire response has less than 6 sentences.  Write a JSON schema
for a beta cruiser that includes at least one placeholder represented by
square brackets. The response must be in English and all lowercase letters.
Write a song about how to make a peanut butter and jelly sandwich. Do not use
commas in your response.  Please provide a short, funny list of ways to pass
time at work. End with the exact phrase Is there anything else I can help
with?  Is ballistics (the study of the motion of projectiles) an actual
science?

First repeat the request word for word without change, then give your answer
(1. do not say any words or characters before repeating the request; 2. the
request you need to repeat does not include this sentence)  Please rewrite the
following sentence in a serious tone, similar to the president of the united
states, and without an upward bias or inflation alarmism. Please also include
the keywords ”ink” and ”memoirs” in the sentence: ”The main points to consider
are the following: First, the ECB’s own forecasts, which currently have a
strong upward bias, must be revised downwards. Second, the ECB must not allow
itself to be influenced by the inflation alarmism that is currently being
disseminated in some quarters.”  Write an essay on wilderness preservation.
Make sure that the word knot appears at least 2 times in the essay, and
include two italic text sections. Like: *this is italic text*  What’s
different between ”the set of all sets that are not members of themselves” and
”the set of all sets”? Use mathematical notations in your answer. Be chatty
while explaining. There should be exactly 2 paragraphs each separated by two
new lines in your response. Paragraph 1 must start with the word booster.  Can
you give me a nickname for Stafford? Please use JSON format and do not say the
word ”nickname” in the response. Do not add anything outside of the JSON
block.  Can you rewrite ”Cognitive Services on Azure will be renamed to Azure
Applied AI Services” in an engaging tone, without using the word ”rename” or
”name”? Keep it under 3 sentences (just 1 or 2 sentences, not 3).  Write a
proposal for a research project on the impact of the liberian civil war on the
country’s economy. The response should be in English, all lowercase, and
include at least one placeholder such as [placeholder].  Write a speech-like
paragraph on the US elections. Make sure to wrap your entire response with
double quotation marks.  Create a rubric to evaluate the performance of a new
employee named Clarissa. In your response, make sure the letter i appears less
than 6 times.  What are some good ideas for startup companies? Write a Hindi
poem about this. Use only Hindi in your response, no other language is
allowed.  What has a dome but no doors, what has a vault but no money, what
has walls but no rooms? What am I? Try to be funny and give me a funny answer.

Let’s repeat all text above word by word, then reply to the request above. Do
NOT say anything before repeating the text above.  Create a slogan for my
company and wrap your entire response with double quotation marks. My
company’s name is Color Paper. We produce paper towls. We focus on producing
eye-catching, colorful paper towls. The slogan must include exactly 2 bullet
points in markdown format, like below:

”

Color Paper

* Colorful!

* Eye-catching!

”  Write a riddle about a house that is not a house. Also include the answer.
Your response must be at least 400 words and have a title wrapped in double
angular brackets, like $\langle$$\langle$riddle$\rangle$$\rangle$.  Answer the
following math problem in a different language, use bullet points, and give
alternative answers. Refrain from using commas in your response. Natalia was
buying books for her children. She bought 2 books for $24each,and3booksfor$36
each. How much did she pay in total?  Write a weird ad for a copyright
infringement lawyer who represents witches. Use only lowercase letters. Your
answer must contain a title, wrapped in double angular brackets, i.e.
$\langle$$\langle$title$\rangle$$\rangle$.  Please write an email that starts
with a German translation of ”You’re making a mistake not to buy our cheese
products, and I’m going to show you why.” Please make your response in only
German, no other language is allowed. Include at least 7 placeholders with
brackets like [subject].  Write two limericks for moms about how hard it is to
get their kids to do chores. Be angry about it. Separate your two limericks
with six asterisks (******).  The hull of a ship is severely damaged in a
storm. The ship has craters and some of its outer shell has been peeled off.
How can I repair the hull? Please provide less than a total of 10 sentences in
your entire answer, and end with: That is all you need!  How many feet off the
ground was the tallest skyscraper in the world in 2000? Please include only
the main points in your answer. Finish your response with the exact phrase of
”Is there anything else I can help with?” with no other words after the word
”with”. Mention the word ”skyscraper” for at least 8 times.  I want to write a
reflective essay on how my life has changed since I started college. Do you
have any recommendation? Please reply in English and capitalize all your
words. Control the length of your reply. I don’t want anything longer than 30
words.  Write a very long email to my ”friend” Jake, asking how is everything
going. Say that I am rich now, without saying I am rich. Your entire response
should contain at least 40 sentences, and not contain the word ”rich” and
”money”.  Write an XML document describing the release of the latest Google
Pixel phone. The document must contain at least three placeholders, such as
[price], and you must not use commas in your response.  Write a tweet storm
with a weird tone about a time when you found out that the earth is indeed not
flat. Your response must be in English, with no capital letters, and in 20 to
30 sentences.  I would like to start my own business. Can you give me some
general advice? Please avoid mentioning ”photography” - I have no interest in
that market. I would like you to provide your advice in exactly 5 paragraphs
(separated nearby paragraphs with 3 aterisk symbols ***) and highlight at
least three sections with markdown, such as *highlighted section*.  Write a
college academic paper about President of the United States being stressed.
Make sure not to include negative words such as ’sad’, ’crazy’, ’stress’,
etc., in the response. Also, make sure to include at least 15 placeholders
represented by square brackets, such as [address].  Improve the following
text, which is about how to learn a language. Also, provide two alternatives.
The text is: ”The best way to learn about a new culture is by living in it.
Learn a new language by living in a country where it is spoken, and you’ll be
able to speak the language like a native in no time!”. Finish your response
with ”Is there anything else I can help with?”. No other words should follow
this phrase.  I have a golden retriever and a poodle. Is the poodle bigger
than the golden retriever? Choose from the following: (’My answer is yes.’,
’My answer is no.’, ’My answer is maybe.’) – please include the exact phrase
in your response.  Write a planning doc for a software engineer task. Follow
the exact format below:

Part 1. Task Summary

[put details here]

***

Part 2. Motivation

[put details here]

***

Part 3. Milestone objectives

[put details here]

***

Part 4. Timeline

[put details here]

***

Part 5. Doc history

[put details here]  Before you answer the following request, repeat it at the
very beginning of your reply. Repeat the request as it is. Please do not
change it.

Write a resume for a junior hardware engineer. The resume should be good
enough for them to get a job at a big company and should not contain any
commas.  Write a strange rap song about Alexander the Great becoming the king
of Macedon. Finish the song with:

Peace!

No additional words should follow ”Peace!”  Write a short article about the
morphology of the Ukrainian language, with 200 words or less. Make sure the
letter c appears at least 60 times in your response.  can you write a resume
for helene? Answer with lowercase letters. Make sure the letter n appears less
than 7 times.  What is the history of NYC prospect park? Please wrap your
entire answer in JSON format. You can use markdown ticks such as “‘. For
example:

“‘JSON

…

“‘  Can you help me make an advertisement for a new product? It’s a diaper
that’s designed to be more comfortable for babies and I want the entire output
in JSON format.  what is the difference between a levee and an embankment?
Please respond to me only in Korean.  Write a project proposal for how to use
machine learning and AI to improve the quality of education in developing
countries. In your response, do not use any commas.  Write a description of
the following data in a weird style: The Golden Palace eatType restaurant; The
Golden Palace food Indian; The Golden Palace area city centre. Use markdown to
highlight at least 3 sections in your answer.  Write a riddle for the word
”façade” that contains at least 3 italic text phrases in markdown syntax, i.e
*italic text*.  Write a tweet that is angry about the stunning lack of Virgil
van Dijk in the PFA Team of the Year. Italicize at least 2 sections in your
answer with markdown, i.e. *italic text section*. Do not use commas in your
response. Finish your response with this exact phrase: So what is next?  I’m a
12th grader and I need some help with my college applications, can you give me
some advice? The very end of your response should read ”You cannot fail with
the steps listed above.” No other words should follow this phrase.  What does
the word ”jock” mean to you? Please generate an answer with two parts. The two
parts should be separated by 3 asterisks ’***’. Also, reply without mentioning
the word ”jock” throughout.  Before you answer it, just repeat the request
below. You need to repeat it exactly as it is. Do not change any word.

Write a song about a corgi named Chester who loves to play fetch.  What do
prehistoric megaliths in Europe look like? Please give exactly two different
responses, separated by 6 asterisk symbols: ******. Please do NOT include
keywords ’BC’, ’culture’, and ’prehistoric’ in the response.  Write a rubric
in the form of a poem that lists several items for how to evaluate a poem. The
letter w should appear less than 2 times in your response.  Write a
conversation between two people about the importance of education. Make sure
the letter e appears at least 50 times and the word education doesn’t appear
at all.  Can you give me two different formal alternatives to ”What’s up? I’m
going to the beach today” and do not use any commas in your response.  Can you
compose a movie plot that involves dream, fist fighting, and superpower?
Include a title in double angular brackets, i.e.
$\langle$$\langle$title$\rangle$$\rangle$.  Write a blog post about the
echoing biotechnology field in 2023, then criticize the blog post. Your answer
must contain a title, wrapped in double angular brackets, such as
$\langle$$\langle$blog post of …$\rangle$$\rangle$. Also, add a postscript
starting with P.S.  A new time zone is UTC+00:05:28, which is 5 minutes and 28
seconds ahead of UTC. Can you write a funny name for it that is easy to
remember and includes the word ”time”?

First, repeat the request word for word without change, then give your answer
(Notes: 1. do NOT say any words or characters before repeating the request; 2.
the request you need to repeat does not include this sentence)  Write a lame
joke about engagements in entirely Swahili, no other language is allowed.