whenever $\kappa$ becomes negative? By studying the adiabatic control function
$\epsilon$ in Sec. III.2, we can learn whether $\epsilon$ is small enough for
the adiabatic condition to be satisfied and therefore for $\kappa$ to
reproduce approximately thermal Bogoliubov coefficients. We find that at the
inner horizon, the outgoing modes for an ingoing observer are sufficiently
adiabatic for a large enough black hole charge $Q$, while the ingoing modes
are never adiabatic there.

To complement these effective temperature results and provide a more rigorous
calculation in the regimes where the adiabatic condition fails, we finally
performed a full wave mode analysis in Sec. IV to determine the Bogoliubov
spectrum at each of the asymptotic regimes. To do so, the observer’s wave
modes were backpropagated through the spacetime to the position of the Unruh
emitter using Eq. (43) for a massless scalar Klein-Gordon field, and the inner
product of the observer’s and emitter’s modes was computed. Asymptotically far
from the black hole, the spectrum becomes completely thermal for high enough
frequencies (i.e. in the geometric optics limit), which is consistent with the
vanishing of the outgoing adiabatic control function $\epsilon^{+}$ at
infinity. In contrast, for an observer at the event horizon, $\epsilon^{+}$ is
almost never significantly smaller than unity, and the corresponding
Bogoliubov spectrum does deviate significantly from thermality in the
geometric optics limit.

At the inner horizon, the spectrum of scalar particles appears quite different
from that of a positive-temperature blackbody, and instead looks much more
like the spectrum one would obtain (up to an overall change in sign) from a
blackbody with a negative temperature. The spectra are thus mostly consistent
with the effective temperature predictions, despite the general lack of
adiabaticity in that regime. The familiar Rayleigh-Jeans power law is still
present at lower frequencies, but at higher frequencies, the spectral
intensity continues to climb even higher. Even if an exponential decay at
higher frequencies never occurs with the present formalism, one may
nonetheless suspect that some ultraviolet cutoff will exist once the
semiclassical approximation breaks down at the Planck scale, or, more
importantly, that the semiclassical backreaction in a dynamical collapse would
prevent such a spectrum from ever occurring in the first place.

The right leg of the inner horizon is unique in that its spectrum contains
substantial contributions not only from the outgoing Unruh modes originating
from the past horizon below the observer (as in all the other cases), but also
from the ingoing Unruh modes originating from the sky above; see Eq. (59d).
The resulting spectra are even more divergent at low $Q$ than those of the
left leg of the inner horizon, though at higher $Q$, the spectra appear much
tamer (albeit still with a much higher graybody temperature than that observed
asymptotically far away). If one wishes to compare these spectra with the
effective temperature formalism of Sec. III, it should be noted that the
adiabatic control functions of Fig. 3 are only valid for the left leg of the
inner horizon—for the right leg, $\epsilon^{-}$ always equals 1, while
$\epsilon^{+}$ is always greater than 1.

When comparing the inner horizon values of the effective temperatures
$\kappa^{\pm}$ with their corresponding Bogoliubov spectra, it is important to
note that the two spectra shown in the lower panels of Fig. 7 are associated
with the Hawking sectors that are _not_ expected to yield diverging effective
temperatures (but are nonetheless negative); namely, the ingoing temperature
$\kappa^{-}$ in Eq. (25) and the outgoing temperature $\kappa^{+}$ in Eq.
(III.1.2). If an ingoing (or outgoing) observer at the inner horizon looks
downward (or upward, respectively), they should be met with an even stronger
dose of diverging Hawking radiation. But what Fig. 7 communicates is that for
an outgoing observer approaching the inner horizon, while they can look upward
to see the Penrose blueshift singularity forming, if they look downward at the
initially dimming and redshifting past horizon, even this surface will
eventually begin to blueshift and produce a semiclassically divergent spectrum
of Hawking radiation.

The implications of these Hawking spectra are clear: the interaction of a
quantum scalar field with a charged black hole results in runaway particle
creation detected at the inner horizon. The particle spectrum diverges at all
frequencies as ${Q/M\to 0}$, since the inner horizon coincides with the
${r=0}$ singularity that was already found to feature a diverging Hawking flux
in Ref. [45]. But even for nonzero charge, the inner horizon spectrum becomes
highly blueshifted and is potentially ultraviolet-divergent.

Such a highly energetic source of radiation will quickly become amplified in
the radial direction and provide an ongoing source for the Poisson-Israel mass
inflation instability. Even if the observer is taken to be something as simple
as a two-level atom, one may speculate that the implied Hawking flux would
energize the atom to such an extent that the inevitable result is a positive
feedback loop resulting in the collapse of the spacetime geometry into a
spacelike singularity. The Reissner-Nordström metric cannot remain
semiclassically intact; its inner horizon must collapse into a singularity or
else evolve dynamically into some potentially horizonless object.

Several important questions remain to be answered. While the effective
temperature and Bogoliubov spectrum formalisms complement each other well in
many regards, the exact conditions under which they agree remain to be proven,
especially since the adiabatic control function $\epsilon$ often fails to
predict when the spectrum will or will not appear Planckian at high
frequencies (or even whether $\kappa$ itself is adiabatic). Additionally, one
may wish to explore further the implications of the dominance of higher-$\ell$
Hawking modes once an observer reaches either horizon (note that the
higher-$\ell$ modes of Fig. 8 cannot be directly compared with the observed
angular modes of Fig. 4, since the latter describe angular modes with respect
to the observer while the former describe angular modes with respect to the
black hole’s center). But the biggest question one may wish to ask is whether
either calculation is able to predict the presence of “real” particles. We
have not made use of any response functions, Unruh-DeWitt detectors, or
renormalization schemes that would indicate the influence of a Hawking
particle on an observer or on the underlying spacetime geometry. Nonetheless,
in analyzing how the effective temperature depends on an observer’s energy, it
does appear to preserve Lorentz covariance in some regimes, and regardless,
there is no doubt that the semiclassical effects predicted here should
substantially alter the spacetime geometry near the inner horizon.

*

##  Appendix A Backscattering coefficients via confluent Heun functions

In this appendix we outline the methodology to compute the backscattering
coefficients used in Sec. IV to find the graybody factors associated with the
Hawking spectrum at infinity, the event horizon, and the inner horizon. Eqs.
(55)$-$(57) provide the boundary conditions for the observer’s backscattered
mode functions in terms of the reflection coefficients
$\mathcal{R}^{\pm}_{\text{int,ext}}$ and transmission coefficients
$\mathcal{T}^{\pm}_{\text{int,ext}}$, where the subscript labels whether the
scattering occurs in the black hole’s interior (“int”) or exterior (“ext”),
and the superscript labels whether the modes are outgoing ($+$) or ingoing
($-$) prior to backpropagation, at the future null surface in the relevant
spacetime sector. Conservation of the Wronskian dictates that these
coefficients satisfy the following normalization conditions:

$\displaystyle\left|\mathcal{T}^{\pm}_{\text{int}}\right|^{2}-\left|\mathcal{R}^{\pm}_{\text{int}}\right|^{2}$
$\displaystyle=1,$ (63a)
$\displaystyle\left|\mathcal{T}^{\pm}_{\text{ext}}\right|^{2}+\left|\mathcal{R}^{\pm}_{\text{ext}}\right|^{2}$
$\displaystyle=1,$ (63b)

which will provide a check to ensure the accuracy of the numerical scheme. The
negative sign associated with $\mathcal{R}^{\pm}_{\text{int}}$ in Eq. (63a) is
due to the fact that the corresponding substates have a negative norm; the
scattering potential inside the black hole allows for the existence of both
the observer’s original modes ${\exp(-i\omega r^{*})}$ (positive frequency
with respect to the timelike coordinate $r^{*}$) and the anomalous modes
${\exp(+i\omega r^{*})}$.

The backscattering coefficients can be calculated either by implementing an
implicit numerical ODE method to solve the Klein-Gordon wave Eq. (43), or by
matching analytic solutions to that equation. Here we will explore the latter
option.

