# Near-horizon microstructure and superradiant instability of black holes

Rong-Zhen Guo<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical
Physics, Institute of Theoretical Physics, Chinese Academy of Sciences,
Beijing 100190, China  School of Physical Sciences, University of Chinese
Academy of Sciences, No. 19A Yuquan Road, Beijing 100049, China  Chen Yuan
<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  Qing-Guo Huang  corresponding author:
<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  School of Fundamental Physics and
Mathematical Sciences Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China  Center for Gravitation and Cosmology, College of Physical
Science and Technology, Yangzhou University, Yangzhou 225009, China

###### Abstract

Ultralight bosons, as important candidates of dark matter, can condense around
spinning black holes (BHs) to form long-lived “boson clouds” due to
superradiance instability. The boson-BH system can be observed through
gravitational wave detection and may become a new window to find traces of
ultralight bosons. In this letter we explore the effects on the superradiant
instability of BHs from the near-horizon microstructure. By introducing the
reflection parameter near a BH horizon, we derived analytical results on the
corrections to both energy levels of bosonic cloud and its characteristic
frequencies of superradiance instability. Our results imply that the evolution
of a boson-BH system and gravitational waves it emits would be influenced by
the near-horizon physics of a BH.

###### pacs:

???

Introduction.  Ultralight bosonic particles, which are predicted in various
beyond Standard Model scenarios Arvanitaki  _et al._ (2010a); Essig  _et al._
(2013); Irastorza and Redondo (2018); Goodsell  _et al._ (2009); Jaeckel and
Ringwald (2010); Graham  _et al._ (2016); Agrawal  _et al._ (2020), could be a
promising dark matter (DM) candidate Preskill  _et al._ (1983); Abbott and
Sikivie (1983); Dine and Fischler (1983); Svrcek and Witten (2006); Arvanitaki
_et al._ (2010a); Arvanitaki and Dubovsky (2011); Essig  _et al._ (2013);
Brito  _et al._ (2015a); Marsh (2016); Hui  _et al._ (2017); Annulli  _et al._
(2020); Chadha-Day  _et al._ (2021). Given that detecting ultralight particles
is challenging in particle physics, using gravitational waves (GWs) to probe
ultralight bosons has aroused much attention recently. The main mechanism
responsible for this method is the superradiant instability in a massive boson
and spinning black hole (BH) system Damour  _et al._ (1976); Zouros and
Eardley (1979); Detweiler (1980); Dolan (2007); Arvanitaki  _et al._ (2010b);
Shlapentokh-Rothman (2014); Pani  _et al._ (2012a, b); Witek  _et al._ (2013);
Brito  _et al._ (2013); Endlich and Penco (2017); East (2017); East and
Pretorius (2017); Baryakhtar  _et al._ (2017a); Cardoso  _et al._ (2018); East
(2018); Frolov  _et al._ (2018); Dolan (2018); Baumann  _et al._ (2019a);
Brito  _et al._ (2020); Tsukada  _et al._ (2021); Yuan  _et al._ (2021). When
the Compton wavelength of bosons is about the order of the BH horizon, the
bosons would extract energy and angular momentum from the BH, forming a
macroscopic rotating boson cloud, which dissipates its energy through the
emission of nearly monochromatic GWs on a very lone timescale. Therefore, the
boson-BH system can act as a continuous GW source and could be detected by GW
observations Arvanitaki  _et al._ (2015, 2017); Baryakhtar  _et al._ (2017b);
Brito  _et al._ (2017a, b); Isi  _et al._ (2019); Ghosh  _et al._ (2019);
Palomba  _et al._ (2019); Sun  _et al._ (2020); Zhu  _et al._ (2020); Brito
_et al._ (2020); Ng  _et al._ (2020); Tsukada  _et al._ (2021); Yuan  _et al._
(2021). Since ultralight bosons such as axions or string axiverse could be
important candidates of dark matter and it is difficult to find its evidence
by traditional collision experiments, superradiant instabilities of
astrophysical BHs have been used to constrain its parameters Kodama and
Yoshino (2012); Brito  _et al._ (2015b); Arvanitaki and Dubovsky (2011); Brito
_et al._ (2017a, b); Tsukada  _et al._ (2019); Zhu  _et al._ (2020); Tsukada
_et al._ (2021).

BH in general relativity could be seen as a special region in the spacetime
surrounded by event horizon, which behaves as a smooth one-way membrane Thorne
_et al._ (1986). Because of Penrose-Hawking singularity theorems Hawking and
Ellis (2011), there exists a physically ill-defined singularity inside its
event horizon. However, this picture is challanged by two aspects. One aspect
is the existence of singularity for the belief that there should be no such
pathological point in the real physical world. To avoid the formation of
singularity, the idea that BH is not the end state of gravitational collapse
arises Barcelo  _et al._ (2008). For example, the end state might be a
gravastar Morris  _et al._ (1988), boson star Schunck and Mielke (2003),
wormhole Morris  _et al._ (1988) or other exotic compact object (ECO) Holdom
and Ren (2017). Another aspect stems from the comprehensive consideration of
both gravity and quantum effects. Several inferences of BH thermodynamics like
Bekenstein-Hawking entropy imply the discrete degrees of freedom living on
event horizon, which lead to the BH area quantization (named Bekenstein-
Mukhanov anstaz) Bekenstein (1974); Bekenstein and Mukhanov (1995).
Furthermore, several candidates of quantum gravity predict different quantum-
corrected BH models (e.g. fuzzball Skenderis and Taylor (2008) in string
theory). Moreover, for solving the BH information paradox Polchinski (2017),
firewall Almheiri  _et al._ (2013) might exist near event horizon. All of the
above would significantly change the physics near the event horizon.

Despite the existence of colorful quantum-corrected BH models or ECOs, these
modifications to classical BHs can be described by a phenomenological model
that introduces reflectivity $\mathcal{R}(\omega)$ near the horizon of a BH.
It is widely used in the study of echoes of gravitational waves Cardoso  _et
al._ (2016, 2019); Mark  _et al._ (2017); Price and Khanna (2017); Fang  _et
al._ (2021). Cardoso et al. Cardoso  _et al._ (2019) has proved BH echoes from
BH area quatization. It has been also used to test possible Lorentz violation
near event horizon Oshita and Afshordi (2019). This effective description
could be understood from two different perspectives. Mathematically, such a
boundary condition is essentially the linearization of much general boundary
conditions Burgess  _et al._ (2018). And from the physical point of view,
because the radial equation is a Schrodinger-like equation, the reflection
parameters $\mathcal{R}(\omega)$ here can be considered to mark the nonzero
probability of particle reflection near the black hole horizon. In addition,
Burgess et.al Rummel and Burgess (2019); Burgess  _et al._ (2018) point that
$\mathcal{R}(\omega)$ is related to the couplings of effective field theory
living in the near-horizon region, thus it could be used to constrain
fundamental physics in a more direct way.

Since superradiance instability could be solved as a boundary value problem in
the framework of BH perturbation theory Brito  _et al._ (2015c), it implies
that superradiance instability could be used to detect possible effects of the
near horizon microstructure, thus provide clues of quantum gravity. Here we
study a scalar field minically coupled to gravity on the Kerr background with
the modified boundary condition given near event horizon and derive analytical
results of the corrections to the bound state frequencies of bosonic clouds.
In this letter we take the units $c=G=1$.

Superradiant instability of black holes.  We consider a scalar field $\psi$
with mass $\mu$ in the vicinity of a Kerr black hole whose metric in standard
Boyer-Lindquist coordinates $\left\\{t,r,\theta,\phi\right\\}$ takes the
following form

$\displaystyle
ds^{2}=-\left(1-\frac{2Mr}{\Sigma}\right)dt^{2}-\frac{2Mar\sin^{2}\theta}{\Sigma}dtd\phi$
$\displaystyle+\frac{\Sigma}{\Delta}dr^{2}+\Sigma
d\theta^{2}+\frac{A\sin^{2}\theta}{\Sigma}d\phi^{2},$ (1)

where $M$ and $J=aM$ are the mass and angular momentum of Kerr black hole, and

$\displaystyle\Sigma$ $\displaystyle=$ $\displaystyle
r^{2}+a^{2}\cos^{2}\theta,$ (2) $\displaystyle\Delta$ $\displaystyle=$
$\displaystyle r^{2}-2Mr+a^{2}\equiv\left(r-r_{+}\right)\left(r-r_{-}\right),$
(3) $\displaystyle A$ $\displaystyle=$
$\displaystyle(r^{2}+a^{2})-a^{2}\Delta\sin^{2}\theta.$ (4)

The event horizon and Cauchy horizon are located at
$r_{+}=M+\sqrt{M^{2}-a^{2}}$ and $r_{-}=M-\sqrt{M^{2}-a^{2}}$, respectively.
The Klein-Gordon field equation $(\nabla^{\nu}\nabla_{\nu}-\mu^{2})\psi=0$ in
such a background reads

$\begin{split}&\left[\frac{\left(r^{2}+a^{2}\right)^{2}}{\Delta}-a^{2}\sin^{2}\theta\right]\frac{\partial^{2}\psi}{\partial
t^{2}}+\frac{4Mar}{\Delta}\frac{\partial^{2}\psi}{\partial t\partial\phi}\\\
&+\left[\frac{a^{2}}{\Delta}-\frac{1}{\sin^{2}\theta}\right]\frac{\partial^{2}\psi}{\partial\phi^{2}}-\frac{\partial}{\partial
r}\left(\Delta\frac{\partial\psi}{\partial r}\right)\\\
&-\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial\psi}{\partial\theta}\right)-\text{$\mu^{2}\Sigma^{2}\psi=0$}\end{split}$
(5)

whose mode solution, Shlapentokh-Rothman (2015), takes the form

$\psi(t,r,\theta,\phi)=e^{-i\omega t}e^{im\phi}Z_{lm}(r)S_{lm}(\theta)$ (6)

with $\omega\in\mathbb{C}$ and azimuthal number $m\in\mathbb{Z}$, and then Eq.
(5) can be separated into

$\displaystyle\frac{1}{\sin\theta}\frac{d}{d\theta}\left(\sin\theta\frac{dS_{lm}}{d\theta}\right)+a^{2}\left(\omega^{2}-\mu^{2}\right)\cos^{2}\theta
S_{lm}$
$\displaystyle-\left(\frac{m^{2}}{\sin^{2}\theta}-\lambda\right)S_{lm}=0,$ (7)

and

$\displaystyle\frac{d}{dr}\left(\Delta\frac{dZ_{lm}}{dr}\right)$
$\displaystyle+$
$\displaystyle\frac{\omega^{2}\left(r^{2}+a^{2}\right)^{2}-4aMm\omega
r+a^{2}m^{2}}{\Delta}Z_{lm}$ (8) $\displaystyle-$
$\displaystyle\left(\mu^{2}r^{2}+a^{2}\omega^{2}+\lambda\right)Z_{lm}=0,$