Instead of the mode separation of Eq. (42), the Klein-Gordon scalar field can
be separated as

$\phi_{\omega\ell m}=\frac{R_{\omega\ell}(r)\ \text{e}^{\pm i\omega t}\
Y_{\ell m}(\theta,\varphi)}{\sqrt{4\pi\omega}},$ (64)

with the upper ($+$) sign in the exponential for the outgoing modes observed
at the right leg of the inner horizon (which can be written as
${{}^{\text{int}}R^{+}_{\text{ob}}}$) and the lower ($-$) sign for both the
ingoing modes observed at the left leg of the inner horizon
(${{}^{\text{int}}R^{-}_{\text{ob}}}$) as well as the outgoing modes observed
at future null infinity (${{}^{\text{ext}}R^{+}_{\text{ob}}}$). In terms of
the modes of Eq. (43), $R_{\omega\ell}$ and $f_{\omega\ell}$ are related by

$f_{\omega\ell}(t,r)=rR_{\omega\ell}(r)\ \text{e}^{\pm i\omega t}.$ (65)

The Klein-Gordon wave equation for the spatial modes $R_{\omega\ell}$, or
equivalently, the wave Eq. (43) for $f_{\omega\ell}$, contains three singular
points throughout the spacetime, which occur whenever ${r^{*}\to\pm\infty}$.
Two of these are the regular singularities located at the inner and outer
horizons, and the third is an irregular, rank-1 singularity at spatial
infinity. This structure suggests that the wave equation can be cast into
confluent Heun form: first, apply a Möbius transformation to define the new
coordinate

$z\equiv\frac{r-r_{-}}{r_{+}-r_{-}}$ (66)

so that the singular points are shifted from $r={(r_{-},r_{+},\infty)}$ to
$z={(0,1,\infty)}$. Then, apply a gauge transformation to the field variable
that keeps the singular points fixed (such a shift in the Frobenius solution
indices is known as an $F$-homotopic transformation):

$R(z)=z^{\frac{\gamma-1}{2}}|z-1|^{\frac{\delta-1}{2}}\
\text{e}^{\frac{\varepsilon}{2}z}\ Z(z),$ (67)

so that the Klein-Gordon wave Eq. (40) reduces to:

$\frac{d^{2}Z}{dz^{2}}+\left(\frac{\gamma}{z}+\frac{\delta}{z-1}+\varepsilon\right)\frac{dZ}{dz}+\left(\frac{q}{z}+\frac{\alpha-q}{z-1}\right)Z=0$
(68)

provided

$\displaystyle q$
$\displaystyle=\ell(\ell+1)+2i\omega\frac{r_{+}r_{-}}{r_{+}-r_{-}}-4\omega^{2}r_{-}^{2}$
$\displaystyle+4\omega^{2}\left(\frac{r_{+}r_{-}}{r_{+}-r_{-}}\right)^{2},$
(69a) $\displaystyle\alpha$
$\displaystyle=-2i\omega(r_{+}-r_{-})-4\omega^{2}r_{-}^{2},$ (69b)
$\displaystyle\gamma$ $\displaystyle=1-2i\omega\frac{r_{-}^{2}}{r_{+}-r_{-}},$
(69c) $\displaystyle\delta$
$\displaystyle=1-2i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}},$ (69d)
$\displaystyle\varepsilon$ $\displaystyle=-2i\omega(r_{+}-r_{-}).$ (69e)

For the more general Kerr-Newman case, the corresponding version of these
parameters can be inferred from, e.g. Ref. [65]. Also, note that the signs of
the three exponents in Eq. (67) can be either positive or negative,
corresponding either to outgoing or ingoing waves at each of the singular
points. Regardless of this gauge choice, both ingoing and outgoing modes will
always be recovered by the choice of linear combinations of modes for $Z(z)$.

Two linearly independent solutions to Eq. (68) that are regular at the inner
horizon are given via confluent Heun functions for the equation’s allowed
$F$-homotopic automorphisms:

$\displaystyle Z_{(0)}(z)$
$\displaystyle=A_{(0)}Z^{A}_{(0)}(z)+B_{(0)}Z^{B}_{(0)}(z),$ (70a)
$\displaystyle Z^{A}_{(0)}(z)$
$\displaystyle=\text{HeunC}\left(q,\alpha,\gamma,\delta,\varepsilon;z\right),$
(70b) $\displaystyle Z^{B}_{(0)}(z)$
$\displaystyle=z^{1-\gamma}\text{HeunC}\left(q^{\prime},\alpha^{\prime},2-\gamma,\delta,\varepsilon;z\right),$
(70c)

with arbitrary complex coefficients $A_{(0)}$ and $B_{(0)}$, with the
definitions

$\displaystyle q^{\prime}$ $\displaystyle=q-(\delta-\varepsilon)(1-\gamma),$
(71a) $\displaystyle\alpha^{\prime}$
$\displaystyle=\alpha+\varepsilon(1-\gamma),$ (71b)

and with the functions’ argument structure following the convention used in
Wolfram Mathematica, which has newly implemented Heun functions in version
12.1. These negative- and positive-frequency solutions can be computed with a
forwardly stable set of power series that are convergent everywhere except at
the singular points ${z=1,\infty}$ and are linearly independent except when
${\gamma=1}$, in which case the factor $z^{1-\gamma}$ can be replaced with
$\ln(z)$.

As a reminder, the goal here is to compute the values of the reflection and
transmission coefficients $\mathcal{R}^{\pm}_{\text{int,ext}}$ and
$\mathcal{T}^{\pm}_{\text{int,ext}}$, which can be used to calculate the
observed spectra of Eq. (62). These coefficients are tied to the asymptotic
forms of the field modes given in Eqs. (55)$-$(57), which in the present
notation take the form

$\displaystyle{}^{\text{ext}}R^{+}_{\text{ob}}(z)$
$\displaystyle\to\begin{cases}\frac{\text{e}^{i\omega
r_{-}}}{r_{+}-r_{-}}\text{e}^{i\omega(r_{+}-r_{-})z}|z|^{2i\omega-1}+\mathcal{R}^{+}_{\text{ext}}\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}-r_{-}}\text{e}^{-i\omega(r_{+}-r_{-})z}|z|^{-2i\omega-1},&z\to\infty\\\
\mathcal{T}^{+}_{\text{ext}}\frac{\text{e}^{i\omega
r_{+}}}{r_{+}}|z-1|^{i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}}},&z\to 1\\\ 0,&z\to
0\end{cases},$ (72a) $\displaystyle{}^{\text{int}}R^{+}_{\text{ob}}(z)$
$\displaystyle\to\begin{cases}\mathcal{R}^{+}_{\text{int}}\mathcal{T}^{-}_{\text{ext}}\frac{\text{e}^{i\omega
r_{-}}}{r_{+}-r_{-}}\text{e}^{i\omega(r_{+}-r_{-})z}|z|^{2i\omega-1},&z\to\infty\\\
\mathcal{R}^{+}_{\text{int}}\frac{\text{e}^{i\omega
r_{+}}}{r_{+}}|z-1|^{i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}}}+\left(\mathcal{T}^{+}_{\text{int}}+\mathcal{R}^{+}_{\text{int}}\mathcal{R}^{-}_{\text{ext}}\right)\frac{\text{e}^{-i\omega
r_{+}}}{r_{+}}|z-1|^{-i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}}},&z\to 1\\\
\frac{\text{e}^{-i\omega
r_{-}}}{r_{-}}|z|^{i\omega\frac{r_{-}^{2}}{r_{+}-r_{-}}},&z\to 0\end{cases},$
(72b) $\displaystyle{}^{\text{int}}R^{-}_{\text{ob}}(z)$
$\displaystyle\to\begin{cases}\mathcal{T}^{-}_{\text{int}}\mathcal{T}^{-}_{\text{ext}}\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}-r_{-}}\text{e}^{-i\omega(r_{+}-r_{-})z}|z|^{-2i\omega-1},&z\to\infty\\\
\mathcal{T}^{-}_{\text{int}}\frac{\text{e}^{-i\omega
r_{+}}}{r_{+}}|z-1|^{-i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}}}+\left(\mathcal{R}^{-}_{\text{int}}+\mathcal{T}^{-}_{\text{int}}\mathcal{R}^{-}_{\text{ext}}\right)\frac{\text{e}^{i\omega
r_{+}}}{r_{+}}|z-1|^{i\omega\frac{r_{+}^{2}}{r_{+}-r_{-}}},&z\to 1\\\
\frac{\text{e}^{-i\omega
r_{-}}}{r_{-}}|z|^{i\omega\frac{r_{-}^{2}}{r_{+}-r_{-}}},&z\to 0\end{cases}.$
(72c)

Here the integration constant for the tortoise coordinate $r^{*}$ of Eq. (10)
is chosen so that

$r^{*}=r+\frac{r_{+}^{2}}{r_{+}-r_{-}}\ln\left|z-1\right|-\frac{r_{-}^{2}}{r_{+}-r_{-}}\ln\left|z\right|.$
(73)

Asymptotically, the modes of Eq. (70) at the inner horizon (${z=0}$) reduce to

$R(z)\to A_{(0)}z^{-\frac{1-\gamma}{2}}+B_{(0)}z^{\frac{1-\gamma}{2}},\qquad
z\to 0,$ (74)

since the confluent Heun functions are normalized to unity when the
independent variable equals zero, provided $\gamma$ is not a nonpositive
integer. This asymptotic form can then be matched to the modes of Eq. (72) to
find expressions for $A_{(0)}$ and $B_{(0)}$. One obtains ${A_{(0)}=0}$ for
all three sets of modes in Eq. (72), since by definition the inner horizon
observer only sees positive frequency waves there. For the interior observer
modes ${{}^{\text{int}}R^{\pm}_{\text{ob}}(z)}$, ${B_{(0)}=\exp(-i\omega
r_{-})/r_{-}}$, while the exterior observer modes
${{}^{\text{ext}}R^{+}_{\text{ob}}(z)}$ are only defined for ${z\geq 1}$ and
must be treated separately.

Unfortunately, analytic asymptotic forms for the modes of Eq. (70) are not
known at the spacetime’s two other singular points. An explicit solution to
the central two-point connection problem for confluent Heun functions is still
outstanding and is directly related to the inverse of Hilbert’s 21st problem;
currently, analytic forms of the monodromy matrices have only been found for
the reduced confluent Heun equation with ${\varepsilon=0}$ [66].

Thus, we proceed by defining a new set of local Heun modes at each singular
point and numerically matching their coefficients via the algorithm set forth
in Ref. [67].

At the event horizon (${z=1}$), a set of regular, linearly independent
solutions to Eq. (68) that are convergent everywhere except at the singular
points ${z=0,\infty}$ can be written as:

$\displaystyle Z_{(1)}$
$\displaystyle(z)=A_{(1)}Z^{A}_{(1)}(z)+B_{(1)}Z^{B}_{(1)}(z),$ (75a)
$\displaystyle Z^{A}_{(1)}$
$\displaystyle(z)=\text{HeunC}\left(q-\alpha,-\alpha,\delta,\gamma,-\varepsilon;1-z\right),$
(75b) $\displaystyle Z^{B}_{(1)}$ $\displaystyle(z)=(1-z)^{1-\delta}$
$\displaystyle\times$
$\displaystyle\text{HeunC}\left(q^{\prime}-\alpha^{\prime},-\alpha^{\prime},2-\delta,\gamma,-\varepsilon;1-z\right),$
(75c)

with arbitrary complex coefficients $A_{(1)}$ and $B_{(1)}$, and with the
definitions

$\displaystyle q^{\prime}$ $\displaystyle=q-\gamma(1-\delta),$ (76a)
$\displaystyle\alpha^{\prime}$ $\displaystyle=\alpha+\varepsilon(1-\delta).$
(76b)

Asymptotically, the modes of Eq. (75) at the event horizon (${z=1}$) reduce to

$R(z)\to\text{e}^{\frac{\varepsilon}{2}}|z-1|^{-\frac{1-\delta}{2}}\left(A_{(1)}+B_{(1)}(1-z)^{1-\delta}\right),$
(77)

which leads to the matching

$\displaystyle{}^{\text{ext}}A_{(1)}^{+}$ $\displaystyle=0,$
$\displaystyle{}^{\text{ext}}B_{(1)}^{+}$
$\displaystyle=\mathcal{T}^{+}_{\text{ext}}\frac{\text{e}^{i\omega(2r_{+}-r_{-})}}{r_{+}};$
(78a) $\displaystyle{}^{\text{int}}A_{(1)}^{+}$
$\displaystyle=\left(\mathcal{T}^{+}_{\text{int}}+\mathcal{R}^{+}_{\text{int}}\mathcal{R}^{-}_{\text{ext}}\right)\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}},$ $\displaystyle{}^{\text{int}}B_{(1)}^{+}$
$\displaystyle=\mathcal{R}^{+}_{\text{int}}\frac{\text{e}^{i\omega(2r_{+}-r_{-})}}{r_{+}};$
(78b) $\displaystyle{}^{\text{int}}A_{(1)}^{-}$
$\displaystyle=\mathcal{T}^{-}_{\text{int}}\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}},$ $\displaystyle{}^{\text{int}}B_{(1)}^{-}$
$\displaystyle=\left(\mathcal{R}^{-}_{\text{int}}+\mathcal{T}^{-}_{\text{int}}\mathcal{R}^{-}_{\text{ext}}\right)\frac{\text{e}^{i\omega(2r_{+}-r_{-})}}{r_{+}};$
(78c)

for each respective set of modes; i.e., the coefficients from Eq. (75) for
${{}^{\text{ext,int}}R^{\pm}_{\text{ob}}}$ are labeled
${{}^{\text{ext,int}}A^{\pm}_{(1)}}$ and ${{}^{\text{ext,int}}B^{\pm}_{(1)}}$.
Eqs. (78) are strictly only valid for ${z<1}$; for the exterior (${z>1}$), an
additional factor of ${\exp[2\pi\omega r_{+}^{2}/(r_{+}-r_{-})]}$ must be
included in the right-hand side of the equations for each of the $B$
coefficients to account for the lack of absolute values in the trailing factor
of Eq. (77).

At some point $z_{*}$ in the interior (we take ${z_{*}=0.5}$ for simplicity),
both Eqs. (70) and (75) provide regular solutions to the wave Eq. (68). One
can convert between them with the respective linear systems

$\displaystyle
Z^{A,B}_{(1)}(z_{*})=C_{A}^{A,B}Z^{A}_{(0)}(z_{*})+C_{B}^{A,B}Z^{B}_{(0)}(z_{*}),$
(79a) $\displaystyle(Z^{A,B}_{(1)})^{\prime}(z)\big{|}_{z_{*}}$
$\displaystyle=C_{A}^{A,B}(Z^{A}_{(0)})^{\prime}(z)\big{|}_{z_{*}}+C_{B}^{A,B}(Z^{B}_{(0)})^{\prime}(z)\big{|}_{z_{*}}.$
(79b)

The functions $Z^{A}_{(0)}(z_{*})$, $Z^{B}_{(0)}(z_{*})$, $Z^{A}_{(1)}(z_{*})$
and $Z^{B}_{(1)}(z_{*})$ can be computed numerically, and therefore the
constants $C^{A}_{A}$, $C^{B}_{A}$, $C^{A}_{B}$, and $C^{B}_{B}$ can also be
computed. Once these constants are known, the total eigenmodes $Z_{(0)}(z)$
and $Z_{(1)}(z)$ can be matched to solve for each of the backscattering
coefficients:

$\displaystyle A_{(0)}$ $\displaystyle=A_{(1)}C_{A}^{A}+B_{(1)}C_{A}^{B},$
(80a) $\displaystyle B_{(0)}$
$\displaystyle=A_{(1)}C_{B}^{A}+B_{(1)}C_{B}^{B}.$ (80b)

Once the backscattering coefficients connecting ${z=0}$ to ${z=1}$ are known,
a similar process will yield the coefficients connecting ${z=1}$ to
${z=\infty}$. As $z$ approaches infinity, the confluent Heun solutions to Eq.
(68) asymptotically (in a sector) take the form

$\displaystyle
Z_{(\infty)}(z)=A_{(\infty)}Z^{A}_{(\infty)}(z)+B_{(\infty)}Z^{B}_{(\infty)}(z),$
(81a) $\displaystyle Z^{A}_{(\infty)}(z)=z^{-\frac{\alpha}{\varepsilon}},$
(81b) $\displaystyle Z^{B}_{(\infty)}(z)=\text{e}^{-\varepsilon
z}z^{\frac{\alpha}{\varepsilon}-\gamma-\delta},$ (81c)

with arbitrary complex coefficients $A_{(\infty)}$ and $B_{(\infty)}$.
Comparison with the asymptotic forms of Eq. (72) reveals the following matched
values for these coefficients:

$\displaystyle{}^{\text{ext}}A_{(\infty)}^{+}$
$\displaystyle=\mathcal{R}^{+}_{\text{ext}}\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}-r_{-}},$ $\displaystyle{}^{\text{ext}}B_{(\infty)}^{+}$
$\displaystyle=\frac{\text{e}^{i\omega r_{-}}}{r_{+}-r_{-}};$ (82a)
$\displaystyle{}^{\text{int}}A_{(\infty)}^{+}$ $\displaystyle=0,$
$\displaystyle{}^{\text{int}}B_{(\infty)}^{+}$
$\displaystyle=\mathcal{R}^{+}_{\text{int}}\mathcal{T}^{-}_{\text{ext}}\frac{\text{e}^{i\omega
r_{-}}}{r_{+}-r_{-}};$ (82b) $\displaystyle{}^{\text{int}}A_{(\infty)}^{-}$
$\displaystyle=\mathcal{T}^{-}_{\text{int}}\mathcal{T}^{-}_{\text{ext}}\frac{\text{e}^{-i\omega
r_{-}}}{r_{+}-r_{-}},$ $\displaystyle{}^{\text{int}}B_{(\infty)}^{-}$
$\displaystyle=0,$ (82c)

where the coefficient notation is the same as in Eq. (78).

For some sufficiently large radial coordinate ${z=z^{*}}$ (we find
heuristically that ${z^{*}=18/(\omega\sqrt{1-Q^{2}})}$ is more than sufficient
to ensure convergence at machine-level precision), both Eqs. (75) and (81)
satisfy the wave Eq. (68), and so the two sets of solutions can be matched.
One has the system

$\displaystyle
Z^{A,B}_{(1)}(z^{*})=D_{A}^{A,B}Z^{A}_{(\infty)}(z^{*})+D_{B}^{A,B}Z^{B}_{(\infty)}(z^{*}),$
(83a) $\displaystyle(Z^{A,B}_{(1)})^{\prime}(z)\big{|}_{z^{*}}$
$\displaystyle=D_{A}^{A,B}(Z^{A}_{(\infty)})^{\prime}(z)\big{|}_{z^{*}}+D_{B}^{A,B}(Z^{B}_{(\infty)})^{\prime}(z)\big{|}_{z^{*}}$
(83b)

to solve for the constants $D^{A}_{A}$, $D^{B}_{A}$, $D^{A}_{B}$, and
$D^{B}_{B}$, which can then be used to solve for the backscattering
coefficients with the system

$\displaystyle A_{(\infty)}$
$\displaystyle=A_{(1)}D_{A}^{A}+B_{(1)}D_{A}^{B},$ (84a) $\displaystyle
B_{(\infty)}$ $\displaystyle=A_{(1)}D_{B}^{A}+B_{(1)}D_{B}^{B}.$ (84b)

Altogether, the relevant backscattering coefficients can be written as follows
(note that multiple variations to the below equations are possible based on
implicit relations between and among the $C$ and $D$ coefficients):

$\displaystyle\mathcal{T}^{+}_{\text{ext}}$
$\displaystyle=\frac{r_{+}}{r_{+}-r_{-}}\frac{1}{D_{B}^{B}}\
\text{e}^{-2i\omega(r_{+}-r_{-})-\pi\omega/\varkappa_{+}},$ (85a)
$\displaystyle\mathcal{T}^{-}_{\text{ext}}$
$\displaystyle=\frac{r_{+}-r_{-}}{r_{+}}\frac{\widetilde{D}}{D_{B}^{B}},$
(85b) $\displaystyle\mathcal{R}^{+}_{\text{ext}}$
$\displaystyle=\frac{D_{A}^{B}}{D_{B}^{B}}\ \text{e}^{2i\omega r_{-}},$ (85c)
$\displaystyle\mathcal{R}^{-}_{\text{ext}}$
$\displaystyle=-\frac{D_{B}^{A}}{D_{B}^{B}}\ \text{e}^{-2i\omega
r_{+}-\pi\omega/\varkappa_{+}},$ (85d)
$\displaystyle\mathcal{T}^{+}_{\text{int}}$
$\displaystyle=-\frac{r_{+}}{r_{-}}\frac{C_{A}^{B}D_{B}^{B}-C_{A}^{A}D_{B}^{A}\
\text{e}^{-4i\omega r_{+}-\pi\omega/\varkappa_{+}}}{\widetilde{C}D_{B}^{B}},$
(85e) $\displaystyle\mathcal{T}^{-}_{\text{int}}$
$\displaystyle=-\frac{r_{+}}{r_{-}}\frac{C_{A}^{B}}{\widetilde{C}},$ (85f)
$\displaystyle\mathcal{R}^{+}_{\text{int}}$
$\displaystyle=\frac{r_{+}}{r_{-}}\frac{C_{A}^{A}}{\widetilde{C}}\
\text{e}^{-2i\omega r_{+}},$ (85g) $\displaystyle\mathcal{R}^{-}_{\text{int}}$
$\displaystyle=\frac{r_{+}}{r_{-}}\frac{C_{A}^{A}D_{B}^{B}-C_{A}^{B}D_{B}^{A}}{\widetilde{C}D_{B}^{B}}\
\text{e}^{-2i\omega r_{+}},$ (85h)

where

$\displaystyle\widetilde{C}$ $\displaystyle\equiv
C_{A}^{A}C_{B}^{B}-C_{A}^{B}C_{B}^{A},$ (86) $\displaystyle\widetilde{D}$
$\displaystyle\equiv D_{A}^{A}D_{B}^{B}-D_{A}^{B}D_{B}^{A}.$ (87)

The resulting numerical values of the backscattering coefficients are used to
calculate the Hawking spectra of Fig. (7).