where $\lambda$ is the separation constant. The radial equation, Eq. (8), is a
confluent Heun equation, and we could rewrite it into the form for convenience
Baumann  _et al._ (2019b):

$\displaystyle\begin{aligned}
0=&\frac{1}{Z_{lm}\Delta}\frac{\mathrm{d}}{\mathrm{d}r}\left(\Delta\frac{\mathrm{d}Z_{lm}}{\mathrm{~{}d}r}\right)-\frac{\lambda}{\Delta}-\left(\mu^{2}-\omega^{2}\right)+\frac{P_{+}^{2}}{\left(r-r_{+}\right)^{2}}\\\
&+\frac{P_{-}^{2}}{\left(r-r_{-}\right)^{2}}-\frac{A_{+}}{\left(r_{+}-r_{-}\right)\left(r-r_{+}\right)}+\frac{A_{-}}{\left(r_{+}-r_{-}\right)\left(r-r_{-}\right)},\end{aligned}$
(9)

where we have defined the following coefficients

$\displaystyle A_{\pm}\equiv
P_{+}^{2}+P_{-}^{2}+\gamma^{2}+\gamma_{\pm}^{2},~{}P_{\pm}\equiv\frac{am-2M\omega
r_{\pm}}{r_{+}-r_{-}}$ (10)

with

$\displaystyle\begin{aligned}
\gamma^{2}&\equiv\frac{1}{4}\left(r_{+}-r_{-}\right)^{2}\left(\mu^{2}-\omega^{2}\right)\\\
\gamma_{\pm}^{2}&\equiv\left[M^{2}\left(\mu^{2}-7\omega^{2}\right)\pm
M\left(r_{+}-r_{-}\right)\left(\mu^{2}-2\omega^{2}\right)\right].\end{aligned}$
(11)

The eigenfunctions $S(\theta)$ in Eq. (7) are the spheroidal harmonics with
its eigenvalues Simone and Will (1992)

$\lambda\equiv\lambda^{lm}={\displaystyle\sum_{j}}f_{2j}^{lm}(ad)^{2j},\
\\{j\in\mathbb{\mathbb{N}}\\}$ (12)

where $(l,m)\in\mathbb{Z}$, $|m|\leq l$, $d^{2}=\omega^{2}-\mu^{2}$ , and
$\lambda\simeq f_{0}^{lm}=l(l+1)$ at the leading order.

In this letter we mainly focus on the ultralight bosons with $\omega\sim\mu\ll
1/M$. Note that $\lambda$ itself is a complex number, and then even if we
choose $\lambda=l(l+1)$, $l$ should be seen as a complex number close to a
integer. We use the matched asymptotic expansion to get the spectrum of
bosonic field analytically Starobinsky (1973); Baumann  _et al._ (2019b).
Radial Teukolsky equation could be solved in “far zone” and “near zone” on the
Kerr background. In the “far zone”, which is dominated by the non-relativistic
effect, the Klein-Gorden equation could be approximated to the Schrodinger
equation with a Newtonian potential, so that the radial equation will have the
same form of the radial hydrogen atom wave function:

$\frac{d^{2}}{dr^{2}}\left(rZ_{lm}\right)+\left[\omega^{2}-\mu^{2}+\frac{2M\mu^{2}}{r}-\frac{l(l+1)}{r^{2}}\right](rZ_{lm})=0.$
(13)

Obviously, the far zone should be defined on $r\gg M$. For the pure outgoing
wave condition at $kr\rightarrow+\infty$, the solution of the above equation
reads

$rZ_{lm}(r)=(2kr)^{l+1}e^{-kr}U(l+1-M\mu^{2}/k,2l+2,2kr)$ (14)

where $U$ is the confluent hypergeometric function, and
$k=\sqrt{\mu^{2}-\omega^{2}}$. Eq. (14) is just the bound state of hydrogen
atoms if $l+1-M\mu^{2}/k\equiv l+1-\nu=-n$ ($n\in\text{$\mathbb{N}$}$). Here
the state with superradiance instability could be approximately described by a
small shift of energy level of hydrogen atoms, namely

$\nu-l-1=n+\delta\nu,\ n\in\mathbb{N},\ \delta\nu\in\mathbb{C},\
\left|\delta\nu\right|\ll n.$ (15)

If we take an asymptotic expansion of Eq. (13) at $2kr\to 0$, using the
property of conflunt hypergeometric functions and gamma functions:

$\displaystyle
U(x,y,\alpha)=\frac{\Gamma(1-y)}{\Gamma(x-y+1)}\sum_{k=0}^{\infty}\frac{(x)_{k}\alpha^{k}}{(y)_{k}k!}$
(16)
$\displaystyle+\frac{\Gamma(y-1)}{\Gamma(x)}\alpha^{1-y}\sum_{k=0}^{\infty}\frac{(x-y+1)_{k}\alpha^{k}}{(2-y)_{k}k!};\
y\notin\bm{Z}$

$\Gamma(\alpha)\propto\frac{(-1)^{n}(1+O(\alpha+n))}{n!(\alpha+n)};\
(\alpha\rightarrow-n)\wedge n\in\mathbb{N}$ (17)

where $(x)_{0}\equiv 1,\
(x)_{k}\equiv\text{$\Gamma\left(a+k\right)/\Gamma\left(a\right)$},$ the
asymptotic expansion of $Z_{lm}(r)$ in Eq. (14) at $M\ll r\ll 1/k$ takes the
form of

$\displaystyle Z_{lm}(r)$
$\displaystyle\sim(-1)^{n}\frac{(2l+1+n)!}{(2l+1)!}(2kr)^{l}+\cdots$ (18)
$\displaystyle+(-1)^{n+1}\delta\nu(2l)!n!(2kr)^{-l-1}+\cdots\ .$

On the other hand, in the “near zone”, which should close to the event horizon
at $r_{+}$, the main term of Radial Teukolsky equation is $-l(l+1)/\Delta$ and
$P_{+}^{2}/\left(r-r_{+}\right)^{2}$, which could be easily read by Eq. (9).
Introduce the rescale radial coordinate:

$\displaystyle z$
$\displaystyle=\left(r-r_{+}\right)/\left(r_{+}-r_{-}\right),$ (19)

then Eq. (8) can be written as

$\displaystyle
z\left(z+1\right)\frac{d}{dz}\left[z(z+1)\frac{dZ_{lm}}{dz}\right]$
$\displaystyle+\left[P_{+}^{2}-l\left(l+1\right)z\left(z+1\right)\right]Z_{lm}=0.$
(20)

Meanwhile, once fixed the radial coordinate, expanding Eq. (9) to the leading
order of $\mu M$ could lead to:

$\displaystyle\begin{aligned}
&\left[\frac{\mathrm{d}^{2}}{\mathrm{~{}d}z^{2}}+\left(\frac{1}{z}+\frac{1}{z+1}\right)\frac{\mathrm{d}}{\mathrm{d}z}-\frac{\ell(\ell+1)}{z(z+1)}\right]Z_{lm}\\\
&+\left[\frac{P_{+}^{2}}{z^{2}}+\frac{P_{+}^{2}}{(z+1)^{2}}-\frac{2P_{+}^{2}}{z}+\frac{2P_{+}^{2}}{z+1}\right]Z_{lm}=0\end{aligned}$
(21)

If we compare Eq. (20) and Eq. (21), we could estimate that Eq. (20) should be
valid on $r_{+}<r\ll\max(l/\omega,l/m)$. So these two regions have the overlap
on $M\ll r\ll\max(l/\omega,l/m)$. In order to match the solution in Eq. (18)
in the overlapping region $M\ll r\ll\max(l/\omega,l/m)$, the solution of Eq.
(20) is given by

$Z_{lm}(z)=\left(\frac{z}{z+1}\right)^{iP_{+}}\left[Aw_{5}\left(z\right)+Bw_{6}\left(z\right)\right]$
(22)

where

$\displaystyle
w_{5}(z)=\left[-\left(z+1\right)\right]^{l}{}_{2}F_{1}(-l,-l+2iP_{+};-2l;\left(z+1\right)^{-1})$
(23)
$\displaystyle\sim\left(-z\right)^{l}\sim\left(\frac{-r}{r_{+}-r_{-}}\right)^{l}$

and

$\displaystyle w_{6}\left(z\right)$
$\displaystyle=\left[-\left(z+1\right)\right]^{-l-1}$ (24)
$\displaystyle\times{}_{2}F_{1}(l+1,l+1+2iP_{+};2l+2;\left(z+1\right)^{-1})$
$\displaystyle\sim\left(-z\right)^{-l-1}\sim\left(\frac{-r}{r_{+}-r_{-}}\right)^{-l-1},$

are two linear independent solution of hypergeometric equation. So that we
find the coefficients $A,B$ are:

$A=\left(-1\right)^{n+l}\left[2k(r_{+}-r_{-})\right]^{l}\frac{\left(2l+n+1\right)!}{\left(2l+1\right)!}$
(25)

$B=\left(-1\right)^{n+l}\left[2k(r_{+}-r_{-})\right]^{-l-1}\delta\text{$\nu$}\left(2l\right)!\left(n\right)!.$
(26)

This is the matching condition gluning two asymptotic regions.

In order to figure out the solution of Eq. (8), we still need to take the
boundary condition near horizon into account. At classical level, we should
only have an ingoing wave at the horizon because the effective potential drops
to zero. Here we suppose that the microstructure due to the quantum effects
modify the classical description of a black hole at horizon scales and such
effects can be encoded in a frequency dependent reflectivity
$\mathcal{R}(\omega)$:

$Z_{lm}(r_{*})\sim
e^{-i\left(\omega-m\Omega_{H}\right)r_{*}}+\mathcal{R}(\omega)e^{i\left(\omega-m\Omega_{H}\right)r_{*}},$
(27)

where

$\displaystyle\Omega_{H}=a/(2Mr_{+}),$ (28)

for $r_{*}\rightarrow-\infty$, where $r_{*}$ is the tortoise coordinate

$r_{*}=r+{2Mr_{+}\over r_{+}-r_{-}}\ln{r-r_{+}\over 2M}-{2Mr_{-}\over
r_{+}-r_{-}}\ln{r-r_{-}\over 2M}.$ (29)

Actually $w_{5}$ and $w_{6}$ can be re-written by another set of function
bases of hypergeometric equation $w_{3}$ and $w_{4}$ for respectively
referring to ingoing and outgoing solutions near horizon, namely

$\displaystyle w_{5}$ $\displaystyle=$ $\displaystyle
T_{53}w_{3}+T_{54}w_{4},$ (30) $\displaystyle w_{6}$ $\displaystyle=$
$\displaystyle T_{63}w_{3}+T_{64}w_{4},$ (31)

where

$\displaystyle w_{3}(z)$ $\displaystyle=$
$\displaystyle{}_{2}F_{1}(-l,l+1;2iP_{+}+1;-z),$ (32) $\displaystyle w_{4}(z)$
$\displaystyle=$
$\displaystyle{}_{2}F_{1}(l+1-2iP_{+},-l-2iP_{+};-2iP_{+}+1;-z),$

and

$\displaystyle T_{53}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{l}\Gamma(1-2iP_{+})\Gamma(2iP_{+})\Gamma(-2l)}{\text{$\Gamma(-l)\Gamma(-l-2iP_{+})\Gamma(1+2iP_{+})$}},$
(34) $\displaystyle T_{54}$ $\displaystyle=$
$\displaystyle-\frac{(-1)^{l}\Gamma(2iP_{+})\Gamma(-2l)}{\Gamma(-l+2iP_{+})\Gamma(-l)},$
(35) $\displaystyle T_{63}$ $\displaystyle=$
$\displaystyle-\frac{(-1)^{l}\Gamma(1-2iP_{+})\Gamma(2iP_{+})\Gamma(2l+2)}{\text{$\Gamma(l+1)\Gamma(l+1-2iP_{+})\Gamma(1+2iP_{+})$}},$
(36) $\displaystyle T_{64}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{l}\Gamma(2iP_{+})\Gamma(2l+2)}{\Gamma(l+1+2iP_{+})\Gamma(l+1)}.$
(37)

And then the boundary condition in Eq. (27) gives

$\mathcal{R}=\frac{AT_{54}+BT_{64}}{AT_{53}+BT_{63}}.$ (38)

Supposing $\mathcal{R}=\mathcal{R}_{0}e^{i\phi_{\mathcal{R}}}$, after a
tedious but straightforward computation, we obtain

$\delta\nu=\delta\bar{\nu}\frac{1-\mathcal{R}_{0}e^{i\left(\phi_{\mathcal{R}}+2\sum_{j=1}^{l}\phi_{j}\right)}}{1+\mathcal{R}_{0}e^{i\left(\phi_{\mathcal{R}}+2\sum_{j=1}^{l}\phi_{j}\right)}},$
(39)

where $\phi_{j}=\arctan(2P_{+}/j)$, and

$\displaystyle\delta\bar{\nu}=iP_{+}\left[2k(r_{+}-r_{-})\right]^{2l+1}\frac{(2l+n+1)!}{n!}$
(40)
$\displaystyle\left[\frac{l!}{(2l)!(2l+1)!}\right]^{2}\prod_{j=1}^{l}(j^{2}+4P_{+}^{2})$

denotes the result without reflection at horizon Pani  _et al._ (2012b).

Considering $M\mu^{2}/\sqrt{\mu^{2}-\omega^{2}}=\nu$ and
$\nu=\nu_{0}+\delta\nu$ with $\nu_{0}=n+l+1$, we get

$\omega=\omega_{0}+\delta\omega,$ (41)

where

$\displaystyle\omega_{0}^{2}$ $\displaystyle=$
$\displaystyle\mu^{2}-\mu^{2}\left({M\mu\over l+n+1}\right)^{2},$ (42)
$\displaystyle\delta\omega$ $\displaystyle\simeq$
$\displaystyle\frac{\delta\nu}{M}\left(\frac{M\mu}{l+n+1}\right)^{3}.$ (43)

Since the real and imaginary parts of $\delta\omega$ denote the level shift of
massive scalar field and the stability/instability property of the
corresponding mode, we decompose $\delta\omega$ on the Cartesian complex plane
as follows

$\displaystyle\text{Re}[\delta\omega]$ $\displaystyle=$
$\displaystyle\frac{\delta\bar{\nu}}{iM}\left(\frac{M\mu}{l+n+1}\right)^{3}\frac{2\mathcal{R}_{0}\sin\phi_{\omega}}{\mathcal{R}_{0}^{2}+1+2\mathcal{R}_{0}\cos\phi_{\omega}},$
(44) $\displaystyle\text{Im}[\delta\omega]$ $\displaystyle=$
$\displaystyle\frac{\delta\bar{\nu}}{iM}\left(\frac{M\mu}{l+n+1}\right)^{3}\frac{1-\mathcal{R}_{0}^{2}}{\mathcal{R}_{0}^{2}+1+2\mathcal{R}_{0}\cos\phi_{\omega}},$
(45)

where $\phi_{\omega}=\phi_{\mathcal{R}}+2\text{$\sum_{j=1}^{l}\phi_{j}$}$.

Discussion. In this letter, we derive analytical formulas of superradiance
instability influenced by the modified boundary conditions near the event
horizon. Our result shows that the behavior of superradiance instability of
ultralight scalar bosons could be affected by the potential near-horizon new
physics. Since the evolution of a boson-BH system is related to the bound
state frequencies, Eq. (44), our results imply that the properties such as the
e-folding time a boson cloud takes to reach the saturation point, the maximum
boson cloud mass, the gravitational wave signal would be influenced by the
near-horizon microstructure. Further astrophysical observation could be used
to detect both the existence of ultralight bosons and the possible breakdown
of classical BHs.

Acknowledgments.  This work is supported by the National Key Research and
Development Program of China Grant No.2020YFC2201502, grants from NSFC (grant
No. 11975019, 11690021, 11991052, 12047503), the Key Research Program of the
Chinese Academy of Sciences (Grant NO. XDPB15), Key Research Program of
Frontier Sciences, CAS, Grant NO. ZDBS-LY-7009, and the science research
grants from the China Manned Space Project with NO. CMS-CSST-2021-B01.