## References

  * Hawking [1974] S. W. Hawking, Nature 248, 30 (1974).
  * Hawking [1975] S. W. Hawking, Communications in Mathematical Physics 43, 199 (1975).
  * Parikh and Wilczek [2000] M. K. Parikh and F. Wilczek, Phys. Rev. Lett. 85, 5042 (2000).
  * Dey  _et al._ [2017] R. Dey, S. Liberati, and D. Pranzetti, Physics Letters B 774, 308 (2017).
  * Giddings [2016] S. B. Giddings, Physics Letters B 754, 39 (2016).
  * Barceló  _et al._ [2011] C. Barceló, S. Liberati, S. Sonego, and M. Visser, Journal of High Energy Physics 2011, 1 (2011).
  * Barceló  _et al._ [2011] C. Barceló, S. Liberati, S. Sonego, and M. Visser, Phys. Rev. D 83, 041501 (2011).
  * Barbado  _et al._ [2011] L. C. Barbado, C. Barceló, and L. J. Garay, Classical and Quantum Gravity 28, 125021 (2011).
  * Penrose [1968] R. Penrose, in  _Battelle Rencontres: 1967 lectures in mathematics and physics_ , edited by C. de Witt-Morette and J. A. Wheeler (W. A. Benjamin, New York, 1968) pp. 121–235.
  * Simpson and Penrose [1973] M. Simpson and R. Penrose, Int. J. Theor. Phys. 7, 183 (1973).
  * Poisson and Israel [1990] E. Poisson and W. Israel, Phys. Rev. D 41, 1796 (1990).
  * Barrabes  _et al._ [1990] C. Barrabes, W. Israel, and E. Poisson, Class. Quant. Grav. 7, L273 (1990).
  * Misner and Sharp [1964] C. W. Misner and D. H. Sharp, Phys. Rev. 136, B571 (1964).
  * Gnedin and Gnedin [1993] M. L. Gnedin and N. Y. Gnedin, Class. Quant. Grav. 10, 1083 (1993).
  * Burko [2002] L. M. Burko, Phys. Rev. D 66, 024046 (2002).
  * Burko [2003] L. M. Burko, Phys. Rev. Lett. 90, 121101 (2003).
  * Hamilton and Avelino [2010] A. J. Hamilton and P. P. Avelino, Physics Reports 495, 1 (2010).
  * Hamilton and Polhemus [2011] A. J. S. Hamilton and G. Polhemus, Phys. Rev. D 84, 124055 (2011).
  * Hamilton [2011a] A. J. S. Hamilton, Phys. Rev. D 84, 124056 (2011a).
  * Hamilton [2011b] A. J. S. Hamilton, Phys. Rev. D 84, 124057 (2011b).
  * Hamilton [2017] A. J. S. Hamilton, Phys. Rev. D 96, 084041 (2017).
  * Chesler  _et al._ [2020] P. M. Chesler, R. Narayan, and E. Curiel, Class. Quant. Grav. 37, 025009 (2020).
  * McMaken and Hamilton [2021] T. McMaken and A. J. S. Hamilton, Phys. Rev. D 103, 084014 (2021).
  * Ori [1991] A. Ori, Phys. Rev. Lett. 67, 789 (1991).
  * Ori [1992] A. Ori, Phys. Rev. Lett. 68, 2117 (1992).
  * Brady and Chambers [1995] P. R. Brady and C. M. Chambers, Phys. Rev. D 51, 4177 (1995).
  * Ori and Flanagan [1996] A. Ori and E. E. Flanagan, Phys. Rev. D 53, R1754 (1996).
  * Ori [1998] A. Ori, Phys. Rev. D 57, 4745 (1998).
  * Ori [1999] A. Ori, Phys. Rev. Lett. 83, 5423 (1999).
  * Dafermos [2005] M. Dafermos, Communications on Pure and Applied Mathematics 58, 445 (2005).
  * Novikov and Starobinskii [1980] I. D. Novikov and A. A. Starobinskii, Sov. Phys. - JETP (Engl. Transl.) 51 (1980).
  * Herman and Hiscock [1994] R. Herman and W. A. Hiscock, Phys. Rev. D 49, 3946 (1994).
  * Sorkin and Piran [2001a] E. Sorkin and T. Piran, Phys. Rev. D 63, 084006 (2001a).
  * Sorkin and Piran [2001b] E. Sorkin and T. Piran, Phys. Rev. D 63, 124024 (2001b).
  * Hong  _et al._ [2010] S. E. Hong, D.-I. Hwang, E. D. Stewart, and D.-H. Yeom, Classical and Quantum Gravity 27, 045014 (2010).
  * Zilberman  _et al._ [2020] N. Zilberman, A. Levi, and A. Ori, Phys. Rev. Lett. 124, 171302 (2020).
  * Hollands  _et al._ [2020] S. Hollands, R. M. Wald, and J. Zahn, Class. Quant. Grav. 37, 115009 (2020).
  * Lanir  _et al._ [2019] A. Lanir, A. Ori, N. Zilberman, O. Sela, A. Maline, and A. Levi, Phys. Rev. D 99, 061502 (2019).
  * Taylor [2020] P. Taylor, Class. Quant. Grav. 37, 045004 (2020).
  * Barceló  _et al._ [2021] C. Barceló, V. Boyanov, R. Carballo-Rubio, and L. J. Garay, Class. Quant. Grav. 38, 125003 (2021).
  * Arrechea  _et al._ [2021a] J. Arrechea, C. Barceló, R. Carballo-Rubio, and L. J. Garay, Class. Quant. Grav. 38, 115014 (2021a).
  * Arrechea  _et al._ [2021b] J. Arrechea, C. Barceló, V. Boyanov, and L. J. Garay, Universe 7 (2021b).
  * Zilberman and Ori [2021] N. Zilberman and A. Ori, Phys. Rev. D 104, 024066 (2021).
  * Klein  _et al._ [2021] C. Klein, J. Zahn, and S. Hollands, Phys. Rev. Lett. 127, 231301 (2021).
  * Hamilton [2018] A. J. S. Hamilton, General Relativity and Gravitation 50, 50 (2018).
  * Crispino  _et al._ [2016] L. C. B. Crispino, A. Higuchi, L. A. Oliveira, and E. S. de Oliveira, The European Physical Journal C 76, 168 (2016).
  * Ong and Good [2020] Y. C. Ong and M. R. R. Good, Phys. Rev. Research 2, 033322 (2020).
  * Curir [1979] A. Curir, Il Nuovo Cimento B (1971-1996) 51, 262 (1979).
  * Cvetič  _et al._ [2018] M. Cvetič, G. W. Gibbons, H. Lü, and C. N. Pope, Phys. Rev. D 98, 106015 (2018).
  * Peltola and Mäkelä [2006] A. Peltola and J. Mäkelä, Internat. J. Modern Phys. D 15, 817 (2006).
  * Park [2008] M.-I. Park, Physics Letters B 663, 259 (2008).
  * Tiandho [2017] Y. Tiandho, J. Math. Fundam. Sci. 49, 294 (2017).
  * Unruh [1976] W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * Hamilton and Lisle [2008] A. J. S. Hamilton and J. P. Lisle, American Journal of Physics 76, 519 (2008).
  * Brynjolfsson and Thorlacius [2008] E. J. Brynjolfsson and L. Thorlacius, Journal of High Energy Physics 2008, 066 (2008).
  * Loranz and Hiscock [1997] D. J. Loranz and W. A. Hiscock, Phys. Rev. D 55, 3893 (1997).
  * Hiscock [1980] W. A. Hiscock, Phys. Rev. D 21, 2057 (1980).
  * Barbado  _et al._ [2016] L. C. Barbado, C. Barceló, L. J. Garay, and G. Jannes, Phys. Rev. D 94, 064004 (2016).
  * Birrell and Davies [1982] N. D. Birrell and P. C. W. Davies,  _Quantum Fields in Curved Space_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 1982).
  * Anderson  _et al._ [2013] P. R. Anderson, R. Balbinot, A. Fabbri, and R. Parentani, Phys. Rev. D 87, 124018 (2013).
  * Balbinot  _et al._ [2019] R. Balbinot, A. Fabbri, R. A. Dudley, and P. R. Anderson, Phys. Rev. D 100, 105021 (2019).
  * Starobinskiǐ and Churilov [1973] A. A. Starobinskiǐ and S. M. Churilov, Zh. Eksp. Teor. Fiz. 65 (1973).
  * Lanir  _et al._ [2018] A. Lanir, A. Levi, A. Ori, and O. Sela, Phys. Rev. D 97, 024033 (2018).
  * Page [1976] D. N. Page, Phys. Rev. D 13, 198 (1976).
  * Vieira  _et al._ [2022] H. S. Vieira, V. B. Bezerra, and C. R. Muniz, The European Physical Journal C 82, 932 (2022).
  * Kazakov [2006] A. Y. Kazakov, Journal of Physics A: Mathematical and General 39, 2339 (2006).
  * Motygin [2018] O. V. Motygin, in  _2018 Days on Diffraction (DD)_ (2018) pp. 223–229.