## References

  * Arvanitaki  _et al._ (2010a) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D81, 123530 (2010a), arXiv:0905.4720 [hep-th] .
  * Essig  _et al._ (2013) R. Essig _et al._ , in  _Community Summer Study 2013: Snowmass on the Mississippi_ (2013) arXiv:1311.0029 [hep-ph] .
  * Irastorza and Redondo (2018) I. G. Irastorza and J. Redondo, Prog.Part.Nucl.Phys. 102, 89 (2018), arXiv:1801.08127 [hep-ph] .
  * Goodsell  _et al._ (2009) M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, JHEP 11, 027 (2009), arXiv:0909.0515 [hep-ph] .
  * Jaeckel and Ringwald (2010) J. Jaeckel and A. Ringwald, Ann.Rev.Nucl.Part.Sci. 60, 405 (2010), arXiv:1002.0329 [hep-ph] .
  * Graham  _et al._ (2016) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016), arXiv:1504.02102 [hep-ph] .
  * Agrawal  _et al._ (2020) P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Phys. Lett. B 801, 135136 (2020), arXiv:1810.07188 [hep-ph] .
  * Preskill  _et al._ (1983) J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie (1983) L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler (1983) M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Svrcek and Witten (2006) P. Svrcek and E. Witten, JHEP 06, 051 (2006), arXiv:hep-th/0605206 [hep-th] .
  * Arvanitaki and Dubovsky (2011) A. Arvanitaki and S. Dubovsky, Phys. Rev. D83, 044026 (2011), arXiv:1004.3558 [hep-th] .
  * Brito  _et al._ (2015a) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015a), arXiv:1501.06570 [gr-qc] .
  * Marsh (2016) D. J. E. Marsh, Phys. Rept. 643, 1 (2016), arXiv:1510.07633 [astro-ph.CO] .
  * Hui  _et al._ (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Annulli  _et al._ (2020) L. Annulli, V. Cardoso, and R. Vicente, Phys. Rev. D 102, 063022 (2020), arXiv:2009.00012 [gr-qc] .
  * Chadha-Day  _et al._ (2021) F. Chadha-Day, J. Ellis, and D. J. E. Marsh, (2021), arXiv:2105.01406 [hep-ph] .
  * Damour  _et al._ (1976) T. Damour, N. Deruelle, and R. Ruffini, Lett.Nuovo Cim. 15, 257 (1976).
  * Zouros and Eardley (1979) T. Zouros and D. Eardley, Annals Phys. 118, 139 (1979).
  * Detweiler (1980) S. L. Detweiler, Phys.Rev. D22, 2323 (1980).
  * Dolan (2007) S. R. Dolan, Phys. Rev. D 76 (2007), 10.1103/PhysRevD.76.084001.
  * Arvanitaki  _et al._ (2010b) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, Phys. Rev. D 81 (2010b), 10.1103/PhysRevD.81.123530.
  * Shlapentokh-Rothman (2014) Y. Shlapentokh-Rothman, Commun.Math.Phys. 329, 859 (2014), arXiv:1302.3448 [gr-qc] .
  * Pani  _et al._ (2012a) P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. Lett. 109, 131102 (2012a).
  * Pani  _et al._ (2012b) P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Phys. Rev. D 86, 104017 (2012b).
  * Witek  _et al._ (2013) H. Witek, V. Cardoso, A. Ishibashi, and U. Sperhake, Phys.Rev. D87, 043513 (2013), arXiv:1212.0551 [gr-qc] .
  * Brito  _et al._ (2013) R. Brito, V. Cardoso, and P. Pani, Phys. Rev. D 88 (2013), 10.1103/PhysRevD.88.023514.
  * Endlich and Penco (2017) S. Endlich and R. Penco, JHEP 05, 052 (2017), arXiv:1609.06723 [hep-th] .
  * East (2017) W. E. East, Phys. Rev. D96, 024004 (2017), arXiv:1705.01544 [gr-qc] .
  * East and Pretorius (2017) W. E. East and F. Pretorius, Phys. Rev. Lett. 119, 041101 (2017), arXiv:1704.04791 [gr-qc] .
  * Baryakhtar  _et al._ (2017a) M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D 96, 035019 (2017a), arXiv:1704.05081 [hep-ph] .
  * Cardoso  _et al._ (2018) V. Cardoso, Ó. J. Dias, G. S. Hartnett, M. Middleton, P. Pani, and J. E. Santos, Journal of Cosmology and Astroparticle Physics 2018, 043 (2018).
  * East (2018) W. E. East, Phys. Rev. Lett. 121, 131104 (2018), arXiv:1807.00043 [gr-qc] .
  * Frolov  _et al._ (2018) V. P. Frolov, P. Krtouš, D. Kubizňák, and J. E. Santos, Phys. Rev. Lett. 120, 231103 (2018), arXiv:1804.00030 [hep-th] .
  * Dolan (2018) S. R. Dolan, Phys. Rev. D98, 104006 (2018), arXiv:1806.01604 [gr-qc] .
  * Baumann  _et al._ (2019a) D. Baumann, H. S. Chia, J. Stout, and L. ter Haar, JCAP 12, 006 (2019a), arXiv:1908.10370 [gr-qc] .
  * Brito  _et al._ (2020) R. Brito, S. Grillo, and P. Pani, Phys. Rev. Lett. 124, 211101 (2020), arXiv:2002.04055 [gr-qc] .
  * Tsukada  _et al._ (2021) L. Tsukada, R. Brito, W. E. East, and N. Siemonsen, Phys. Rev. D 103, 083005 (2021), arXiv:2011.06995 [astro-ph.HE] .
  * Yuan  _et al._ (2021) C. Yuan, R. Brito, and V. Cardoso, Phys. Rev. D 104, 044011 (2021), arXiv:2106.00021 [gr-qc] .
  * Arvanitaki  _et al._ (2015) A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D 91, 084011 (2015), arXiv:1411.2263 [hep-ph] .
  * Arvanitaki  _et al._ (2017) A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky, and R. Lasenby, Phys. Rev. D 95, 043001 (2017), arXiv:1604.03958 [hep-ph] .
  * Baryakhtar  _et al._ (2017b) M. Baryakhtar, R. Lasenby, and M. Teo, Phys. Rev. D 96, 035019 (2017b), arXiv:1704.05081 [hep-ph] .
  * Brito  _et al._ (2017a) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. Lett. 119, 131101 (2017a), arXiv:1706.05097 [gr-qc] .
  * Brito  _et al._ (2017b) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. D 96, 064050 (2017b), arXiv:1706.06311 [gr-qc] .
  * Isi  _et al._ (2019) M. Isi, L. Sun, R. Brito, and A. Melatos, Phys. Rev. D 99, 084042 (2019), [Erratum: Phys.Rev.D 102, 049901 (2020)], arXiv:1810.03812 [gr-qc] .
  * Ghosh  _et al._ (2019) S. Ghosh, E. Berti, R. Brito, and M. Richartz, Phys. Rev. D 99, 104030 (2019), arXiv:1812.01620 [gr-qc] .
  * Palomba  _et al._ (2019) C. Palomba _et al._ , Phys. Rev. Lett. 123, 171101 (2019), arXiv:1909.08854 [astro-ph.HE] .
  * Sun  _et al._ (2020) L. Sun, R. Brito, and M. Isi, Phys. Rev. D 101, 063020 (2020), [Erratum: Phys.Rev.D 102, 089902 (2020)], arXiv:1909.11267 [gr-qc] .
  * Zhu  _et al._ (2020) S. J. Zhu, M. Baryakhtar, M. A. Papa, D. Tsuna, N. Kawanaka, and H.-B. Eggenstein, Phys. Rev. D 102, 063020 (2020), arXiv:2003.03359 [gr-qc] .
  * Ng  _et al._ (2020) K. K. Y. Ng, M. Isi, C.-J. Haster, and S. Vitale, Phys. Rev. D 102, 083020 (2020), arXiv:2007.12793 [gr-qc] .
  * Kodama and Yoshino (2012) H. Kodama and H. Yoshino,  _Proceedings, 2011 Asia Pacific School/Workshop on Cosmology and Gravitation: Shanghai, China, February 10-14, 2011_ , Int. J. Mod. Phys. Conf. Ser. 7, 84 (2012), arXiv:1108.1365 [hep-th] .
  * Brito  _et al._ (2015b) R. Brito, V. Cardoso, and P. Pani, Class. Quant. Grav. 32, 134001 (2015b), arXiv:1411.0686 [gr-qc] .
  * Tsukada  _et al._ (2019) L. Tsukada, T. Callister, A. Matas, and P. Meyers, Phys. Rev. D 99, 103015 (2019), arXiv:1812.09622 [astro-ph.HE] .
  * Thorne  _et al._ (1986) K. S. Thorne, R. H. Price, and D. A. Macdonald, eds.,  _BLACK HOLES: THE MEMBRANE PARADIGM_ (1986).
  * Hawking and Ellis (2011) S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2011).
  * Barcelo  _et al._ (2008) C. Barcelo, S. Liberati, S. Sonego, and M. Visser, Phys. Rev. D77, 044032 (2008), arXiv:0712.1130 [gr-qc] .
  * Morris  _et al._ (1988) M. S. Morris, K. S. Thorne, and U. Yurtsever, Phys. Rev. Lett. 61, 1446 (1988).
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Holdom and Ren (2017) B. Holdom and J. Ren, Phys. Rev. D95, 084034 (2017), arXiv:1612.04889 [gr-qc] .
  * Bekenstein (1974) J. D. Bekenstein, Lett. Nuovo Cim. 11, 467 (1974).
  * Bekenstein and Mukhanov (1995) J. D. Bekenstein and V. F. Mukhanov, Phys. Lett. B360, 7 (1995), arXiv:gr-qc/9505012 [gr-qc] .
  * Skenderis and Taylor (2008) K. Skenderis and M. Taylor, Phys. Rept. 467, 117 (2008), arXiv:0804.0552 [hep-th] .
  * Polchinski (2017) J. Polchinski, in  _Proceedings, Theoretical Advanced Study Institute in Elementary Particle Physics: New Frontiers in Fields and Strings (TASI 2015): Boulder, CO, USA, June 1-26, 2015_ (2017) pp. 353–397, arXiv:1609.04036 [hep-th] .
  * Almheiri  _et al._ (2013) A. Almheiri, D. Marolf, J. Polchinski, and J. Sully, JHEP 02, 062 (2013), arXiv:1207.3123 [hep-th] .
  * Cardoso  _et al._ (2016) V. Cardoso, E. Franzin, and P. Pani, Phys. Rev. Lett. 116, 171101 (2016), [Erratum: Phys. Rev. Lett.117,no.8,089902(2016)], arXiv:1602.07309 [gr-qc] .
  * Cardoso  _et al._ (2019) V. Cardoso, V. F. Foit, and M. Kleban, JCAP 1908, 006 (2019), arXiv:1902.10164 [hep-th] .
  * Mark  _et al._ (2017) Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, Phys. Rev. D96, 084002 (2017), arXiv:1706.06155 [gr-qc] .
  * Price and Khanna (2017) R. H. Price and G. Khanna, Class. Quant. Grav. 34, 225005 (2017), arXiv:1702.04833 [gr-qc] .
  * Fang  _et al._ (2021) Y. Fang, R.-Z. Guo, and Q.-G. Huang, (2021), arXiv:2108.04511 [astro-ph.CO] .
  * Oshita and Afshordi (2019) N. Oshita and N. Afshordi, Phys. Rev. D99, 044002 (2019), arXiv:1807.10287 [gr-qc] .
  * Burgess  _et al._ (2018) C. P. Burgess, R. Plestid, and M. Rummel, JHEP 09, 113 (2018), arXiv:1808.00847 [gr-qc] .
  * Rummel and Burgess (2019) M. Rummel and C. P. Burgess, (2019), arXiv:2001.00041 [gr-qc] .
  * Brito  _et al._ (2015c) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys 906, 1501 (2015c).
  * Shlapentokh-Rothman (2015) Y. Shlapentokh-Rothman (2015) pp. 289–345, arXiv:1302.6902 [gr-qc] .
  * Baumann  _et al._ (2019b) D. Baumann, H. S. Chia, J. Stout, and L. T. Haar, Journal of Cosmology and Astroparticle Physics 2019 (2019b), 10.1088/1475-7516/2019/12/006, arXiv:1908.10370 .
  * Simone and Will (1992) L. E. Simone and C. M. Will, Class. Quant. Grav. 9, 963 (1992).
  * Starobinsky (1973) A. Starobinsky, Sov. Phys. JETP 37, 28 (1973).


figure table

# Pharmacokinetic parameters quantification in DCE-MRI for prostate cancer

Jhonalbert Aponte  Servicio de Radioterapia Avanzada, Fundación Arturo López
Pérez, Santiago de Chile, Chile  Álvaro Ruiz  Servicio de Radioterapia
Avanzada, Fundación Arturo López Pérez, Santiago de Chile, Chile  Jacksson
Sánchez  Faculty of Science and Technology, Physics Department, Universidad
Nacional Pedro Henríquez Ureña, Santo Domingo, Dominican Republic  Miguel
Martín-Landrove  Centre for Medical Visualization, National Institute for
Bioengineering, Universidad Central de Venezuela, Caracas, Venezuela  Centro
de Diagnóstico Docente Las Mercedes, Caracas, Venezuela

###### Abstract

Tumor vascularity detection and quantification are of high relevance in the
assessment of cancer lesions not only for disease diagnostics but for therapy
considerations and monitoring. The present work addressed the quantification
of pharmacokinetic parameters derived from the two-compartment Brix model by
analyzing and processing Dynamic Contrast-Enhanced Magnetic Resonance Images
(DCE-MRI) of prostate cancer lesions. The 3D image sets were acquired at
regular time intervals, covering all the phases implied in contrast injection
(wash-in and wash-out phases), and the standardized image intensity is
determined for each voxel, conforming to a 4D data set. Previous voxel
classification was carried out by the three-time-point method proposed by
Degani et al. (1997) and Furman-Haran et al. (1998) to identify regions of
interest. Relevant pharmacokinetic parameters, such as $k_{el}$, the vascular
elimination rate, and $k_{ep}$, the extravascular transfer rate, are extracted
by a novel interpolation method applicable to compartment models. Parameter
distribution maps were obtained for either pathological or unaffected
glandular regions indicating that a three-compartment model, including fast
and slow exchange compartments, provides a more suitable description of the
contrast kinetics. Results can be applied to prostate cancer diagnostic
evaluation and therapy follow-up.

###### keywords:

DCE-MRI, Levenberg-Marquardt, prostate cancer, tumor vascularity, two-
compartment pharmacokinetic model

*Miguel Martín-Landrove<EMAIL_ADDRESS>
##  1 Introduction

Dynamic Contrast-Enhanced MRI or DCE-MRI has been applied extensively to
diagnose and quantify several pathologies associated with cancer [1, 2, 3, 4,
5, 6]. It essentially consists of the controlled intravenous delivery of a
contrast agent, typically a Gadolinium compound, followed by the acquisition
of a volumetric $T_{1}$-weighted MRI. The resulting image dataset can be
analyzed either qualitatively [7, 8, 9, 10, 11, 12] or quantitatively [4, 5,
6, 13, 14]. Several compartmental models are proposed for the pharmacokinetics
of the contrast agent [14, 15, 16, 17, 18]. In the case of prostate DCE-MRI,
pharmacokinetic two-compartment models are applicable [5, 15, 16]. The present
work is addressed to the quantification of pharmacokinetic parameters, as
derived from the two-compartment Brix-Tofts model [15, 16] and multi-
compartmental models [19, 20, 21, 22], through the analysis and processing of
DCE-MRI of prostate cancer lesions, using the Levenberg-Marquardt algorithm
[23, 24] and modified de Prony method [25, 26].

The work is organized as follows, in Section 2, general aspects of the two-
compartment Brix-Tofts model, three-compartment models, image acquisition and
qualitative analysis of the image data set, and quantitative evaluation of the
pharmacokinetic parameters are discussed; in Section 3, the results are
discussed and some conclusions are presented.

##  2 Materials and methods

###  2.1 Two-compartment Brix-Tofts model

Essentially, the model establishes that the contrast intake occurs in two
connected environments or compartments: one, which is called the central
compartment related to the capillary and vascular space, and the second, the
peripheral compartment, generally related to extravascular and extracellular
space. The model can be depicted schematically as shown in Figure 1, after
[15],

Figure 1: Schematic representation of e two-compartment Brix-Tofts model

The parameters of the model $k_{in}$, $k^{trnas}$, $k_{ep}$, $k_{el}$, $v_{p}$
and $v_{e}$ describe the pharmacokinetics in the system, i.e., concentration
exchange velocities and volumes of the compartments, respectively. The
differential equations that describe the pharmacokinetics are,

$\displaystyle\frac{dM_{1}}{dt}=$ $\displaystyle
k_{in}-(k^{trans}+k_{el})M_{1}+k_{ep}M_{2}$ (1)
$\displaystyle\frac{dM_{2}}{dt}=$ $\displaystyle k^{trans}M_{1}-k_{ep}M_{2}$

where $M_{1}$ and $M_{2}$ are the contrast agent total mass in compartments 1
and 2 respectively. Assuming that the transfer velocity between both
compartments is equal,

$k^{trans}v_{p}=k_{ep}v_{e}$ (2)

and $v_{e}\ll v_{p}$, equation 1 can be written,

$\frac{dC_{1}}{dt}=\frac{k_{in}}{v_{p}}-k_{el}C_{1},\\\
\frac{dC_{2}}{dt}=\frac{v_{p}}{v_{e}}k^{trans}C_{1}-k_{ep}C_{2}$ (3)

where $C_{1}$ and $C_{2}$ are the contrast concentrations in the compartments.
Because the contrast agent is delivered as a bolus, the set defined in
equation 3 has a solution,

$\displaystyle C_{1}(t)=$
$\displaystyle\frac{k_{in}}{v_{p}k_{el}}(e^{k_{el}t^{\prime}}-1)e^{-k_{el}t}$
(4) $\displaystyle C_{2}(t)=$
$\displaystyle\frac{k_{in}k^{trans}}{v_{e}}\times$
$\displaystyle\times\left[v(e^{k_{el}t^{\prime}}-1)e^{-k_{el}t}-u(e^{k_{ep}t^{\prime}}-1)e^{-k_{ep}t}\right]$

where coefficients $u$ and $v$ are given by,

$\displaystyle u=$ $\displaystyle\left[k_{ep}(k_{ep}-k_{el})\right]^{-1}$ (5)
$\displaystyle v=$ $\displaystyle\left[k_{el}(k_{ep}-k_{el})\right]^{-1}$

In the solution set, equation 4, the time parameter $t^{\prime}$ defines the
different phases for the progression of the contrast agent. During the wash-in
phase, $0\leq t\leq\tau$, where $\tau$ is the contrast infusion time,
$t^{\prime}$ is taken equal to $t$. In the wash-out phase, $t>\tau$,
$t^{\prime}$ is set equal to $\tau$.

###  2.2 Fast and Slow Exchange Compartments

The Brix-Tofts model represents a simple approach to a real system. Several
authors [19, 20, 21, 22] have proposed multi-compartmental models to describe
the kinetics of the contrast in the tissue. Among these models, the simplest
one corresponds to the assumption that the system can be described by a three-
compartment model [20, 22] with the following compartmental equations,

$\displaystyle\frac{dC_{p}}{dt}=$
$\displaystyle\frac{k_{in}}{v_{p}}-(k_{s}^{trans}+k_{f}^{trans}+k_{el})C_{p}+$
(6)
$\displaystyle+\frac{k_{eps}v_{es}}{v_{p}}C_{s}+\frac{k_{epf}v_{ef}}{v_{p}}C_{f}$
$\displaystyle\frac{dC_{s}}{dt}=$
$\displaystyle\frac{k_{s}^{trans}v_{p}}{v_{es}}C_{p}-k_{eps}C_{s}$
$\displaystyle\frac{dC_{f}}{dt}=$
$\displaystyle\frac{k_{f}^{trans}v_{p}}{v_{ef}}C_{p}-k_{epf}C_{f}$

Under the same assumptions made for the Brix-Tofts model, equations (6) have
the solution,

$\displaystyle
C_{p}(t)=\frac{k_{in}}{v_{p}k_{el}}(e^{k_{el}t^{\prime}}-1)e^{-k_{el}t}$
$\displaystyle C_{s}(t)=\frac{k_{in}k_{s}^{trans}}{v_{es}}\times$
$\displaystyle\qquad\quad\times\left[v_{s}(e^{k_{el}t^{\prime}}-1)e^{-k_{el}t}-u_{s}(e^{k_{eps}t^{\prime}}-1)e^{-k_{eps}t}\right]$
$\displaystyle C_{f}(t)=\frac{k_{in}k_{f}^{trans}}{v_{ef}}\times$
$\displaystyle\qquad\quad\times\left[v_{f}(e^{k_{el}t^{\prime}}-1)e^{-k_{el}t}-u_{f}(e^{k_{epf}t^{\prime}}-1)e^{-k_{epf}t}\right]$
(7)

with,

$\displaystyle u_{s,f}=$
$\displaystyle\left[k_{ep;s,f}(k_{ep;s,f}-k_{el})\right]^{-1}$ (8)
$\displaystyle v_{s,f}=$
$\displaystyle\left[k_{el}(k_{ep;s,f}-k_{el})\right]^{-1}$

and similarly, in equations (2.2) during the wash-in phase, $0\leq t\leq\tau$,
where $\tau$ is the contrast infusion time, $t^{\prime}$ is taken equal to
$t$, and in the wash-out phase, $t>\tau$, $t^{\prime}$ is set equal to $\tau$.
The contrast concentration in the tissue is the sum of $C_{s}$ and $C_{f}$.

###  2.3 Image acquisition and quantitative analysis

$T_{1}$-weighted MRI intensity is related to the extravascular and
extracellular contrast concentration $C_{2}$ through the following relation
[5],

$S_{t}\approx S_{0}(1+FC_{2}(t))$ (9)

which holds if certain conditions apply for the MRI sequence parameters, i.e.,
$TR\alpha C_{2}\ll 1$ and $TE\beta C_{2}$, where $\alpha$ and $\beta$
determine the enhancement of longitudinal and transversal relaxation rates,
respectively, due to the contrast agent, and $S_{0}$ is the image intensity
without contrast. In general, DCE - MRI protocols are fine-tuned to fulfill
the necessary conditions for equation 9 to apply, and in such a case,
equations for the concentration $C_{2}$ can be obtained from the image data
set as,

$\displaystyle
C_{2}=A\left[v(1-e^{-k_{el}t})-u(1-e^{-k_{ep}t})\right],t\leq\tau$
$\displaystyle
C_{2}=A\left[v(e^{k_{el}\tau}-1)e^{-k_{el}\tau}-u(e^{k_{ep}\tau}-1)e^{-k_{ep}t}\right],t>\tau$
(10)

Taking into account that $u$ and $v$ are given by equation 5, only four
parameters are needed to fit the model to experimental data: $k_{e}l$,
$k_{e}p$, $\tau$ and $A$. Prostate DCE-MRI was obtained from the Collection
Prostate-Diagnosis at The Cancer Imaging Archive (TCIA), National Cancer
Institute [27].

The image data set was previously registered to the image for the initial
condition, i.e., $T_{1}$-weighted MRI with no contrast, to diminish
physiological and involuntary patient movements. 3D rigid image registration
was performed using MATLAB Image Processing Toolbox, with Mutual Information
[28] as a similarity metric for the Regular Step Gradient Descent [29, 30]
optimization algorithm.

####  2.3.1 Levenberg-Marquardt fit.

Quantification of pharmacokinetic parameters was also performed with MATLAB,
using least square fitting by the Levenberg-Marquardt algorithm [23, 24]. As
previously mentioned, only four parameters are needed to fit the model to
experimental data; if it is assumed that only the ratio between $u$ and $v$ is
used instead of their actual expressions, equation 5 and the proportionality
coefficient $A$ is incorporated to either $u$ or $v$, then the number of
parameters is maintained, and equation 2.3 can be rewritten as,

$\displaystyle
C_{2}=A^{\prime}\left[\frac{k_{ep}}{k_{el}}(1-e^{-k_{el}t})-(1-e^{-k_{ep}t})\right],t\leq\tau$
$\displaystyle
C_{2}=A^{\prime}\left[\frac{k_{ep}}{k_{el}}(e^{k_{el}\tau}-1)e^{k_{el}t}-(e^{k_{ep}\tau}-1)e^{-k_{ep}t}\right],t>\tau$
(11)

On the other hand, if the two-exponential dependence of the model is
preserved, i.e., it is a two-compartment model, some freedom can be gained for
the fitting of experimental data and only five parameters are needed,

$\displaystyle C_{2}=A_{1}(1-e^{-k_{el}t})-A_{2}(1-e^{-k_{ep}t}),t\leq\tau$
$\displaystyle
C_{2}=A_{1}(e^{k_{el}\tau}-1)e^{-k_{el}t}-A_{2}(e^{k_{ep}\tau}-1)e^{-k_{ep}t},t>\tau$
(12)

An example of the application of the model with a different number of fitting
parameters is shown in Figure 2. As expected, the model with five parameters,
yields the best fit, while the four-parameter model imposes too many
restrictions between the parameters avoiding data fitting. Nevertheless, the
four-parameter model could be used to confirm the Brix-Tofts model’s exact
validity within the data noise limits.

Figure 2: Different data fitting results for the same pharmacokinetic model,
fitted curve is indicated in red. (a) Four parameters, (b) Five parameters. In
both cases, the Levenberg Marquardt algorithm was used

####  2.3.2 Modified de Prony method fit. Two Compartments Model

A modified version of the de Prony method [25] was implemented for multi-echo
$T_{2}$-weighted MRI [26] to quantify transversal rate distributions in solid
brain tumors. Essentially, the method assumes that the signal or image
intensity is described by a superposition of a finite number of decaying
exponential functions and it is sampled at regularly fixed time intervals as
is the case for DCE-MRI. If the Brix-Tofts model is assumed and following
equation 2.3.1, voxel intensity is given by,

$\displaystyle p_{i}=(A_{1}-A_{2})-A_{1}X_{1}^{i}+A_{2}X_{2}^{i},1\leq i\leq
n,i\delta t\leq\tau$ $\displaystyle
p_{n+j}=B_{1}X_{1}^{n+j}-B_{2}X_{2}^{n+j},1\leq j\leq m,(n+j)\delta t>\tau$
(13)

where the following definitions apply,

$\displaystyle X_{1}$ $\displaystyle=e^{-k_{el}\delta t}$ (14) $\displaystyle
X_{2}$ $\displaystyle=e^{-k_{ep}\delta t}$ $\displaystyle A_{1}$
$\displaystyle=F\frac{k_{in}k^{trans}}{v_{e}}\left[k_{el}(k_{ep}-k_{el})\right]^{-1}$
$\displaystyle A_{2}$
$\displaystyle=F\frac{k_{in}k^{trans}}{v_{e}}\left[k_{ep}(k_{ep}-k_{el})\right]^{-1}$
$\displaystyle B_{1}$ $\displaystyle=(e^{k_{el}\tau}-1)A_{1}$ $\displaystyle
B_{2}$ $\displaystyle=(e^{k_{ep}\tau}-1)A_{2}$