11institutetext: Department of Electronics and Communication Engineering,
Heritage Institute of Technology, Kolkata, West Bengal
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Longitudinal Volumetric Study for the Progression of Alzheimer’s Disease
from Structural MR Images

Prayas Sanyal  Srinjay Mukherjee  Arkapravo Das  Anindya Sen

###### Abstract

Alzheimer’s Disease (AD) is primarily an irreversible neurodegenerative
disorder affecting millions of individuals today. The prognosis of the disease
solely depends on treating symptoms as they arise and proper caregiving, as
there are no current medical preventative treatments. For this purpose, early
detection of the disease at its most premature state is of paramount
importance. This work aims to survey imaging biomarkers corresponding to the
progression of Alzheimer’s Disease (AD). A longitudinal study of structural MR
images was performed for given temporal test subjects selected randomly from
the Alzheimer’s Disease Neuroimaging Initiative (ADNI) database. The pipeline
implemented includes modern pre-processing techniques such as spatial image
registration, skull stripping, and inhomogeneity correction. The temporal data
across multiple visits spanning several years helped identify the structural
change in the form of volumes of cerebrospinal fluid (CSF), grey matter (GM),
and white matter (WM) as the patients progressed further into the disease.
Tissue classes are segmented using an unsupervised learning approach using
intensity histogram information. The segmented features thus extracted provide
insights such as atrophy, increase or intolerable shifting of GM, WM and CSF
and should help in future research for automated analysis of Alzheimer’s
detection with clinical domain explainability.

###### Keywords:

Alzheimer’s Disease, Feature Extraction, MRI, Longitudinal Analysis, ADNI

##  1 Introduction

Reports show that 8.8 million Indians above the age of 60 live with dementia
(an umbrella term for Alzheimer’s Disease), amounting to approximately 7.4% of
the total population. It is projected that the number of dementia afflicted
may increase to 16.9 million by 2036 with the rapid rise of correlated risk
factors such as diabetes and hypertension, both under-treated in rural
settings [1, 2]. Common AD symptoms include memory impairment, loss of
cognitive function, depression and paranoia. Even though these symptoms
inevitably worsen over time, those providing care feel more in control with an
early diagnosis and seek assistance from support groups along with others who
have experienced similar misfortunes [3]. A review published in 2013 states
that the annual overhead cost of looking after a person with Alzheimer’s in
India can reach heights of 2,02,450 in urban areas and 66,025 in rural areas,
including intangible costs such as loss of productivity, medication cost,
consultation and hospitalization [4].

German physicist Alois Alzheimer laid the foundations of most of the modern
understanding of the underlying causes of this disorder when he noticed the
common biomarkers through histological techniques on the brains of his
diseased patients [5]. Recent technologies exist, such as MRI, which helps in
visualising atrophy or change in tissue volume of such biomarkers over a
period of time, or fMRI, which can identify high-resolution activation of
different brain regions during various cognitive tasks [6]. Fig. 1 shows the
example of an AD patient’s raw unprocessed brain slice in different planes
obtained from MR imaging in a magnetization-prepared rapid gradient-echo (MP-
RAGE) sequence.

|  |
---|---|---

(a)

(b)

(c)

Figure 1: MP-RAGE T1 Weighted MR slices: (a) Axial plane (b) Sagittal plane
(c) Coronal plane

###  1.1 Related Works

Multiple large-scale longitudinal studies have been carried out to investigate
the causal connection between various AD indicators and the development of
cognitive decline and dementia (e.g., ADNI, Australian Imaging) [7, 8]. Due to
the dearth of publicly available datasets corresponding to the Indian
populace, ADNI database was extensively used for our experimentation. Over the
past two decades, several techniques have been applied to measuring brain
atrophy in structural MRI, such as the boundary shift integral technique [9,
10], which uses a linear registration to align base and follow-up scans to
track the change in brain boundary; or tensor-based morphometry (TBM) [11]
which utilises non-linear image registration to map the volume change between
the baseline and follow up images. Numerous machine-learning feature
extraction techniques have also been studied for distinctive biomarker
detection in the advent of AD, but most lack reproducibility and
generalizability in the current scenario [12, 13, 14].

##  2 Methodology

We put forward a methodology to extract imaging features from MP-RAGE MRI
images to include medical domain explainability for future group pattern
analysis works. For each patient, a standard pipeline was implemented in the
open-source R software (integrated with fslr [15]). Fig 2 represents a
flowchart to elucidate the steps such as database formation, pre-processing of
the images, and feature extraction. We have only analysed the volumetric
changes of three particular tissues, namely, cerebrospinal fluid (CSF), grey
matter (GM) and white matter (WM), over the course of multiple visits for a
patient diagnosed with AD. These volumetric data thus obtained are considered
to be the features of Alzheimer’s Disease.

Figure 2: Flowchart for Structural Analysis in MP-RAGE MR images

###  2.1 Database Used

A representative dataset was created by selecting MP-RAGE structural MRI scans
of AD patients at random after gaining access to the Alzheimer’s Disease
Neuroimaging Initiative (ADNI - adni.loni.usc.edu), which was launched in 2003
and led by Principal Investigator Michael W. Weiner, MD. Patients were chosen
such that they had at least three visits to track the change in the features
of their disease progression.

###  2.2 Pre-processing

This section explains in detail the pre-processing paradigm implemented to
align all the different visits. A wrapper script thus created with each of the
steps enables us to easily process incoming files and perform the collection
of transformations to further analyze our data.

####  2.2.1 File Preparation

The primary acquisition plane for each of the scans was sagittal, and the MP-
RAGE scans downloaded from the ADNI website were in Digital Imaging and
Communications in Medicine (DICOM) format. The DICOM format files (256
separated files for each slice) were first converted into a NIfTI stack (1
single file) to represent each visit for easier processing and visualization.
In DICOM image files, the data comprises 2D layers, whereas NIfTI
(Neuroimaging Informatics Technology Initiative) is a file format where images
and other data are stored in a 3D structure. To access the DICOM files and
convert them into a NIfTI file format, the $oro.dicom$ and $oro.nifti$
packages in R are used, which was introduced by Brandon Witcher et al. [16].

####  2.2.2 Inhomogeneity Correction

A significant issue faced while analyzing structural properties of the brain
is the non-uniformity of intensity in structural MR images. Any low-frequency
intensity which can be non-uniformly present in the data is referred to as
inhomogeneity or bias field. A standard inhomogeneity correction algorithm
known as the N4ITK or Improved Bias Correction [17] was used to remove low-
frequency intensities. The image formation model of N3 is given by the
following equation:

$v(x)=u(x)f(x)+n(x)$

where $v$ is the input image, $u$ is the uncorrupted actual image, $x$ is any
location in the image, $f$ is the bias field, and $n$ is noise(which is
assumed to be independent and gaussian). As the data is log-transformed and
assumes a noise-free scenario ($n(x)=0$), the previous equation becomes

$log(v(x))=log(u(x))+log(f(x))$

The N4 algorithm then uses a robust B-spline approximation algorithm of the
bias field and iterates until a convergence criterion is met where the
modified bias field is approximately similar to the last performed iteration.
After the updated transformation, the image is outputted back into the
original file using inverse log transform. To visualise the effects of bias
field correction Fig 3 shows a sample patient scan in all three planes by
applying N4 bias correction using ANTsR library in R.

(a)

(b)

(c)

Figure 3: Inhomogeneity correction: (a) Source Scan (b) Bias Field Corrected
Scan (c) Bias field corrected scan subtracted from source scan

Even though visually, there might not be a stark difference in the corrected
scan, white and grey matter intensities are uniform in distribution compared
to the source image. Fig 3(c) helps to visualize the difference where the blue
gradient represents higher differences. It becomes imperative how there was a
slowly varying field existing in the source image, which was removed in the
corrected version.

####  2.2.3 Skull Stripping

This pre-processing technique removes the skull voxels from the T1-weighted
MRI scan. To do this operation on the N4 bias-corrected image, FSL software
[18] was used, which is integrated into R using the fslR package. Fslr
performs required operations on NIfTI image objects using FSL commands and
returns them as R objects [15]. The function $fslr::fslbet$ is readily
available, which takes in the N4 corrected NIfTI image and calls the FSL BET
function to perform skull stripping. Brain Extraction Tool (BET) [19] uses a
deformable model that evolves to fit into the brain tissue surface locally
adaptive model forces. To form a brain tissue mask which can be overlayed into
the input scan, an NIfTI array is created and filled with binary ones using
dimensions information from the header file. Secondly, the areas of the skull-
stripped image containing the brain tissue were determined, or every pixel
with a positive intensity value, because the skull stripping algorithm BET
sets everything that is not brain to binary zeroes. A centre of gravity (CoG)
of the extracted tissue was determined during the first iteration, and the
algorithm was repeated by feeding the CoG as additional information to improve
the skull stripping paradigm as one single iteration was observed to contain
deformities. Bias correction, explained in the previous subsection, is
essential for this step as the skull stripping algorithm heavily depends on
the tissue intensities and a low variation in the field can hinder the output.
The skull-stripped image thus obtained is shown in Fig 4.

|
---|---

(a)

(b)

Figure 4: Skull Stripping: (a) Brain tissue mask overlay (b) Extracted brain
tissue

####  2.2.4 Image Registration to Base Scan

The AD-affected patient scans selected with three or four visits need to be
registered to make image locations such as voxels, and tissues have similar
interpretations and align them spatially to compare changes. Each subsequent
visit of each patient was co-registered to their source visit after performing
inhomogeneity correction and skull stripping. Very few degrees of freedom are
needed to register follow-up scans to the base scan as they belong to the same
modality, so a rigid or affine registration is sufficient. Rigid registration
is a linear registration technique containing six degrees of freedom and is
represented by the equation: $T_{rigid}(v)=R_{v}+t$, where t is a translation
vector $t=(t_{x},t_{y},t_{z})$ and R is a rotation matrix. All of the above
operations were performed using the ANTsR package available in R.

###  2.3 Tissue Segmentation and Volume Extraction

After applying the necessary pre-processing steps described in the former
subsections on the MR scans, a three-class tissue segmentation is performed
using the FAST algorithm on the FSLR package integrated into R [15, 20]. The
FAST algorithm, which uses the pipeline as described by Zhang et al., puts
forward a hidden Markov random field model (HMRF) and an expectation
maximization algorithm to perform probabilistic tissue segmentation. The
intensity histogram is important for observing the working of the algorithm,
so a sample intensity histogram is shown in Fig 5 of an Alzheimer’s affected
patient’s MRI scan after preprocessing(no bias field, noise removed, skull
stripped).

Figure 5: Intensity histogram of a sample AD patient scan after pre-processing

If the voxel frequency vs. intensity plot is observed, the three different
tissue classes are apparent as three Gaussian peaks. The lowest intensity peak
corresponds to Cerebrospinal Fluid (CSF), whereas the subsequent peaks
determine the intensity peaks of grey matter (GM) and white matter (WM). A
probability distribution for each tissue class is assigned such that
$\mathcal{P}$(CSF) is near zero intensity levels; $\mathcal{P}$(GM) is of low
to moderate levels of intensity; $\mathcal{P}$(WM) is of moderate to high
levels of intensity. Spatial neighbourhood information or the HMRF model is
used to segment each tissue class as simple classifiers like K-means do not
use spatial entropy information, and we need a proper balance. The log
probability model can be represented as $logp(intensity)+\beta logp(MRF)$,
where $\beta$ is user-adjustable. Output files of the FAST algorithm provide
three separate probability maps for corresponding tissue classes: $pve_{0}$,
$pve_{1}$, and $pve_{2}$ for CSF, GM and WM, respectively. The $ortho2$
function in R is used to overlay the probability maps on the pre-processed T1
weighted image, as shown in Fig 6.

(a)

(b)

(c)

Figure 6: Segmentation of tissue classes of sample AD patient slice:
(a)Cerebrospinal fluid (b) Grey Matter (c) White Matter

The changes in volumes of each patient for the tissue classes across three or
four visits are observed and reported. These changes can be considered as
relevant imaging features to understand the prognosis of Alzheimer’s disease.

##  3 Experimental Results

The above pipeline was iterated for each of the visits of the patients chosen
from the ADNI database. $0.4$ was empirically chosen as the threshold for the
output tissue probability maps after multiple experiments to calculate the
volumes of each tissue class for the respective scans to extract features.
First, the voxel dimension was calculated with a function called $voxdim$
available in $oro.nifti$ package. Then the number of voxels above the chosen
threshold is calculated corresponding to the probability maps for each tissue
class, given by $t_{vox}$. Finally, the volume is calculated in $mL$ with the
equation $Volume=vox_{dim}\times(t_{vox}\div 1000)$, where $vox_{dim}$ is the
voxel dimensions in $(ml)$ and the respective tissue voxels are $t_{vox}$
which is normalized to $1000$. The changes in volume in each tissue class were
plotted in Fig 8 with each line depicting a patient after max scaling the
data. Twenty-six patients’ volumetric data were observed, but after several
permutations, fifteen of them were chosen to be plotted to show a majority
trend (either increase in volume or decrease). The data is made available upon
request to the authors for privacy concerns of ADNI guidelines.

(a) Grey matter volume changes

(b) Cerebrospinal fluid volume changes

Figure 7: Line plot for volume changes of each tissue class (a) White matter
volume changes Figure 8: Line plot for volume changes of each tissue class
(cont.)

The mean volume and standard deviation for each chosen visit corresponding to
each tissue class were calculated and reported in Table 1.

Table 1: Volumes of Each Tissue Class (Mean $\pm$ Standard Deviation) Volume ($mL$) | Visit 1 | Visit 2 | Visit 3 | Visit 4
---|---|---|---|---
Grey Matter | 551.34 $\pm$ 61.56 | 538.01 $\pm$ 54.35 | 526.23 $\pm$ 54.01 | 497.42 $\pm$ 79.36
Cerebrospinal Fluid | 383.67 $\pm$ 61.02 | 382.39 $\pm$ 64.50 | 391.09 $\pm$ 63.60 | 381.29 $\pm$ 42.86
White Matter | 554.27 $\pm$ 69.86 | 553.63 $\pm$ 58.93 | 554.81 $\pm$ 58.27 | 549.02 $\pm$ 83.78

##  4 Discussion

We have created an easy-to-reproduce preprocessing pipeline and a wrapper
script in R software to determine the volume from structural MRI scans of AD-
affected patients for large-scale longitudinal analysis. Due to limited access
to computational resources and an absence of Alzheimer’s datasets
corresponding to the Indian landscape, we have considered MR scans of twenty-
six patients from the ADNI database and compared their volume changes across
three different tissue classes, namely cerebrospinal fluid (CSF), white matter
(WM) and grey matter (GM). Of the twenty-six patients, fifteen were chosen,
following a common pattern and their corresponding volumes were plotted across
the three or four visits. The few inconsistencies can be attributed to several
factors, such as responsiveness to symptom treatment, inaccurate MR scan
acquisition, or incorrect segmentation due to spurious noise and other
correlated abnormalities. In general, dementia is considered a grey matter
degenerative disease as the common symptoms are congruent to a loss of grey
matter volume [21]. Alzheimer’s affected patients show a massive loss in grey
matter volume when compared to normal control counterparts. Our experimented
volumes are consistent with clinical findings as we observe a general decline
in GM volume, as shown in Fig 7(a) over a patient’s prognosis across several
visits. Similarly, when Fig 7(b) is observed, it is imperative that an
intolerable CSF shifting should also be considered as a common biomarker for
Alzheimer’s Disease [22]. Fig 8(a) shows the white matter volume changes in AD
patients, and even though a majority trend could not be identified from the
line plot, specific abnormalities such as hyperintensities and lesion
formation have been analyzed for their correlation with Alzheimer’s in
multiple studies [23, 24].