As it can be seen in the definitions shown in equation 14, the parameters are
not independent if the Brix-Tofts model strictly applies and only a set of
four independent parameters is left as discussed in the previous section. In
particular, there is some relationship between the parameters,

$\frac{A_{1}}{A_{2}}=\frac{k_{ep}}{k_{el}}$ (15)

Because the shortest sampling interval in DCE-MRI depends on the volume
acquisition time, the number of points and the precise definition of $\tau$,
which is required for the application of this method, is limited.
Nevertheless, as a first approach, each section of the voxel intensity
evolution, i.e., wash-in or wash-out, can be analyzed separately and the
resulting exponents compared. Let us consider the case of the wash-in section,
which according to equation 2.3.2 contains a time-independent term and can be
eliminated by taking differences between consecutive points,
$q_{i}=p_{i}-p_{i+1}$, for $1\leq i\leq n-1$. The original system, equation
2.3.2 is transformed to [26],

$\left[\begin{array}[]{c}q_{3}\\\
q_{4}\end{array}\right]=\left[\begin{array}[]{cc}q_{2}&-q_{1}\\\
q_{3}&-q_{2}\end{array}\right]\left[\begin{array}[]{c}X_{1}+X_{2}\\\
X_{1}X_{2}\end{array}\right]$ (16)

with solutions $Z_{1}^{*}\equiv X_{1}+X_{2}$, $Z_{2}^{*}\equiv X_{1}X_{2}$.
The solution for $X_{1}$ and $X_{2}$ are obtained by finding the roots of a
second-degree polynomial,

$X^{2}-Z_{1}^{*}X+Z_{2}^{*}=0$ (17)

with the condition that the roots must be real and $0<X<1$. Exponents and
coefficients are determined straightforwardly by substitution in equations
2.3.2 and 14. An analogous procedure could be performed for the wash-out
section and the resulting exponents and coefficients compared to those
determined for the wash-in section. In principle, parameters extracted from
each trend, i.e., wash-in and wash-out, should be equal if a strict Brix-Tofts
model applies. This procedure is somewhat cumbersome and could be affected by
the SNR in the image intensity and patient movement so the matching of
parameters seems to be unlikely. The method is then limited only to the
analysis of the wash-in or wash-out data [31].

A possible solution to this fact is to use the whole set of data points,
assuming a common $\tau$, which can be modeled as the following discretized
set of equations,

$\displaystyle p_{i}=A_{1}(1-X_{1}^{i})-A_{2}(1-X_{2}^{i}),1\leq i\leq
n,i\delta t\leq\tau$ $\displaystyle p_{n+j}=A_{1}(X_{1}^{j}-X_{1}^{n+j})-$
$\displaystyle\qquad\qquad-A_{2}(X_{2}^{j}-X_{2}^{n+j}),1\leq j\leq
m,(n+j)\delta t>\tau$ (18)

which can be written in matrix form as

$\left(\begin{array}[]{c}p_{1}\\\ \vdots\\\ p_{i}\\\ \vdots\\\ p_{n}\\\
p_{n+1}\\\ \vdots\\\ p_{n+j}\\\ \vdots\\\
p_{n+m}\end{array}\right)=\left(\begin{array}[]{cc}(1-X_{1}^{1})&-(1-X_{2}^{1})\\\
\vdots&\vdots\\\ (1-X_{1}^{i})&-(1-X_{2}^{i})\\\ \vdots&\vdots\\\
(1-X_{1}^{n})&-(1-X_{2}^{n})\\\
(X_{1}^{1}-X_{1}^{n+1})&-(X_{2}^{1}-X_{2}^{n+1})\\\ \vdots&\vdots\\\
(X_{1}^{j}-X_{1}^{n+j})&-(X_{2}^{j}-X_{2}^{n+j})\\\ \vdots&\vdots\\\
(X_{1}^{m}-X_{1}^{n+m})&-(X_{2}^{m}-X_{2}^{n+m})\end{array}\right)\left(\begin{array}[]{c}A_{1}\\\
A_{2}\end{array}\right)$ (19)

with $X_{1}>X_{2}$.

Equation 19 resembles a modified version of Vandermonde matrix in terms of
geometric progressions and can be solved by a nonnegative least-square method
for the coefficients $A_{1}$ and $A_{2}$, if a pair $X_{1}$, $X_{2}$ is given.
It is necessary to search the space $(X_{1},X_{2})$ for an optimal solution of
equation 19. To do so these values are selected initially at random within the
square of side 1 and subjected to the condition $X_{1}>X_{2}$, as shown in
Figure 3. The best solution for $A_{1}$, $A_{2}$, $X_{1}$, and $X_{2}$ is
obtained by following the steepest-descent optimization method to minimize the
residuals of equation 19, A possible path is shown in Figure 3.

Figure 3: Location of exponential parameters for the set given by equation 19.
Dashed areas represent possible solutions. In red, the iterative path followed
by a steepest-descent method to obtain an optimal solution for $A_{1}$,
$A_{2}$ of equation 19.

One of the advantages of the Vandermonde matrix formulation is that it allows
for an irregular sampling of the kinetics data as compared to the modified de
Prony formulation which requires regular sampling. This fact allows for the
implementation of DCE-MRI protocols suited to the particular requirements of
pharmacokinetics.

####  2.3.3 Modified de Prony method fit. Three Compartments Model

In the case of a three compartments model, the total concentration in the
tissue is

$C\equiv C_{s}+C_{f}$ (20)

so equations (2.3.1) have to be modified to

$\displaystyle C=$ $\displaystyle A_{1}(1-e^{-k_{el}t})-$ (21)
$\displaystyle\qquad-A_{2}(1-e^{-k_{eps}t})-A_{3}(1-e^{-k_{epf}t}),t\leq\tau$
$\displaystyle C=$ $\displaystyle A_{1}(e^{k_{el}\tau}-1)e^{-k_{el}t}-$
$\displaystyle\qquad\qquad-A_{2}(e^{k_{eps}\tau}-1)e^{-k_{eps}t}-$
$\displaystyle\qquad\qquad\qquad\qquad-
A_{3}(e^{k_{epf}\tau}-1)e^{-k_{epf}t},t>\tau$

and the straightforward modifications to equations (2.3.2) and (19) with the
condition $X_{1}>X_{2}>X_{3}$, which in matrix form is,