The mean and standard deviation of each volume calculated and reported in
Table 1 shed light on the behaviour of volume changes. As discussed above,
grey matter is observed to follow a general decline across the four visits and
can be heavily correlated with Alzheimer’s Disease. However, white matter and
cerebrospinal fluid volumes stay consistent with marginal change in the mean
volume but exhibit rapid shifts in the line plots.

This presented work can be considered as a first step for future large-scale
analysis of the exhaustive MR scans of the ADNI dataset with more
experimentation on different thresholds and more accurate pre-processing and
segmentation techniques. Instead of black-box automated analysis from raw MR
scans or intensive annotation methods, this added step of calculating
additional imaging features related to the disease adds a layer of
explainability with medical domain knowledge.

##  5 Acknowledgements

All the scans collected for analysis in preparation for this paper were
accessed from the Laboratory of Neuro-imaging at the University of South
California, which distributes ADNI data. ADNI is funded by the National
Institute of Aging, the National Institute of Biomedical Imaging and
Bioengineering, and through generous contributions from the following: AbbVie,
Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon
Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir,
Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company;
EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech,
Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy
Research & Development, LLC.; Johnson & Johnson Pharmaceutical Research &
Development LLC.; Lumosity; Lundbeck; Merck & Co., Inc.; Meso Scale
Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis
Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda
Pharmaceutical Company; and Transition Therapeutics. The study also deeply
benefits from the open-source course on Introduction to Neurohacking in R
offered by Johns Hopkins University.

## References

  * [1] Lee, J, Meijer, E, Langa, KM, et al. Prevalence of dementia in India: National and state estimates from a nationwide study. Alzheimer’s Dement. 2023; 19: 2898–2912. 10.1002/alz.12928
  * [2] Ravindranath, V., and Sundarakumar, J. S. (2021). Changing demography and the challenge of dementia in India. Nature Reviews Neurology, 17(12), 747-758. 10.1038/s41582-021-00565-x
  * [3] Kumar A, Sidhu J, Goyal A, et al. Alzheimer Disease. [Updated 2022 Jun 5]. In: StatPearls [Internet]. Treasure Island (FL): StatPearls Publishing
  * [4] Rao, G., and Bharath, S. (2013). Cost of dementia care in India: Delusion or reality? Indian Journal of Public Health, 57(2), 71. 10.4103/0019-557X.114986
  * [5] Stelzmann RA, Schnitzlein N, Murtagh FR : An English translation of Alzheimer’s paper, “über eine eigenartige Erkankung der Hirnrinde” Clinical Anatomy (1995)
  * [6] Ogawa S, Lee TM, Kay AR, et al. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci U S A. 1990; 87:9868.
  * [7] Bondi, M. W., Edmonds, E. C., and Salmon, D. P. (2017). Alzheimer’s Disease: Past, Present, and Future. Journal of the International Neuropsychological Society: JINS, 23(9-10), 818. 10.1017/S135561771700100X
  * [8] Jack, C. R., Knopman, D. S., Jagust, W. J., Shaw, L. M., Aisen, P. S., Weiner, M. W., … Trojanowski, J. Q. (2010). Hypothetical model of dynamic biomarkers of the Alzheimer’s pathological cascade. The Lancet Neurology, 9(1), 119–128. 10.1016/s1474-4422(09)70299-6
  * [9] Leung K.K., Ridgway G.R., Ourselin S., Fox N.C. Consistent multi-time-point brain atrophy estimation from the boundary shift integral. Neuroimage. 2012; 59:3995–4005. 10.1016/j.neuroimage.2011.10.068
  * [10] Freeborough P.A., Fox N.C. The boundary shift integral: an accurate and robust measure of cerebral volume changes from registered repeat MRI. IEEE Trans. Med. Imaging. 1997; 16:623–629 10.1109/42.640753
  * [11] Hua X., Hibar D.P., Ching C.R.K., Boyle C.P., Rajagopalan P., Gutman B.A., Leow A.D., Toga A.W., Jack C.R., Harvey D., Weiner M.W., Thompson P.M: Unbiased tensor-based morphometry (TBM): Improved robustness and sample size estimates for Alzheimer’s disease clinical trials. Neuroimage. 2013; 66:648–661 10.1016/j.neuroimage.2012.10.086
  * [12] Pan, D., Zeng, A., Jia, L., Huang, Y., Frizzell, T., and Song, X. (2020). Early Detection of Alzheimer’s Disease Using Magnetic Resonance Imaging: A Novel Approach Combining Convolutional Neural Networks and Ensemble Learning. Frontiers in Neuroscience, 14. 10.3389/fnins.2020.00259
  * [13] Christian S., Antonio C., Petronilla B., Gilardi M. C., Aldo Q., Isabella C.: Magnetic resonance imaging biomarkers for the early diagnosis of Alzheimer’s disease: a machine learning approach. Frontiers in Neuroscience. 2015 9:307 10.3389/fnins.2015.00307
  * [14] Moradi E, Pepe A, Gaser C, Huttunen H, Tohka J; Alzheimer’s Disease Neuroimaging Initiative. Machine learning framework for early MRI-based Alzheimer’s conversion prediction in MCI subjects. Neuroimage. 2015 Jan 1; 104:398-412. 10.1016/j.neuroimage.2014.10.002
  * [15] Muschelli, John & Sweeney, Elizabeth & Lindquist, Martin & Crainiceanu, Ciprian. (2015). Fslr: Connecting the FSL software with R. R Journal. 7. 163-175. 10.32614/RJ-2015-013.
  * [16] Whitcher B, Schmid VJ, Thornton A (2011). ”Working with the DICOM and NIfTI Data Standards in R.” Journal of Statistical Software, 44(6), 1–28. 10.18637/jss.v044.i06.
  * [17] Tustison, N. J., Avants, B. B., Cook, P. A., Zheng, Y., Egan, A., Yushkevich, P. A., & Gee, J. C. (2010). N4ITK: Improved N3 Bias Correction. IEEE Transactions on Medical Imaging, 29(6), 1310. 10.1109/TMI.2010.2046908
  * [18] Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM (2012). “FSL.” Neuroimage, 62(2), 782–790.
  * [19] Smith, S. M. (2002). Fast robust automated brain extraction. Human Brain Mapping, 17(3), 143–155. 10.1002/hbm.10062
  * [20] Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm. IEEE Trans Med Imaging. 2001 Jan;20(1):45-57. 10.1109/42.906424
  * [21] Wu, Zhanxiong, Yun Peng, Ming Hong, and Yingchun Zhang. ”Gray Matter Deterioration Pattern During Alzheimer’s Disease Progression: A Regions-of-Interest Based Surface Morphometry Study.” Frontiers in Aging Neuroscience 13, (2021): 593898. 10.3389/fnagi.2021.593898.
  * [22] Paterson, R.W., Slattery, C.F., Poole, T. et al. Cerebrospinal fluid in the differential diagnosis of Alzheimer’s disease: clinical utility of an extended panel of biomarkers in a specialist cognitive clinic. Alz Res Therapy 10, 32 (2018). 10.1186/s13195-018-0361-3
  * [23] Kao, Hui, et al. ”White Matter Changes in Patients with Alzheimer’S Disease and Associated Factors.” Journal of Clinical Medicine, vol. 8, no. 2, 2019, 10.3390/jcm8020167
  * [24] Gootjes L., Teipel S., Zebuhr Y., Schwarz R., Leinsinger G., Scheltens P., Möller H.-J., Hampel H. Regional distribution of white matter hyperintensities in vascular dementia, alzheimer’s disease and healthy aging. Dement. Geriatr. Cogn. Disord. 2004;18:180–188. 10.1159/000079199