$\left(\begin{array}[]{c}p_{1}\\\ \vdots\\\ p_{i}\\\ \vdots\\\ p_{n}\\\
p_{n+1}\\\ \vdots\\\ p_{n+j}\\\ \vdots\\\
p_{n+m}\end{array}\right)=\left(\begin{array}[]{ccc}(1-X_{1}^{1})&-(1-X_{2}^{1}&-(1-X_{3}^{1})\\\
\vdots&\vdots\\\ (1-X_{1}^{i})&-(1-X_{2}^{i})&-(1-X_{3}^{i})\\\
\vdots&\vdots\\\ (1-X_{1}^{n})&-(1-X_{2}^{n})&-(1-X_{3}^{n})\\\
(X_{1}^{1}-X_{1}^{n+1})&-(X_{2}^{1}-X_{2}^{n+1})&-(X_{3}^{1}-X_{3}^{n+1})\\\
\vdots&\vdots\\\
(X_{1}^{j}-X_{1}^{n+j})&-(X_{2}^{j}-X_{2}^{n+j})&-(X_{3}^{j}-X_{3}^{n+j})\\\
\vdots&\vdots\\\
(X_{1}^{m}-X_{1}^{n+m})&-(X_{2}^{m}-X_{2}^{n+m})&-(X_{3}^{m}-X_{3}^{n+m})\end{array}\right)\left(\begin{array}[]{c}A_{1}\\\
A_{2}\\\ A_{3}\end{array}\right)$ (22)

###  2.4 Tissue classification by semi-qualitative methods

Tissue classification has been accomplished in DCE-MRI data in different ways
[7, 8, 9, 10, 11, 12], exploiting properties of the data time evolution, and
has been extensively used clinically. Among them, the so-called three time
point method [7, 8, 10, 12, 32] has an extended use. The method consists of
selecting three temporal points where the time evolution of the contrast agent
uptake is represented appropriately. These points are commonly selected as the
time for the initial image set, corresponding to an absence of contrast, an
intermediate time point, associated with the infusion time $\tau$, and a third
point at the end of the contrast washout region. Image intensity differences
are calculated for each voxel and are used for tissue classification according
to the scheme shown in Figure 4.

Figure 4: Three-time points method for tissue classification in DCE-MRI [8,
9].

Color intensity is determined by the image difference $C_{1}$ also shown in
Figure 4. The selection of the color scheme is made to correlate with
histological measurements [9] assuming that ‘red’ means pathological tissue
with a Type III kinetics, ‘blue’ means not affected tissue with a Type I
kinetics, and ‘green’ applies to those tissues with uncertain conditions, with
a Type II kinetics. If this color code scheme is applied to the Brix-Tofts
pharmacokinetic model it is possible to obtain a color map for the
pharmacokinetic parameters. In the present work, semi-qualitative analysis is
used to classify quantitative results.

##  3 Results and Discussion

###  3.1 Semi-qualitative results

Results obtained with the application of the three-time point method are
partially shown in Figure 5. The prostate was manually segmented and only
points within the ROI were evaluated. There was a gradual color transition for
all the cases, i.e., red to green to blue, from the tumor lesion’s interior to
its periphery, as shown in some examples at the bottom rows in Figure 5.
Spurious points are neglected based on their intensity, i.e., $C_{i}$ value,
and local neighborhood.

Figure 5: Application of the three-time point method. The top row, left side,
superposition of the classified image, and the original one, right side, is
the distribution of points according to the $C_{2}/C_{1}$ ratio. The bottom
rows are examples of the qualitative classification made by the three-time
point method.

Some semi-quantitative information can be extracted from the qualitative
analysis. The semi-quantitative parameters $C_{1}$ and $C_{2}$, defined in
Figure 4 can be used to determine semi-quantitative evaluations that help to
characterize the pathology or changes during a therapy follow-up [33, 34, 35,
36]. An example of this approach is shown in the right-hand side of Figure 5,
where a distribution is obtained for the $C_{2}\textfractionsolidus C_{1}$
ratio.

###  3.2 Quantitative results

####  3.2.1 Levenberg-Marquardt results.

The quantitative analysis of the image data set was performed in a previous
work [37, 31] using only a five-parameter fit of the Brix-Tofts model to the
measured data since as can be seen in Figure 2, it was substantially better
than the four-parameter implementation. Figure 6 it is shown an example of the
results for different points within the prostate, demonstrating the
discrimination capacity of the method to establish significant differences
between the pharmacokinetic parameters $k_{el}$ and $k_{ep}$, on a voxel by
voxel basis.

Figure 6: Results of the Levenberg-Marquardt five parameters fitting at two
different points: top, a point located within the cancerous lesion, a red
region in Figure 5 with fitted parameters $k_{ep}=0.1321s^{-1}$ and
$k_{el}=0.0034s^{-1}$; bottom, located in a region outside the cancerous
lesion, blue-green (Figure 5) with parameters, $k_{ep}=0.0123s^{-1}$, and
$k_{el}=0.00009s^{-1}$

The procedure was extended to determine the distribution of these parameters
over a region of interest like the one represented in Figure 5, as shown in
Figure 7.

Figure 7: Parameter distributions obtained from Levenberg-Marquardt fit

Figure 7 reflects an interesting result that possibly requires further
analysis. It is the apparent difference in point clustering, particularly for
those points associated with tumor lesions (red) and undefined and possibly
infiltrated tissue (green). This difference can be used to further classify
the tumor lesions or their evolution during therapy follow-up. Even though the
Levenberg-Marquardt fit is a very reliable method it depends strongly on the
actual parameter structure of the pharmacokinetic model as previously
discussed and shown in Figure 2, making it unsuitable for models that include
more than two compartments. In the following, we will limit our discussion to
a simpler approach based on the modified de Prony method.

####  3.2.2 Modified de Prony method results. Two Compartments Model

The analysis of the quantitative results based on a Levenberg-Marquardt
indicates the validity of the Brix-Tofts two-compartment model in the sense
that it applies to a specific voxel. When the whole set of voxels is analyzed,
as shown in Figure 7, it is evident that there are two distinct regions for
the kinetic parameters that can be understood as the existence of fast
exchange compartments ($k_{ep}$ high) and slow compartments ($k_{ep}$ low),
allowing for the proposition of a three-compartment model as previously
discussed. The application of the two-compartment model using the modified de
Prony algorithm proposed in this work reveals the same behavior as shown by
the results in Figure 8 and Figure 9

Figure 8: Two compartment model distributions $(k_{ep},k_{el})$. Left, color
map; Middle, $k_{ep}$ distribution; Right, $k_{el}$ distribution. Colors are
given according to the 3-time-point classification Figure 9: Distributions
$(k_{ep},k_{el})$. (a) 2D histograms; (b) $(k_{ep},k_{el})$ plot. Colors are
given according to the 3-time-point classification

In particular, in Figure 8, the distribution for the $k_{ep}$ parameter is bi-
modal with a high exchange region mostly populated by voxels associated with
the presence of tumor (red and green colors) while the slow exchange region is
predominantly composed of voxels that belong to benign tissue (blue color).
This feature can be greatly emphasized if the logarithm of the quotient
$\frac{A_{1}}{A_{2}}$ is plotted against the logarithm of
$\frac{k_{ep}}{k_{el}}$, as shown in Figure 10.

Figure 10: Top: Data trends for $\frac{C_{el}}{C_{ep}}$ (or
$\frac{A_{1}}{A_{2}}$, according to equation 14) and $\frac{k_{ep}}{k_{el}}$
showing how these ratios are correlated. Red straight lines represent the
condition $\frac{C_{el}}{C_{ep}}=\frac{k_{ep}}{k_{el}}$ which corresponds to a
pure two-compartment model. Fast exchange compartments (red dots) are
differentiated from slow exchange compartments (blue dots). Black lines
identify trends in either the fast or slow compartments. Bottom rows: Examples
of tissue classification according to the type of exchange compartments. (F)
Fast exchange compartments ($k_{ep}$ high). (S) Slow exchange compartments
($k_{ep}$ low). Colors are given according to the 3-time-point classification.
Notice that the majority of the fast exchange component is associated with
tumor activity, while the slow exchange component is associated with benign
tissue

If only one type of compartment is available for the tissue, corresponding to
a pure Brix-Tofts model, the plot should behave as a straight line derived
from equation 15 and shown as a red line in Figure 10. The real trend of the
data supports what is observed in the $k_{ep}$ distribution shown in Figure 8
with evident bi-modal behavior. All these observations suggest that the
description of the system should be addressed with a three or more-compartment
model. The spatial distribution of the fast and slow exchange compartments is
shown at the bottom of Figure 10. Specific values of the kinetic parameters
for the different kinetic types as classified by the three-time point method
are summarized in Table 1 and Figure 11.

Table 1: Fast and Slow Exchange Kinetic Parameters Classification | $k_{epf}(min^{-1})$ | $k_{eps}(min^{-1})$ | $k_{el}(min^{-1})$
---|---|---|---
Type I | $5.72\pm 0.32$ | $0.63\pm 0.17$ | $0.04\pm 0.05$
Type II | $5.15\pm 0.25$ | $0.96\pm 0.31$ | $0.05\pm 0.05$
Type III | $5.12\pm 0.24$ | $0.78\pm 0.22$ | $0.12\pm 0.08$
Average values for the kinetic parameters in the regions classified by the
three-time point method. Type I corresponds to benign tissue and Type III to
definitely malignant tissue. Type II corresponds to an intermediate condition,
possibly related to tumor invasion.
Figure 11: Data points in $(k_{epf},k_{eps})$ space. Big circles represent
average values from Table 1. Color-shaded areas represent points within a
standard deviation. Colors are given according to the 3-time-point
classification

It is important to remark that fast exchange compartments are mostly
associated with Type II and Type III kinetics while Type I kinetics
corresponds to slow exchange compartments. This allows for an effective
$\langle k_{ep}\rangle$ defined as a weighted average between both compartment
types. The result is summarized in Table 2.

Table 2: Proportion of fast exchange compartments and average $\langle k_{ep}\rangle$ Classification | $f_{fast}$ | $\langle k_{ep}\rangle(min^{-1})$
---|---|---
Type I | $0.20\pm 0.10$ | $1.75\pm 0.60$
Type II | $0.70\pm 0.22$ | $4.00\pm 0.99$
Type III | $0.77\pm 0.20$ | $4.21\pm 0.97$
Figure 12: Average value $k_{ep}$ for different kinetic types. Colors are
given according to the 3-time-point classification

####  3.2.3 Modified de Prony method results. Three Compartments Model

The previous section discussed the analysis with a pure Brix-Tofts model (two
compartments system) with a clear indication of the presence of fast and slow
exchange compartments. With the assumption of a three-compartment model, as
stated in section 2.3.3, a more detailed picture of the kinetic parameters is
obtained, as shown in Figure 13, with some additional structure for the fast
exchange parameter, $k_{epf}$. Also, in Figure 14 it is shown the distribution
of kinetic parameters in the $(k_{el},k_{eps})$ and $(k_{el},k_{epf})$ spaces.
Nevertheless, comparing the $k_{epf}$ distribution, Figure 13b with the
$k_{ep}$ distribution in Figure 8, they are both bi-modal distributions and
very similar indeed, so it can be concluded that the addition of an extra set
of parameters to the model does not improve the overall picture increasing
notoriously computational times.

Figure 13: Distributions of the kinetic parameters using a three compartments
model. (a) $k_{eps}$ for slow exchange,(b) $k_{epf}$ for fast exchange and (c)
$k_{el}$. Colors are given according to the 3-time-point classification Figure
14: Distributions of the kinetic parameters using a three compartments model
in parameter space. (a) $(k_{eps},k_{el})$ for slow exchange and (b)
$(k_{epf},k_{el})$ for fast exchange. Colors are given according to the
3-time-point classification

## Conclusions

In the present work, it has been shown a general methodology to classify and
quantify tumor lesions in the prostate. The classification of the tissues was
performed successfully by the implementation of a suitable three-time point
method algorithm. This classification was used as a reference for the analysis
of quantitative results. The quantitative analysis was successfully performed
by a novel point-wise method based on the two-compartment Brix-Tofts model and
validated with the Levenberg-Marquardt optimization method. The results
confirmed the existence of fast exchange compartments associated with tumor
activity and slow exchange compartments associated with benign tissue. In any
case, the combined qualitative and quantitative analysis can be used to
establish differences for a patient undergoing therapy or pathology
progression. In particular, some results suggest the future possibility of
combining the qualitative and quantitative analysis into a single
representation for prostate cancer evaluation and therapy follow-up.

## Acknowledgment

The authors would like to thank the National Insitute for Bioengineering,
INABIO, at Universidad Central de Venezuela, Venezuela, Universidad Nacional
Pedro Henríquez Ureña, Dominican Republic and Fundación Arturo López Pérez,
Chile, for providing the environment for the realization of this work. Also,
we would like to express our gratitude to the scientific community at the
Physics and Mathematics in Biomedicine Consortium which provided helpful
discussions during the realization of this work.

## References

  * [1] P. L. Choyke, A. J. Dwyer, and M. V. Knopp, “Functional tumor imaging with dynamic contrast-enhanced magnetic resonance imaging,” J. Magn. Reson. Imaging 17, 509–520 (2003).  [doi:10.1002/jmri.10304].
  * [2] R. Alonzi, A. R. Padhan, and C. Allen, “Dynamic contrast enhanced mri in prostate cancer,” Eur. J. Radiol 63, 335–350 (2007).  [doi:10.1016/j.ejrad.2007.06.028].
  * [3] M. Haris, R. K. Gupta, A. Singh, et al., “Differentiation of infective from neoplastic brain lesions by dynamic contrast-enhanced mri,” Neuroradiology 50, 531–540 (2008).  [doi:10.1007/s00234-008-0378-6].
  * [4] C. Foottit, G. O. Cron, M. J. Hogan, et al., “Determination of the venous output function from mr signal phase: Feasibility for quantitative dce-mri in human brain,” Magn. Reson. Med. 63, 772–781 (2010).  [doi:10.1002/mrm.22253].
  * [5] S. Verma, B. Turkbey, N. Muradyan, et al., “Overview of dynamic contrast- enhanced mri in prostate cancer diagnosis and management,” Am. J. Roentgenol. 198, 1277–1288 (2012).  [doi:10.2214/AJR.12.8510].
  * [6] B.-B. Chen and T.-F. Shih, “Dce-mri in hepatocellular carcinomaclinical and therapeutic image biomarker,” World J. Gastroenterol. 20, 3125–3134 (2014).  [doi:10.3748/wjg.v20.i12.3125].
  * [7] H. Degani, V. Gusis, D. Weinstein, et al., “Mapping pathophysiological features of breast tumors by mri at high spatial resolution,” Nat. Med. 3, 780–782 (1997).  [doi:10.1038/nm0797-780].
  * [8] E. A. M. Hauth, C. Stockamp, S. Maderwald, et al., “Evaluation of the three-time-point method for diagnosis of breast lesions in contrast-enhanced mr mammography,” Clin. Imaging 30, 160–165 (2006).  [doi:10.1016/j.clinimag.2005.11.005].
  * [9] R. E. Lenkinski, “The role of magnetic resonance imaging in prostate cancer imaging and staging,” in Proceedings 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2109–2112 (2011).  [doi:10.1109/ISBI.2011.5872829].
  * [10] J. U. Fluckiger, M. C. Schabel, and E. V. R. DiBella, “The effect of temporal sampling on quantitative pharmacokinetic and three-time-point analysis of breast dce-mri,” Magn. Reson. Imaging 30, 934–943 (2012).  [doi:10.1016/j.mri.2012.02.011].
  * [11] C. Lavini, M. Buiter, and M. Maas, “Use of dynamic contrast enhanced time intensity curve shape analysis in mri: theory and practice,” Rep. Med. Imag. 6, 71–82 (2013).
  * [12] E. Furman-Haran, M. S. Feinberg, D. Badhiki, et al., “Standardization of radiological evaluation of dynamic contrast enhanced mri: Application in breast cancer diagnosis,” Technol. Cancer Res. Treat. 13, 445–454 (2014).  [doi:10.7785/tcrt.2013.600263].
  * [13] S. Sourbron, M. Ingrisch, A. Siefert, et al., “Quantification of cerebral blood flow, cerebral blood volume, and blood–brain-barrier leakage with dce-mri,” Magn. Reson. Med. 62, 205–217 (2009).  [doi:10.1002/mrm.22005].
  * [14] T. S. Koh, S. Bisdas, D. M. Koh, et al., “Fundamentals of tracer kinetics for dynamic contrast-enhanced mri,” J. Magn. Reson. Imaging 34, 1262–1276 (2011).  [doi:10.1002/jmri.22795].
  * [15] G. Brix, W. Semmler, R. Port, et al., “Pharmacokinetic parameters in cns gd-dtpa enhanced mr imaging,” J. Comput. Assist. Tomogr. 15, 621–628 (1991).  [doi:10.1097/00004728-199107000-00018].
  * [16] P. S. Tofts and A. G. Kermode, “Measurement of the blood-brain barrier permeability and leakage space using dynamic mr imaging. 1.fundamental concepts,” Magn. Reson. Med. 17, 357–367 (1991).  [doi:10.1002/mrm.1910170208].
  * [17] P. S. Tofts, “Modeling tracer kinetics in dynamic gd-dtpa mr imaging,” J. Magn. Reson. Imaging 7, 91–101 (1997).  [doi:10.1002/jmri.1880070113].
  * [18] R. Fusco, M. Sansone, S. Maffei, et al., “Dynamic contrast-enhanced mri in breast cancer: A comparison between distributed and compartmental tracer kinetic models,” J. Biomed. Graph. Comput. 2, 23–36 (2012).  [doi:10.5430/jbgc.v2n2p23].
  * [19] A. Ambikapathi, T.-H. Chan, C.-H. Lin, et al., “Convex-optimization-based compartmental pharmacokinetic analysis for prostate tumor characterization using dce-mri,” IEEE Trans. Biomed. Eng. 63, 707–720 (2016).  [doi:10.1109/TBME.2015.2469601].
  * [20] L. Chen, P. L. Choyke, T.-H. Chan, et al., “Tissue-specific compartmental analysis for dynamic contrast-enhanced mr imaging of complex tumors,” IEEE Trans. Med. Imag. 30, 2044–2058 (2011).  [doi:10.1109/TMI.2011.2160276].
  * [21] L. Chen, T.-H. Chan, P. L. Choyke, et al., “Camcm: a signal deconvolution tool for in vivo dynamic contrast-enhanced imaging of complex tissues,” Bioinformatics 27, 2607–2609 (2011).  [doi:10.1093/bioinformatics/btr436].
  * [22] Y. Wang, J. Xuan, R. Srikanchana, et al., “Modeling and reconstruction of mixed functional and molecular patterns,” Int. J. Biomed. Imag. 2006.29707, 1–9 (2006).  [doi:10.1155/IJBI/2006/29707].
  * [23] K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
  * [24] D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” SIAM J. Appl. Math. 11, 431–441 (1963).  [doi:10.1137/0111030].
  * [25] G. R. de Prony, “Essai experimental et analytique : sur les lois de la dilatabilite des fluides elastique et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alkool, a differentes temperatures,” Journal de l’E´cole Polytechnique 1, 24 (1795).
  * [26] M. Martín-Landrove, G. Figueroa, M. Paluszny, et al., “A quasi- analytical method for relaxation rate distribution determination of t2-weighted mri in brain,” in Proceedings of the 29th Annual Conference of the IEEE EMBS, 1318–1321 (2007).  [doi:10.1109/IEMBS.2007.4352540].
  * [27] K. Clark, B. Vendt, K. Smith, et al., “The cancer imaging archive (tcia): Maintaining and operating a public information repository,” J. Digit. Imaging 26, 1045–1057 (2013).  [doi:10.1007/s10278-013-9622-7].
  * [28] D. Mattes, D. R. Haynor, H. Vesselle, et al., “Nonrigid multimodality image registration,” in Proceedings of the SPIE, Medical Imaging 2001, 1609–1620 (2001).  [doi:10.1117/12.431046].
  * [29] J. Nocedal and S. Wright, Numerical Optimization. Springer Series in Operation Research and Financial Engineering, Springer-Verlag (2006).
  * [30] I. M. J. van der Bom, S. Klein, M. Staring, et al., “Evaluation of optimization methods for intensity based 2d-3d registration in x-ray guided interventions,” in Proceedings of the SPIE, Medical Imaging 2011, 7962231–79622315 (2011).  [doi:10.1117/12.877655].
  * [31] J. Aponte, “Pharmacokinetic parameters quantification of dynamic contrast enhanced mri of prostate using a variation of the prony method,” Master’s thesis, Universidad Central de Venezuela, Caracas, Venezuela (2016).
  * [32] M. Sansone, R. Fusco, and A. Petrillo, “A geometrical perspective on the 3tp method in dce-mri,” Biomedical Signal Processing and Control 16, 32–39 (2015).  [doi:10.1016/j.bspc.2014.10.004].
  * [33] S. Wu, B. F. Kurland, W. A. Berg, et al., “Signal enhancement ratio (ser) quantified from breast dce-mri and breast cancer risk,” in Proceedings of the SPIE,Medical Imaging 2015, 94140M1–94140M6 (2015).  [doi:10.1117/12.2082466].
  * [34] J. Xiao, H. Rahbar, D. S. Hippe, et al., “Dynamic contrast-enhanced breast mri features correlate with invasive breast cancer angiogenesis,” npj Breast Cancer 7, 421–429 (2021).  [doi:10.1038/s41523-021-00247-3].
  * [35] J. Zeng, Q. Cheng, D. Zhang, et al., “Diagnostic ability of dynamic contrast-enhanced magnetic resonance imaging for prostate cancer and clinically significant prostate cancer in equivocal lesions: A systematic review and meta-analysis,” Frontiers in Oncology 11, 6206281–62062812 (2021).  [doi:10.3389/fonc.2021.620628].
  * [36] J. Zhu, J. Yun, K. Wang, et al., “Assessment of early response to lung cancer chemotherapy by semiquantitative analysis of dynamic contrast-enhanced mri,” Disease Markers 2022, 26692811–26692817 (2022).  [doi:10.1155/2022/2669281].
  * [37] A. Ruiz, “Pharmacokinetic parameters quantification in dinamic contrast enhanced mri for prostate cancer,” Master’s thesis, Universidad Central de Venezuela, Caracas, Venezuela (2017).

###### List of Figures

  1. 1 Schematic representation of e two-compartment Brix-Tofts model
  2. 2 Different data fitting results for the same pharmacokinetic model, fitted curve is indicated in red. (a) Four parameters, (b) Five parameters. In both cases, the Levenberg Marquardt algorithm was used
  3. 3 Location of exponential parameters for the set given by equation 19. Dashed areas represent possible solutions. In red, the iterative path followed by a steepest-descent method to obtain an optimal solution for $A_{1}$, $A_{2}$ of equation 19.
  4. 4 Three-time points method for tissue classification in DCE-MRI [8, 9].
  5. 5 Application of the three-time point method. The top row, left side, superposition of the classified image, and the original one, right side, is the distribution of points according to the $C_{2}/C_{1}$ ratio. The bottom rows are examples of the qualitative classification made by the three-time point method.
  6. 6 Results of the Levenberg-Marquardt five parameters fitting at two different points: top, a point located within the cancerous lesion, a red region in Figure 5 with fitted parameters $k_{ep}=0.1321s^{-1}$ and $k_{el}=0.0034s^{-1}$; bottom, located in a region outside the cancerous lesion, blue-green (Figure 5) with parameters, $k_{ep}=0.0123s^{-1}$, and $k_{el}=0.00009s^{-1}$
  7. 7 Parameter distributions obtained from Levenberg-Marquardt fit
  8. 8 Two compartment model distributions $(k_{ep},k_{el})$. Left, color map; Middle, $k_{ep}$ distribution; Right, $k_{el}$ distribution. Colors are given according to the 3-time-point classification
  9. 9 Distributions $(k_{ep},k_{el})$. (a) 2D histograms; (b) $(k_{ep},k_{el})$ plot. Colors are given according to the 3-time-point classification
  10. 10 Top: Data trends for $\frac{C_{el}}{C_{ep}}$ (or $\frac{A_{1}}{A_{2}}$, according to equation 14) and $\frac{k_{ep}}{k_{el}}$ showing how these ratios are correlated. Red straight lines represent the condition $\frac{C_{el}}{C_{ep}}=\frac{k_{ep}}{k_{el}}$ which corresponds to a pure two-compartment model. Fast exchange compartments (red dots) are differentiated from slow exchange compartments (blue dots). Black lines identify trends in either the fast or slow compartments. Bottom rows: Examples of tissue classification according to the type of exchange compartments. (F) Fast exchange compartments ($k_{ep}$ high). (S) Slow exchange compartments ($k_{ep}$ low). Colors are given according to the 3-time-point classification. Notice that the majority of the fast exchange component is associated with tumor activity, while the slow exchange component is associated with benign tissue
  11. 11 Data points in $(k_{epf},k_{eps})$ space. Big circles represent average values from Table 1. Color-shaded areas represent points within a standard deviation. Colors are given according to the 3-time-point classification
  12. 12 Average value $k_{ep}$ for different kinetic types. Colors are given according to the 3-time-point classification
  13. 13 Distributions of the kinetic parameters using a three compartments model. (a) $k_{eps}$ for slow exchange,(b) $k_{epf}$ for fast exchange and (c) $k_{el}$. Colors are given according to the 3-time-point classification
  14. 14 Distributions of the kinetic parameters using a three compartments model in parameter space. (a) $(k_{eps},k_{el})$ for slow exchange and (b) $(k_{epf},k_{el})$ for fast exchange. Colors are given according to the 3-time-point classification

###### List of Tables

  1. 1 Fast and Slow Exchange Kinetic Parameters
  2. 2 Proportion of fast exchange compartments and average $\langle k_{ep}\rangle$