# Is there a trade-off relation between efficiency and power in a collisional
Penrose process in an extreme Reissner-Nordström spacetime?

Naoki Tsukamoto1<EMAIL_ADDRESS>1Department of General Science and
Education, National Institute of Technology, Hachinohe College, Aomori
039-1192, Japan

###### Abstract

We investigate the power of a collisional Penrose process with an unbound
energy extraction from an extreme Reissner-Nordström black hole. This process
takes infinite time in a time coordinate at a constant radial coordinate
outside of the black hole. For black holes as a power plant, the power of the
process for an observer far away from the black hole can be useful. We define
the power as energy gain from the extreme Reissner-Nordström black hole
divided by the time interval of the process in a coordinate time; we estimate
the upper bound of the power in a near-horizon limit, while the efficiency of
the process can be arbitrary large in that limit. Thus, we conclude that there
is no trade-off relation between the efficiency and power in the collisional
Penrose process in extreme Reissner-Nordström spacetime.

##  I Introduction

In 1969, Penrose suggested a process, called the Penrose process, to extract
rotational energy from a Kerr black hole by dropping a particle with negative
energy in the black hole Penrose:1969pc . Denardo and Ruffini pointed out that
the electromagnetic counterpart of the Penrose process occurs in Reissner-
Nordström spacetime Denardo:1973pyo , and Blandford and Znajek considered the
electromagnetic extraction of rotational energy from the Kerr black hole in an
astrophysical situation Blandford:1977ds . It is considered that the
electromagnetic extraction of energy from the black holes can explain high
energy jets near them.

Piran, Shaham, and Katz Piran:1975  and T. Piran and J. Shaham Piran:1977dm
investigated a particle collision near the Kerr black hole. They pointed out
that the center-of-mass energy of two particles can be arbitrarily high if the
Kerr black hole has an extreme event horizon and one of two particles has a
critical angular momentum. Particle collision with infinite center-of-mass
energy is called the Bañados-Silk-West process since it was rediscovered by
Bañados, Silk, and West in 2009 Banados:2009pr .

After rediscovering the BSW process, several authors critically investigated
the process with arbitrarily high center-of-mass energy: (i) The Kerr black
hole cannot have an extreme rotation in an astrophysical situation
Berti:2009bk , (ii) gravitational radiation and its backreaction constrain the
center-of-mass energy for the particle collision Berti:2009bk , (iii) to
obtain the unbounded center-of-mass energy, an infinite proper time of a
falling particle to reach the extreme event horizon is required
Jacobson:2009zg , and (iv) self-gravity of falling objects constrains the
center-of-mass energy for the collision Kimura:2010qy ; Ogasawara:2018gni .

Reference Banados:2009pr  inspired several authors to investigate the details
of the collisional Penrose process Piran:1977dm . The upper limit of energy
extraction by the collision Penrose process after the BSW collision near the
extreme Kerr black hole is very modest Bejger:2012yb ; Harada:2012ap .
Schnittman found that a collisional Penrose process with energy gain can be
more than $10$ times the energy of incident particles Schnittman:2014zsa ;
Leiderschneider:2015kwa ; Ogasawara:2015umo ; Harada:2016eff ;
Zaslavskii:2016unn . A collisional Penrose process after the head-on collision
of two particles near an extreme Kerr black hole Berti:2014lva ;
Leiderschneider:2015ika ; Zaslavskii:2015ema ; Ogasawara:2016yfk , collisional
Penrose processes with spinning particles Maeda:2018hfi ; Okabayashi:2019wjs ,
and collisional Penrose processes in wormhole spacetimes Tsukamoto:2015hta ;
Zaslavskii:2018kix  and in an overspinning Kerr spacetime Patil:2015fua  were
investigated.

Zaslavskii found the electromagnetic counterpart of BSW collision near an
extreme charged Reissner-Nordström black hole Zaslavskii:2010aw . Zaslavskii
pointed out that energy extraction from an extreme Reissner-Nordström black
hole in a collisional Penrose process after the BSW collision can be unbound
Zaslavskii:2012ax ; Nemoto:2012cq  while the energy extraction from the
extreme Kerr black hole is very modest 111Bhat, Dhurandhar, and Dadhich found
no upper limit on the efficiency of a Penrose process with electromagnetic
interaction in a Kerr-Newman black hole spacetime in Ref. Dhurandhar_1985 ..
The Reissner-Nordström spacetime is more tractable than the Kerr black hole
due to spherical symmetry of the spacetime. A finite center-of-mass energy of
BSW collisions of two shells including their self-gravity was shown in Ref.
Kimura:2010qy , and an upper bound of energy extraction from the extreme
Reissner-Nordström black hole, by fully taking into account the self-gravity
of the colliding shells, was obtained in Ref. Nakao:2017xwe  222Multiple BSW
collisions and multiple Penrose processes in Reissner-Nordström spacetime
Kokubu:2020jvd ; Kokubu:2021cwj  and a BSW collision in higher-dimensional
Reissner-Nordström spacetime Tsukamoto:2013dna  were also investigated..

In this paper, inspired by a trade-off relation between efficiency and power
of a heat engine by Shiraishi, Saito, and Tasaki Shiraishi , we consider the
power of the collisional Penrose process, which gives infinite efficiency in
the extreme Reissner-Nordström spacetime after the BSW collision. Is there a
trade-off relation between efficiency and power of the collisional Penrose
process? To answer this question, we define the power of the collisional
Penrose process as the energy extraction divided by a coordinate time, and we
discuss the maximum of the power in the process.

This paper is organized as follows. In Sec. II, we review the motion of a
charged particle in an extreme Reissner-Nordström spacetime. In Sec. III, we
review the energy extraction from the extreme Reissner-Nordström black hole in
a collision Penrose process and we investigate the power of the process. In
Sec. IV, we conclude with our results. In this paper, we use geometrical units
in which the light speed and Newton’s constant are unity.

##  II Motion of a charged particle in a Reissner-Nordström spacetime

A line element and a vector potential in a Reissner-Nordström spacetime are
expressed by

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2}),$ (1) $\displaystyle A_{\mu}dx^{\mu}=-\frac{Q}{r}dt,$ (2)

where $f(r)$ is given by

$\displaystyle f(r)\equiv 1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}},$ (3)

and $Q$ and $M$ are electrical charge and mass, respectively. We find a black
hole spacetime with an event horizon at $r=r_{\mathrm{H}}\equiv
M+\sqrt{M^{2}-Q^{2}}$ for $\left|Q\right|\leq M$ and a spacetime with naked
singularity for $\left|Q\right|>M$. We assume an extreme charge $Q=M>0$ since
we are interested in a collisional Penrose process with an unbound energy
extraction Zaslavskii:2012ax ; Nemoto:2012cq .

From the Hamiltonian equation, the four-momentum $p^{\mu}$ of a particle with
an electrical charge $q$ is expressed as

$\displaystyle p^{\mu}$ $\displaystyle=$ $\displaystyle\frac{\partial
H}{\partial\pi_{\mu}}$ (4) $\displaystyle=$ $\displaystyle\pi^{\mu}-qA^{\mu},$

where $H$ is the Hamiltonian of the charged particle given by

$\displaystyle
H\equiv\frac{1}{2}g^{\mu\nu}(\pi_{\mu}-qA_{\mu})(\pi_{\nu}-qA_{\nu}),$ (5)

and $\pi^{\mu}$ is the canonical momentum of the charged particle conjugate to
the coordinates $x^{\mu}$. We assume that charged particles have vanishing
angular momentum $L\equiv\pi_{\phi}=0$ and that they only move in a radial
direction on an equatorial plane $\theta=\pi/2$.

From the $t$ component of the four-momentum $p^{\mu}=dx^{\mu}/d\lambda$, where
$\lambda$ is an affine parameter, and from Eq. (4), we get

$\displaystyle\frac{dt}{d\lambda}=\frac{1}{f}\left(E-\frac{qM}{r}\right),$ (6)

where $E\equiv-\pi_{t}$ is the conserved energy of the charged particle. The
particle should satisfy a forward-in-time condition $dt/d\lambda\geq 0$. The
condition is expressed by

$\displaystyle E-\frac{qM}{r}\geq 0$ (7)

and it yields, for $r=r_{\mathrm{H}}=M$,

$\displaystyle E=q.$ (8)

We call a charged particle with $E=q$ a critical particle.

From $g_{\mu\nu}p^{\mu}p^{\nu}=-m^{2}$, where $m$ is the mass of the charged
particle, we obtain

$\displaystyle\left(\frac{dr}{d\lambda}\right)^{2}+V(r)=0,$ (9)

where $V(r)$ is the effective potential of the radial motion of the charged
particle given by

$\displaystyle V(r)\equiv-\left(E-\frac{qM}{r}\right)^{2}+m^{2}f.$ (10)

The charged particle can exist only in a region where the effective potential
$V(r)$ is nonpositive. The radial component of the four-momentum of a particle
can be written as $p^{r}=\sigma\sqrt{-V}$, where $\sigma=-1$ ($\sigma=1$) for
an ingoing (outgoing) particle.

##  III Particle collision and energy extraction from a black hole

In this section, we review energy extraction from the extreme Reissner-
Nordström black hole in a collision Penrose process Zaslavskii:2012ax ;
Nemoto:2012cq  as well as investigate power in this process. We consider that
particles $1$ and $2$ collide at $r=r_{\mathrm{c}}\equiv M(1+\epsilon)$, where
$0<\epsilon\ll 1$, and particles $3$ and $4$ are produced after the collision.
We set $\sigma_{1}=\sigma_{2}=-1$. Here and hereinafter, physical values with
subscripts $1$, $2$, $3$, and $4$ denote physical values of particles $1$,
$2$, $3$, and $4$, respectively. The center-of-mass energy $E_{\mathrm{CM}}$
of particles $1$ and $2$ at the collision is given by

$\displaystyle E_{\mathrm{CM}}^{2}$ $\displaystyle\equiv$ $\displaystyle-
g_{\mu\nu}(p^{\mu}_{1}+p^{\mu}_{2})(p^{\nu}_{1}+p^{\nu}_{2})$ (11)
$\displaystyle=$ $\displaystyle
m^{2}_{1}+m^{2}_{2}+\frac{2}{f(r_{\mathrm{c}})}\left[\left(E_{1}-\frac{q_{1}M}{r_{\mathrm{c}}}\right)\left(E_{2}-\frac{q_{2}M}{r_{\mathrm{c}}}\right)\right.$
$\displaystyle\left.+\sqrt{-V_{1}(r_{\mathrm{c}})}\sqrt{-V_{2}(r_{\mathrm{c}})}\right].$

If one of the particles is critical and the other is not critical, the center-
of-mass energy diverges in a near-horizon limit $\epsilon\rightarrow 0$
Zaslavskii:2010aw . For simplicity, we assume that particle $1$ is critical,
$E_{1}=q_{1}$, and particle $2$ has no charge, $q_{2}=0$. The center-of-mass
energy is obtained as

$\displaystyle E_{\mathrm{CM}}^{2}$ $\displaystyle=$ $\displaystyle
m^{2}_{1}+m^{2}_{2}+\frac{2(1+\epsilon)}{\epsilon}\left[E_{1}E_{2}\right.$
(12)
$\displaystyle\left.-\sqrt{E_{1}^{2}-m_{1}^{2}}\sqrt{E_{2}^{2}-m_{2}^{2}\left(\frac{\epsilon}{1+\epsilon}\right)^{2}}\right]$
$\displaystyle\sim$ $\displaystyle\frac{2A_{1}E_{2}}{\epsilon},$

where $A_{1}$ is defined by $A_{1}\equiv E_{1}-\sqrt{E_{1}^{2}-m_{1}^{2}}$.

The conservation law of the charges before and after the particle collision is
expressed by

$\displaystyle q_{1}+q_{2}=q_{3}+q_{4}.$ (13)

The conservation law of the four-momentum of the particles at the moment of
the collision is expressed by

$\displaystyle p^{\mu}_{1}+p^{\mu}_{2}=p^{\mu}_{3}+p^{\mu}_{4}.$ (14)

The $t$ component of the conservation laws of the four-momentum (14) and the
charges (13) gives the conservation laws of conserved energy,

$\displaystyle E_{1}+E_{2}=E_{3}+E_{4}.$ (15)

We assume that particle $3$ is a near-critical particle with
$q_{3}=E_{3}(1+\delta_{3}\epsilon)$, where $0<\delta_{3}<1$, and that
particles $3$ and $4$ are ingoing particles with $\sigma_{3}=\sigma_{4}=-1$
immediately after the creation of the particles Zaslavskii:2012ax ;
Nemoto:2012cq . From the radial component of Eq. (14), we obtain

$\displaystyle
A_{1}+E_{3}(\delta_{3}-1)=-\sqrt{E_{3}^{2}(1-\delta_{3})^{2}-m_{3}^{2}}.$ (16)

From the square of Eq. (16), we get

$\displaystyle\delta_{3}=1-\frac{m_{3}^{2}+A_{1}^{2}}{2A_{1}E_{3}}.$ (17)

We assume that particle $3$ is reflected at a turning point $r=r_{-}$, and we
assume that $E_{3}>m_{3}$ so that particle $3$ goes to spatial infinity. For
simplicity, we also assume that particle $4$ does not interact with particle
$3$ after their particle production and particle $4$ falls into the event
horizon of the black hole 333One may be concerned that particle $4$ collides
with particle $3$. By taking nonzero conserved angular momenta of particles
$3$ and $4$ into account, the collision can be avoided.. From
$V_{3}(r_{-})=0$, we obtain the turning point $r=r_{-}$ as

$\displaystyle
r_{-}=M\left(1+\frac{E_{3}\delta_{3}\epsilon}{E_{3}-m_{3}}\right).$ (18)

From $r_{-}\leq r_{\mathrm{c}}$, the condition

$\displaystyle E_{3}(1-\delta_{3})\geq m_{3}$ (19)

must be satisfied. If the inequality (19) holds, the inside of the square root
in (16) is positive, and the lower bound of the mass of particle $3$,

$\displaystyle m_{3}\geq A_{1}$ (20)

must hold so that the left-hand side of Eq. (16) is negative. For simplicity,
we assume $m_{0}\equiv m_{1}=m_{2}$ and $E_{1}=E_{2}\gtrsim m_{0}$. As
discussed in Ref. Nemoto:2012cq , the mass and conserved energy of particle
$3$ can be $m_{3}\sim E_{\mathrm{CM}}\sim m_{0}/\sqrt{\epsilon}$ and
$E_{3}\sim m_{0}/\epsilon$, respectively.

The collisional Penrose process is expressed in Fig. 1: Particle $1$ at an
initial position $r=r_{\mathrm{i}}$ falls toward the extreme charged black
hole and it collides with particle $2$ at $r=r_{\mathrm{c}}$ after the
coordinate time of $\Delta t_{1}$. Particles $3$ and $4$ are produced from the
collision. Particle $3$ is reflected at $r=r_{-}$, and it escapes to
$r=r_{\mathrm{i}}$ after the coordinate time of $\Delta t_{3}$. On the other
hand, particle $4$ reaches the extreme event horizon at $r=r_{\mathrm{H}}$.

Figure 1: Collisional Penrose process. Particle $1$ with the critical
condition $E_{1}=q_{1}$ at an initial position $r=r_{\mathrm{i}}$ moves toward
the extreme Reissner-Nordström black hole with $M=Q$ to collide with particle
$2$ with no charge, $q_{2}=0$, at $r=r_{\mathrm{c}}=M(1+\epsilon)$ after the
coordinate time of $\Delta t_{1}$. The collision between particles $1$ and $2$
produces particle $3$, with mass $m_{3}\sim m_{0}/\sqrt{\epsilon}$ and
conserved energy $E_{3}\sim m_{0}/\epsilon$, and particle $4$. Particle $3$ is
reflected at $r=r_{-}=M\left[1+E_{3}\delta_{3}\epsilon/(E_{3}-m_{3})\right]$,
and it reaches $r=r_{\mathrm{i}}$ after the coordinate time of $\Delta t_{3}$.
Particle $4$ falls into the extreme event horizon at $r=r_{\mathrm{H}}=M$.

The time intervals $\Delta t_{1}$ and $\Delta t_{3}$ in the coordinate time
are given by

$\displaystyle\Delta t_{1}$ $\displaystyle=$
$\displaystyle-\frac{E_{1}}{\sqrt{E_{1}^{2}-m_{0}^{2}}}\int^{r_{\mathrm{c}}}_{r_{\mathrm{i}}}\frac{dr}{\left(1-\frac{M}{r}\right)^{2}}$
(21) $\displaystyle\sim$
$\displaystyle\frac{E_{1}}{\sqrt{E_{1}^{2}-m_{0}^{2}}}\left(\frac{M}{\epsilon}+r_{\mathrm{i}}\right)$

and

$\displaystyle\Delta t_{3}$ $\displaystyle=$
$\displaystyle-\int^{r_{-}}_{r_{\mathrm{c}}}\frac{E_{3}-\frac{q_{3}M}{r}}{\left(1-\frac{M}{r}\right)^{2}\sqrt{-V_{3}}}dr$
(22)
$\displaystyle+\int^{r_{\mathrm{i}}}_{r_{-}}\frac{E_{3}-\frac{q_{3}M}{r}}{\left(1-\frac{M}{r}\right)^{2}\sqrt{-V_{3}}}dr$
$\displaystyle\sim$
$\displaystyle\left(\frac{2}{\delta_{3}}-2\right)\frac{M}{\epsilon}+r_{\mathrm{i}},$

respectively.

We define the power of the collisional Penrose process as $\Delta E/\Delta t$,
where $\Delta E$ and $\Delta t$ are defined by $\Delta E\equiv
E_{3}-E_{1}-E_{2}\sim m_{0}/\epsilon$ and $\Delta t\equiv\Delta t_{1}+\Delta
t_{3}$, respectively. From Eqs. (21) and (22), the power of the collisional
Penrose process is given by

$\displaystyle\frac{\Delta E}{\Delta
t}\sim\frac{m_{0}}{\left(\frac{E_{1}}{\sqrt{E_{1}^{2}-m_{0}^{2}}}+\frac{2}{\delta_{3}}-2\right)M+\left(\frac{E_{1}}{\sqrt{E_{1}^{2}-m_{0}^{2}}}+1\right)r_{\mathrm{i}}\epsilon}$

The power is plotted in Fig. 2.

Figure 2: Power of the collisional Penrose process given by Eq. (III). We set
$M=10$ km, $m_{0}=10^{-8}$ km, $E_{1}=1.1m_{0}$, and $r_{\mathrm{i}}=10^{7}$
km.

If $\epsilon\ll M/r_{\mathrm{i}}<1$ is satisfied, by using $E_{3}\sim
m_{0}/\epsilon$, $m_{3}\sim m_{0}/\sqrt{\epsilon}$, and Eq. (17), we get
$\delta_{3}\sim 1-m_{0}/(2A_{1})\sim 1-(e_{1}+\sqrt{e_{1}^{2}-1})/2$, where
$e_{1}$ is a specific conserved energy $e_{1}\equiv E_{1}/m_{0}$ of particle
$1$, and the power is given by

$\displaystyle\frac{\Delta E}{\Delta
t}\sim\frac{1-e_{1}^{2}+(2-e_{1})\sqrt{e_{1}^{2}-1}}{e_{1}^{2}+2e_{1}-2+e_{1}\sqrt{e_{1}^{2}-1}}\frac{m_{0}}{M}.$
(24)

We notice that the collisional Penrose process has the upper bound of the
power, the value of which is estimated to be as given in Eq. (24) in the near-
horizon limit $\epsilon\rightarrow 0$.

##  IV Conclusion

We have defined the power of the collisional Penrose process as the energy
gain from the extreme Reissner-Nordström black hole divided by the time
interval of the process in a coordinate time. The upper bound of the power can
be estimated as given in Eq. (24) in the near-horizon limit
$\epsilon\rightarrow 0$. On the other hand, the efficiency of the collisional
Penrose process defined by $\eta\equiv E_{3}/(E_{1}+E_{2})$ has an arbitrarily
high value $\eta\sim 1/(2\epsilon)$ in the near-horizon limit
$\epsilon\rightarrow 0$ under our assumptions as discussed in Ref.
Nemoto:2012cq . Therefore, we conclude that there is no trade-off relation
between efficiency and power in the collisional Penrose process in the extreme
Reissner-Nordström spacetime under our treatment. One may, however, find a
trade-off relation between efficiency and power if the effect of the self-
gravity of falling particles is taken into account. The power of the
collisional Penrose process including self-gravity of the particles is left
for future work.

We now comment on the upper bound of the power in the collisional Penrose
process. Zaslavskii categorized the scenarios of the collisional Penrose
process with a single collision in the vicinity of the extreme event horizon
to find unbounded energy extraction Zaslavskii:2012ax . All the scenarios,
except for the one shown on Fig. 1, give vanishing power of the process in the
near-horizon limit. However, one may find greater power than Eq. (III) if some
of the conditions on the extreme charge of the black hole, the critical charge
of particles, and the near-horizon collision are violated. Notice also that
the upper bound of the power (III) is not exact due to approximations and
simplifications. In addition, the power of the collisional Penrose processes
can be enhanced if we consider the collision of spinning particles
Armaza:2015eha ; Maeda:2018hfi ; Okabayashi:2019wjs . On the other hand, if
the self-gravity of falling objects is taken into account, the power would be
suppressed Nakao:2017xwe .

The secondary collisional Penrose process of a head-on collision Piran:1977dm
; Zaslavskii:2013nra ; Tsukamoto:2019ihj  between particle $3$, which is
reflected near the black hole, and an additional falling particle may increase
the power of the total process even though the self-gravity of particles $3$
and $4$ would decrease the power in the secondary collision. Investigating the
details of multiple collisional processes Kokubu:2020jvd ; Kokubu:2021cwj  can
also be an interesting future work to find a larger power than Eq. (III).

One may consider that the definition of the power of the process using the
proper time of the particles is more desirable than one using the coordinate
time. The effect of the difference of the definitions on the upper limit of
the power is also left to future work.

## Acknowledgements

The author thanks an anonymous referee for useful comments.

## References

  * (1) R. Penrose, Riv. Nuovo Cim. 1, 252-276 (1969).
  * (2) G. Denardo and R. Ruffini, Phys. Lett. B 45, 259-262 (1973).
  * (3) R. D. Blandford and R. L. Znajek, Mon. Not. Roy. Astron. Soc. 179, 433-456 (1977)
  * (4) T. Piran, J. Shaham, and J. Katz, Astrophys. J. 196 (1975) L107.
  * (5) T. Piran and J. Shaham, Phys. Rev. D 16 (1977), 1615-1635.
  * (6) M. Banados, J. Silk and S. M. West, Phys. Rev. Lett. 103 (2009), 111102.
  * (7) E. Berti, V. Cardoso, L. Gualtieri, F. Pretorius, and U. Sperhake, Phys. Rev. Lett. 103, 239001 (2009).
  * (8) T. Jacobson and T. P. Sotiriou, Phys. Rev. Lett. 104, 021101 (2010).
  * (9) M. Kimura, K. i. Nakao and H. Tagoshi, Phys. Rev. D 83 (2011), 044013.
  * (10) K. Ogasawara and N. Tsukamoto, Phys. Rev. D 99, 024016 (2019).
  * (11) M. Bejger, T. Piran, M. Abramowicz, and F. Hakanson, Phys. Rev. Lett. 109, 121101 (2012).
  * (12) T. Harada, H. Nemoto and U. Miyamoto, Phys. Rev. D 86, 024027 (2012) [erratum: Phys. Rev. D 86, 069902 (2012)].
  * (13) J. D. Schnittman, Phys. Rev. Lett. 113, 261102 (2014).
  * (14) E. Leiderschneider and T. Piran, Phys. Rev. D 93, 043015 (2016).
  * (15) K. Ogasawara, T. Harada and U. Miyamoto, Phys. Rev. D 93, 044054 (2016).
  * (16) T. Harada, K. Ogasawara, and U. Miyamoto, Phys. Rev. D 94, 024038 (2016).
  * (17) O. B. Zaslavskii, Phys. Rev. D 94, 064048 (2016).
  * (18) E. Berti, R. Brito and V. Cardoso, Phys. Rev. Lett. 114, 251103 (2015).
  * (19) E. Leiderschneider and T. Piran, [arXiv:1501.01984 [gr-qc]].
  * (20) O. B. Zaslavskii, EPL 111, 50004 (2015).
  * (21) K. Ogasawara, T. Harada, U. Miyamoto, T. Igata, and M. Patil, Phys. Rev. D 95, 124019 (2017).
  * (22) K. I. Maeda, K. Okabayashi, and H. Okawa, Phys. Rev. D 98, 064027 (2018).
  * (23) K. Okabayashi and K. i. Maeda, PTEP 2020, 013E01 (2020).
  * (24) N. Tsukamoto and C. Bambi, Phys. Rev. D 91, 104040 (2015).
  * (25) O. B. Zaslavskii, Phys. Rev. D 98, 104030 (2018).
  * (26) M. Patil, T. Harada, K. i. Nakao, P. S. Joshi, and M. Kimura, Phys. Rev. D 93, 104015 (2016).
  * (27) O. B. Zaslavskii, JETP Lett. 92, 571-574 (2010).
  * (28) O. B. Zaslavskii, Phys. Rev. D 86 124039 (2012).
  * (29) H. Nemoto, U. Miyamoto, T. Harada, and T. Kokubu, Phys. Rev. D 87 127502 (2013).
  * (30) M. Bhat, S. Dhurandhar, and N. Dadhich, J. Astrophys. Astr. 6, 85-100 (1985).
  * (31) K. i. Nakao, H. Okawa and K. i. Maeda, PTEP 2018 013E01 (2018).
  * (32) T. Kokubu, Phys. Rev. D 102, 064032 (2020).
  * (33) T. Kokubu, S. L. Li, P. Wu, and H. Yu, Phys. Rev. D 104, 104047 (2021).
  * (34) N. Tsukamoto, M. Kimura, and T. Harada, Phys. Rev. D 89, 024020 (2014).
  * (35) N. Shiraishi, K. Saito, and H. Tasaki, Phys. Rev. Lett. 117, 190601 (2016).
  * (36) C. Armaza, M. Bañados, and B. Koch, Class. Quant. Grav. 33, 105014 (2016).
  * (37) O. B. Zaslavskii, Phys. Rev. D 88, 044030 (2013).
  * (38) N. Tsukamoto and T. Kokubu, Phys. Rev. D 101, 044030 (2020).


Choice of training label matters: how to best use deep learning for
quantitative MRI parameter estimation

Sean C. Epstein1*, Timothy J. P. Bray2, Margaret Hall-Craggs2, Hui Zhang1

1 Department of Computer Science & Centre for Medical Image Computing,
University College London, London, United Kingdom

2 Centre for Medical Imaging, University College London, London, United
Kingdom

<EMAIL_ADDRESS>
## Abstract

Deep learning (DL) is gaining popularity as a parameter estimation method for
quantitative MRI. A range of competing implementations have been proposed,
relying on either supervised or self-supervised learning. Self-supervised
approaches, sometimes referred to as unsupervised, have been loosely based on
auto-encoders, whereas supervised methods have, to date, been trained on
groundtruth labels. These two learning paradigms have been shown to have
distinct strengths. Notably, self-supervised approaches offer lower-bias
parameter estimates than their supervised alternatives. This result is
counterintuitive – incorporating prior knowledge with supervised labels
should, in theory, lead to improved accuracy. In this work, we show that this
apparent limitation of supervised approaches stems from the naïve choice of
groundtruth training labels. By training on labels which are deliberately not
groundtruth, we show that the low-bias parameter estimation previously
associated with self-supervised methods can be replicated – and improved on –
within a supervised learning framework. This approach sets the stage for a
single, unifying, deep learning parameter estimation framework, based on
supervised learning, where trade-offs between bias and variance are made by
careful adjustment of training label.

##  1 Introduction

Magnetic resonance imaging (MRI) is widely regarded as the premier clinical
imaging modality, in large part due to the unparalleled range of contrast
mechanisms available to it. Conventional MRI exploits this contrast in a
purely qualitative manner: images provide only relative information, such that
voxel intensities are only meaningful in the context of their neighbours. In
contrast, quantitative MRI (qMRI) provides quantitative images, where voxel
intensities are directly, and meaningfully, related to underlying tissue
properties. Compared to conventional MRI, this approach promises increased
reproducibility, interpretability, and tissue insight, at the cost of time-
intensive image acquisition and post-processing [1].

One of the biggest time and resource bottlenecks in post-processing is
parameter estimation, whereby a signal model is fit to the intensity variation
across multiple MR images acquired at different experimental settings. Each
voxel requires its own independent model fit: solving for the signal model
parameters that best described the single voxel’s data. The computational cost
of this curve-fitting process, which scales with both voxel number and model
complexity, has become a bottleneck for modern qMRI experiments.

Accelerating curve fittings with deep learning (DL) was first proposed more
than 30 years ago[2], but has only recently gained popularity within the qMRI
community[3, 4, 5, 6, 7]. Just like traditional methods, DL relies on model
fitting, but the model being fit is a fundamentally different one. Instead of
fitting a qMRI signal model to a single voxel of interest (i.e. curve
fitting), DL methods fit (“train”) a deep neural network (DNN) model to an
ensemble of training voxels. This model maps a single voxel’s signal to its
corresponding qMRI parameters; the unknowns in its fitting are network
weights, rather than qMRI parameters. Once this DNN model has been fit to
(“trained on”) the training data, parameter estimation is reduced to simply
applying it to new unseen data, one voxel at a time. This approach offers two
broad advantages over traditional fitting: (1) computational cost is
amortised: despite being more computationally expensive than one-voxel signal
model fitting, DL training only needs to be performed once, for any number of
voxels; once trained, networks provide near-instantaneous parameter estimates
on new data, and (2) computational cost is front-loaded: model training can be
performed away from the clinic, before patient data is acquired.

To date, most DL qMRI fitting methods have been implemented within a
supervised learning framework [3, 4, 8, 5, 7, 9, 10, 11]. This approach trains
DNNs to predict groundtruth qMRI model parameters from noisy qMRI signals.
When compared to conventional fitting, this approach has been found to produce
high bias, low variance parameter estimates [12, 11].

An alternative class of DL methods has also been proposed, sometimes referred
to as unsupervised learning[6, 13], but more accurately described as self-
supervised [14]. In this framework, training labels are not explicitly
provided, but are instead extracted by the network from its training input.
This label generation is designed such that the network learns to predict
signal model parameters corresponding to noise-free signals that most-closely
approximate noisy inputs. This self-supervised approach has been found to
produce similar results to conventional non-DL fitting, i.e. lower bias and
higher variance than its groundtruth-labelled supervised alternative [6, 12].

From an information theoretic standpoint, the comparison between supervised
and self-supervised performance raises an obvious unanswered question. How can
it be that supervised methods, which provide strictly more information during
training than their self-supervised counterparts, produce more biased
parameter estimates?

In this work we answer this question by showing that this apparent limitation
of supervised approaches stems purely from the selection of groundtruth
training labels. By training on labels which are deliberately not groundtruth,
we show that the low-bias parameter estimation previously associated with
self-supervised methods can be replicated – and improved on – within a
supervised learning framework.

This approach sets the stage for a single, unifying, deep learning parameter
estimation framework, based on supervised learning, where trade-offs between
bias and variance can be made, on an application-specific basis, by careful
adjustment of training label.

The rest of the paper is organized as follows: Section 2 describes existing DL
parameter estimation approaches, our proposed method, and how they relate to
each other; Section 3 describes the evaluation of our method and its
comparison to the state of the art; Section 4 contains our findings; and
Section 5 summarizes the contribution and discusses future work.

##  2 Theory

Quantitative MRI extracts biomarkers $y$ from MR data $x$, producing
quantitative spatial maps. We here describe existing voxelwise approaches to
this problem (conventional fitting and DL alternatives) as well as our
proposed novel method.

###  2.1 Conventional iterative fitting

This method, which relies on maximum likelihood estimation (MLE), extracts
biomarkers by performing a voxelwise model fit every time new data is
acquired. An appropriate signal model $M$ is required, parameterised by
$n_{y}$ parameters of interest; for each combination of $y$, the probability
of observing the acquired data $x$ is known as the likelihood $L$ of those
parameters:

$L(x,z|y,\epsilon)=\prod_{i=1}^{n_{z}}P(x_{i},z_{i}|y,\epsilon)$ (1)

for $n_{z}$ acquisitions from sampling scheme $z$ and noise model $\epsilon$.
The model parameters $\hat{y}$ which maximise the likelihood $L$ are assumed
to best represent the tissue contained within the voxel of interest:

$\hat{y}=\operatorname*{arg\,max}_{y}L(x,z|y,\epsilon)$ (2)

Under a Gaussian noise model, this likelihood maximization reduces to the
commonly-used non-linear least squares (NLLS):

$\hat{y}=\operatorname*{arg\,min}_{y}\sum_{i=1}^{n_{z}}\lVert
M(z_{i}|y)-x_{i}\rVert^{2}$ (3)

under the assumption of signal model $M$ associated with groundtruth
biomarkers $y_{gt}$, sampling scheme $z$, and noise $\epsilon$:

$x=M(z|y_{gt})+\epsilon$ (4)

Each of these optimisations has $n_{y}$ unknowns, which are solved for
independently across different voxels; the computational cost scales linearly
with the number of voxels $n_{v}$.

Developments in qMRI acquisition and analysis have led to increased (a) image
spatial resolution (i.e. greater $n_{v}$) and (b) model complexity (i.e.
greater $n_{y}$), such that conventional MLE fitting has become increasingly
computationally expensive.

###  2.2 Existing deep learning methods

Deep learning approaches address this by reframing $n_{v}$ independent
problems into a single global model fit: learning the function $\mathcal{F}$
that maps any $x$ to its corresponding $y_{gt}$:

$y_{gt}=\mathcal{F}(x)$ (5)

Deep neural networks aim to approximate this function by composing a large but
finite number of building-block functions, parametrised by $n_{p}$ network
parameters $p$ (“weights”):

$\hat{y}=\hat{\mathcal{F}}(x|\hat{p})$ (6)

In this context, model fitting (“training”), is performed over network
parameters $p$ and involves maximising $\hat{\mathcal{F}}$’s mean performance
over a large set of training examples; the trained network is defined by the
best-fit parameters $\hat{p}$. This fitting problem, whilst more
computationally expensive to solve than any individual voxel ($n_{v}=1$) MLE,
is only tackled once; once $\hat{\mathcal{F}}$ is learnt, it can be applied at
negligible cost to new, unseen, data. This promise of rapid, zero-cost
parameter estimation has led to the development of two broad classes of DL-
based parameter estimation methods.

SupervisedGT methods approximate $\mathcal{F}$ by minimising the difference
between a large number of noise-free training labels (groundtruth parameter
values) and corresponding network outputs (noise-free parameter estimates);
training loss is calculated in the parameter space $Y$:

$\text{{Supervised\textsubscript{GT} training
loss}}=\sum_{i=1}^{n_{train}}\lVert W\cdot(\hat{y}_{i}-y_{gt,i})\rVert^{2}$
(7)

where $n_{train}$ is the number of training samples and $W$ is a tunable
weight matrix which accounts for magnitude differences in signal model
parameters. $W$ is generally a diagonal matrix, with each diagonal element
$W_{ii}$ corresponding to the relative weighting of qMRI parameter $y_{i}$;
setting $W$ as the identity matrix equally weights all parameters in the
training loss.

These methods produce higher bias, lower variance parameter estimation than
conventional MLE fitting [12, 11] and, by adjusting $W$, can be tailored to
selectively boost estimation performance on a subset of the parameter space
$Y$.

In contrast, Self-supervised methods compute training loss within the signal
space $X$, by minimising the difference between network inputs (noisy signals)
and a filtered representation of network outputs (noise-free signal
estimates):

$\text{{Self-supervised training loss}}=\sum_{i=1}^{n_{train}}\lVert
M(z|\hat{y}_{i})-x_{i}\rVert^{2}$ (8)

These methods, which perform similarly to conventional MLE fitting, produce
lower bias, higher variance parameter estimation than SupervisedGT[12, 6].
Unlike SupervisedGT, the relative loss weighting of different signal model
parameters is dictated by sampling scheme $z$.

Under Gaussian noise conditions, single-voxel Self-supervised loss (i.e.
minimising the sum of squared differences between a noisy signal and its
noise-free signal estimate) is indistinguishable from the corresponding
objective function in conventional fitting.

In contrast, under the Rician noise conditions encountered in MRI acquisition,
Self-supervised training loss no longer matches conventional fitting. Indeed,
the sum of squared errors between noisy signals and noise-free estimates is
not an accurate difference metric in the presence of Rician noise.

To summarise: existing supervised DL techniques are associated by high
estimation bias, low variance, and end-user flexibility; in contrast, self-
supervised methods have lower bias, higher variance, but are limited by the
fact their loss is calculated in the signal space $X$.

###  2.3 Proposed deep learning method

In light of this, we propose SupervisedMLE, a novel parameter estimation
method which combines the advantages of SupervisedGT and Self-supervised
methods. This method is contrasted to existing techniques in Fig 1.

This method mimics Self-supervised’s low-bias performance by learning a
regularised form of conventional MLE, but does so in the parameter space $Y$,
within a supervised learning framework. This addresses the limitations of
Self-supervised: Rician noise modelling is incorporated, and parameter loss
weighting is not limited by sampling scheme $z$.

Our method learns $\hat{\mathcal{F}}$ by training on noisy signals paired with
pre-computed conventional MLE labels. These labels act as proxies for the
groundtruth parameters we wish to estimate:

$\text{{Supervised\textsubscript{MLE} training loss}
}=\sum_{i=1}^{n_{train}}\lVert W\cdot(\hat{y}_{i}-y_{MLE,i})\rVert^{2}$ (9)

Our method offers one final advantage over Self-supervised approaches. In
addition to the parameter estimation improvements relating to noise model
correction and parameter loss weighting, it naturally interfaces with
SupervisedGT. In so doing, it presents the opportunity to combine low-bias and
low-variance methods into a single, tunable hybrid approach, by a simple
weighted sum of each method’s loss function:

$\text{{Hybrid training loss}
}=\alpha\cdot\text{{Supervised\textsubscript{MLE} loss}
}+(1-\alpha)\cdot\text{{Supervised\textsubscript{GT} loss}}$ (10)

Fig 1: Comparison between our proposed method (SupervisedMLE) and existing
supervised and self-supervised approaches.

##  3 Experimental evaluation

Three classes of network were investigated and compared: SupervisedGT, Self-
supervised, and SupervisedMLE, as described in Fig 1. Additionally, to control
for differences in loss function weighting between supervised and unsupervised
methods, Self-supervised was converted into supervised form by training
SupervisedMLE on Gaussian-model based MLE labels. All models are summarised in
Table 1.

Network name |  Loss space |  Label |  Label noise model
---|---|---|---
SupervisedGT |  $Y$ |  Groundtruth |  N/A
Self-supervised |  $X$ |  N/A |  N/A
SupervisedMLE, Rician |  $Y$ |  MLE |  Rician
SupervisedMLE, Gaussian |  $Y$ |  MLE |  Gaussian
Table 1: Summary of evaluated parameter estimation networks. $Y$ denotes
parameter space; $X$ denotes signal space.

All networks were trained and tested on the same datasets; differences in
performance can be attributed solely to differences in loss function
formulation and training label selection.

###  3.1 Signal model

The intravoxel incoherent motion (IVIM) model [15] was investigated as an
exemplar 4-parameter non-linear qMRI model which poses a non-trivial model
fitting problem and is well-represented in the DL qMRI literature [4, 6, 13,
16, 17]:

$S(b|S_{0},f,D_{slow},D_{fast})=S_{0}(fe^{-b(D_{fast}+D_{slow})}+(1-f)e^{-bD_{slow}})$
(11)

where $S$ corresponds to the signal model $M$, $b$ corresponds to the sampling
scheme $z$, and $[S_{0},f,D_{slow},D_{fast}]$ corresponds to the parameter-
vector $y$. In physical terms, $S_{0}$ is an intensity normalisation factor,
$f$ is a perfusion fraction, $D_{slow}$ is the diffusivity of non-perfusing
water, and $D_{fast}$ is the pseudo-diffusivity of perfusing water.

###  3.2 Network architecture

Network architecture was harmonised across all network variants, and
represents a common choice in the existing qMRI literature [6]: 3 fully
connected hidden layers, each with a number of nodes matching the number of
signal samples $z$ (i.e. b-values), and an output layer with a number of nodes
matching the number of model parameters. Wider (150 nodes per layer) and
deeper (10 hidden layers) networks were investigated and found to have
equivalent performance, during both training and testing, at the cost of
increased training time. All networks were implemented in Pytorch 1.9.0 with
exponential linear unit activation (ELU) functions [18]; ELU performance is
similar to ReLU, but is more robust to poor network weight initialisation.

###  3.3 Training data

Training datasets were generated at SNR $=[15,30]$ to investigate parameter
estimation performance at both high and low noise levels. At each SNR, 100,000
noise-free signals were generated from uniform IVIM parameter distributions
($S_{0}\in[0.8,1.2]$, $f\in[0.1,0.5]$, $D_{slow}\in[0.4,3.0]10^{-3}mm^{2}/s$,
$D_{fast}\in[10,150]10^{-3}mm^{2}/s$, representing realistic tissue values),
sampling them with a real-world acquisition protocol [19]
($b=[0,10,20,30,50,80,100,200,400,800]$ $s/mm^{2}$), and adding Rician noise.
Training data generative parameters were drawn from uniform, rather than in-
vivo, parameter distributions to minimise bias in network parameter
estimation[11]. Data were split 80/20 between training and validation. MLE
labels were calculated using a bound-constrained non-linear fitting algorithm,
implemented with scipy.optimize.minimize, using either Rician log-likelihood
or sum of squared errors as fitting objective function. This algorithm was
initialised with groundtruth values to improve fitting robustness and avoid
local minima. Training/validation samples associated with ’poor’ MLE labels
(defined as lying on the boundary of the bound-constrained estimation space)
were held out during training and ignored during validation.

###  3.4 Network training

Network training was performed using an Adam optimizer (learning rate = 0.001,
betas = (0.9, 0.999), weight decay=0) as follows: SupervisedGT (at SNR 30) was
trained 16 times on the same data, each time initialising with different
network weights, to improve robustness to local minima during training. From
this set of trained networks, a single SupervisedGT network was selected on
the basis of validation loss. The trained weights of this selected network
were subsequently used to initialise all other networks; in this way, any
differences in network performance could be solely attributed to differences
in training label selection and training loss formulation. In the case of
supervised loss formulations, the inter-parameter weight vector $W$ was chosen
as the inverse of each parameter’s mean value over the training set, to obtain
equal loss weighting across all four IVIM parameters.

###  3.5 Testing data

Networks were tested on both synthetic and real qMRI data. The synthetic
approach offers (a) known parameter groundtruths to assess estimation against,
(b) arbitrarily large datasets, and (c) tunable data distributions, but is
based on possibly simplified qMRI signals. This approach was used to assess
parameter estimation performance in a controlled, rigorous manner; real data
was subsequently used to validate the trends observed in silico.

Synthetic data was generated with sampling, parameter distributions, and noise
levels matching those used in network training. The IVIM parameter space in
which the networks were trained was uniformly sub-divided 10 times in each
dimension, to analyse estimation performance as a function of parameter value.
At each point in the parameter space, 500 corresponding noisy signals were
generated and used to test network performance, accounting for variation under
noise repetition.

Real data was acquired from the pelvis of a healthy volunteer on a wide-bore
3.0T clinical system (Ingenia, Philips, Amsterdam, Netherlands), 5 slices, 224
x 224 matrix, voxel size = 1.56 x 1.56 x 7mm, TE = 76ms, TR = 516ms, scan time
= 44s per 10 b-values listed in § 3.3. For the purposes of assessing parameter
estimation methods, we obtained gold standard voxelwise IVIM parameter
estimates from a supersampled dataset (16-fold repetition of the above
acquisition, within a single scanning session, generating 160 b-values, total
scan time = 11m44s). Conventional MLE was performed on this supersampled data
to produce best-guess “groundtruth” parameters. During testing, the
supersampled dataset was split into 16 distinct 10 b-value acquisitions, each
corresponding to a single realistic clinical acquisition. All images were
visually confirmed to be free from motion artefacts. The mismatch in parameter
distributions between this in-vivo data (highly non-uniform) and the
previously-described synthetic data (uniform by construction) limited the
scope for validating our in-silico results. To address this, a final synthetic
testing dataset was generated from the in-vivo MLE-derived “groundtruth”
parameters, and was used for direct comparison between real and simulated
data.

###  3.6 Evaluation metrics

Parameter estimation performance was evaluated using 3 key metrics: (1) mean
bias with respect to groundtruth, (2) mean standard deviation under noise
repetition, and (3) root mean squared error (RMSE) with respect to
groundtruth. RMSE is the most commonly used metric to evaluate estimation
performance [6, 4], but is limited in its ability to disentangle accuracy and
precision; to this end, mean bias and standard deviation were used as more
specific measures of network performance.

It is important to note that all methods were assessed with respect to
groundtruth qMRI parameters, even those trained on MLE labels. For these
methods, the training and validation loss (MLE-based) differed from the
reported testing loss (groundtruth-based).

##  4 Results & discussion

This section summarises our main findings and discusses the advantages offered
by the parameter estimation method we propose.

###  4.1 Comparison of parameter estimation methods

The relative performance of all previously-discussed parameter estimation
methods is summarised in Figures 2 and 3. These figures show the bias,
variance (represented by its square root: standard deviation), and RMSE of
parameter estimates with respect to groundtruth values, reported for each
model parameter as a function of its value over the synthetic test dataset;
each plotted point represents an average over 500 noise instantiations and a
marginalisation over all non-visualised parameters. Marginalisation was
required for visualisation of a 4-dimensional parameter space, but was
confirmed to be representative of the entire, non-marginalised space, as
discussed in §4.5.

Fig 2: Parameter estimation performance at low SNR (15) as a function of
groundtruth parameter $Y$. Performance summarised by bias & RMSE with respect
to groundtruth and standard deviation with respect to noise repetition.
Conventional MLE fitting is provided as a non-DNN reference standard. For the
sake of visualisation, each plotted point represents marginalisation over all
non-specified ${Y}$ dimensions. Fig 3: Parameter estimation performance,
visualised as in Figure 2, but for high SNR (30) data.

In keeping with previously reported results, we show a bias/variance trade-off
between different parameter estimation methods. Conventional MLE fitting is
provided as a reference (plotted in black). Approaches which, on a theoretical
level, approximate conventional MLE (Self-supervised and SupervisedMLE,
plotted in red), are generally associated with low bias, high variance, and
high RMSE, whereas groundtruth-labelled supervised methods (plotted in blue)
exhibit lower variance and RMSE at the cost of increased bias.

Increases in bias, if consistent across parameter space $Y$, do not
necessarily reduce sensitivity to differences in underlying tissue properties.
However, we show that $Supervised_{GT}$ is associated with bias that varies
significantly as a function of groundtruth parameter values. This results in a
reduction in information content, visualised as the gradient of the bias plots
(top row) in Fig 2. The more negative the gradient, the more parameter
estimates are concentrated in the centre of the parameter estimation space
$\hat{X}$, and the lower the ability of the method to distinguish differences
in tissue properties. This information loss can be seen in Fig 4, which
compares $Supervised_{GT}$ to conventional MLE fitting, and shows the
compression in $\hat{X}$ over the groundtruth parameter-space $X$.

Fig 4: Comparison between SupervisedGT and reference conventional MLE fitting,
expressed in terms of estimation bias and information compression at low SNR
(15). Arrows represent the mean mapping from $Y$ to $\hat{Y}$, averaged over
noise, as a function of parameter space $Y$. For the sake of visualisation,
each plotted point represents marginalisation over all non-specified ${Y}$
dimensions.

###  4.2 Validation against clinical data

The above trends, found in simulation, were also observed in real-world data.
Fig 5 shows the bias, variance, and RMSE of parameter estimates with respect
to “groundtruth” values (obtained from the supersampled dataset described in
§3.5). The $x$ axes of these plots correspond to these reference values. To
aid visualisation, 10 uniform bins were constructed along each parameter
dimension, into which clinical voxels were assigned based on their
“groundtruth” parameter values. Fig 5 plots the mean bias, standard deviation,
and RMSE associated with each bin as a function of the bin’s central value,
together with the distribution of voxels across the 10 bins.

By calculating the variance of the 16 $b=0$ images, the SNR of this clinical
dataset was found to be $\sim$15; Fig 2 is therefore the relevant point of
comparison. It can be readily seen that the trends observed in simulated data,
described in §4.1, are replicated for $f<0.40$, $D_{slow}<1.5$, and the entire
range of $D_{fast}$, namely the regions of parameter-space which are well-
represented in the real-world data. Fig 6 confirms that divergence outside of
these ranges is due to under-representation in the in vivo test data; the
apparent divergences can be replicated in-silico by matching real-world
parameter distributions.

Fig 7 contains exemplar parameter maps from the clinical test data, and shows
the real-world implications of the trends summarised in Figures 2 and 5:
$Supervised_{GT}$’s low-variance, low-RMSE parameter estimation results in
artificially smooth IVIM maps biased towards mean parameter values.

Fig 5: In vivo parameter estimation performance of networks trained on low SNR
(15) synthetic data, as a function of supersampling-derived reference
parameter values. The first three rows summarise performance by showing bias &
RMSE with respect to reference value and standard deviation with respect to
noise repetition, marginalised over all non-specified ${Y}$ dimensions. The
bottom row shows the distribution of reference parameter values across the
parameter range being visualised. Fig 6: Parameter estimation performance of
networks trained on low SNR (15) synthetic data, tested on a synthetic dataset
matching the distribution of in vivo reference parameter values. The first
three rows summarise performance by showing bias & RMSE with respect to
groundtruth value and standard deviation with respect to noise repetition,
marginalised over all non-specified ${Y}$ dimensions. The bottom row shows the
distribution of groundtruth parameter values across the parameter range, which
matches the in vivo dataset by construction. Fig 7: Parameter estimation
performance of networks on real-world test data, visualised as spatial maps.
Groundtruth maps are taken as the maximum likelihood parameter estimates
associated with the complete 160 b-value dataset, whereas network predictions
are obtained from a single 10 b-value subsample.

###  4.3 Advantages offered by our method

Our proposed method occupies the low-bias side of the bias-variance trade-off
discussed in §3.5, and offers four broad advantages over the competing method
in this space (Self-supervised): (1) flexibility in choosing inter-parameter
loss weighting $W$, (2) incorporation of non-Gaussian (e.g. Rician) noise
models, (3) compatibility with complex, non-differentiable signal models $M$,
and (4) ability to interface with low-variance methods, to produce a hybrid
approach tunable to the needs of the task at hand. These advantages are
analysed in turn.

####  4.3.1 Choice of inter-parameter loss weighting $W$

By computing loss in parameter-space $Y$, our method has total flexibility in
adjusting the relative contribution of different $y$ in the training loss
function. In contrast, since Self-supervised calculates training loss in $X$,
the relative weighting depends on the acquisition protocol $z$. Fig 8 compares
our method - weighted so as to not discriminate between different model
parameters - with variants designed to overweight single parameters by a
factor of $10^{6}$. The potential advantages offered by this selective
weighting are seen in the estimation $D_{fast}$, where this approach leads to
a small increase in both precision and accuracy. This parameter-specific
weighting is not accessible within a Self-supervised framework.

In light of the differences arising from inter-parameter loss weighting, for
subsequent analysis we use SupervisedMLE, Gaussian as a proxy for Self-
supervised; both methods encode the same regularised MLE fitting, but differ
in their inter-parameter weighting.

Fig 8: Comparison between SupervisedMLE, Rician, as described above, and
variants which differ in their inter-parameter loss weighting $W$, at high SNR
(30). Each column compares SupervisedMLE, Rician to a different network
variant, uniquely trained to overweight the single relevant signal model
parameter. For the sake of visualisation, each plotted point represents
marginalisation over all non-specified ${Y}$ dimensions.

####  4.3.2 Incorporation of Rician noise modelling

By pre-computing MLE labels using conventional parameter estimation methods,
we are able to incorporate accurate Rician noise modelling. Comparison between
SupervisedMLE, Rician and SupervisedMLE, Gaussian shows the effect of the
choice of noise model; these differences are most pronounced at low SNR (Fig
2) and high $D_{slow}$, when the Gaussian approximation of Rician noise is
known to break down. In this regime, our method gives less biased, more
informative $D_{slow}$ estimates, replicating conventional MLE performance at
a fraction of the computational cost. At high $D_{slow}$, our method has a
flatter, more information-rich, $D_{slow}$ bias curve than all other DL
methods. This information loss is further visualised in Fig 9, which shows the
compression in $D_{slow}$ estimates $\hat{X}$ over the groundtruth parameter-
space $X$. As expected, this compression is most apparent at high values of
$D_{slow}$, when the signal is more likely to approach the Rician noise floor.

Fig 9: Comparison of the information content captured by SupervisedMLE
methods, as a function of the noise model used in computing MLE labels, at low
SNR (15). Arrows represent the mean mapping from $Y$ to $\hat{Y}$, averaged
over noise, as a function of parameter space $Y$. For the sake of
visualisation, each plotted point represents marginalisation over all non-
specified $Y$ dimensions.

####  4.3.3 Compatibility with complex signal models

An additional advantage of computing training loss in parameter-space $Y$ is
that DNN networks are signal model agnostic: network training does not require
explicit calculation of $M$. This approach is advantageous when working with
complex signal models, as made clear by comparison with Self-supervised
methods. In contrast with our proposed approach, Self-supervised methods embed
$M$ between network output and training loss (see Fig 1). During DNN training,
network parameters $p$ are updated by computing partial derivatives of the
training loss; this process requires the loss to be expressed in a
differentiable form. By embedding $M$ in the loss formulation, Self-supervised
methods thus limit themselves to signal models that can be expressed in an
explicitly differentiable form. In contrast, our method does not rely on
computing $M$ during training, and is therefore compatible with a wider range
of complex qMRI signal models.

####  4.3.4 Tunable network approach

As discussed above, we show a clear bias/variance trade-off between different
parameter estimation methods. The optimal choice of method depends on the task
at hand [20], and may not lie at either extreme of this trade-off. Therefore,
it would be advantageous to be able to combine low-bias and low-variance
methods into a single, hybrid approach, with performance tunable by the
relative contribution of each constituent method. Our proposed method, which
interfaces naturally with $Supervised_{GT}$, offers exactly that. An example
of this approach is shown in Fig 10: training loss has been weighted equally
between groundtruth and MLE labels, and, as expected, the resulting network
performance lies in the middle ground between these two extremes.

Fig 10: Proof of concept of a hybrid parameter estimation method, formed by
training a supervised network with an equally-weighted sum of SupervisedMLE,
Rician and SupervisedGT loss functions, at low SNR (15). For the sake of
visualisation, each plotted point represents marginalisation over all non-
specified $Y$ dimensions.

####  4.3.5 Comparison with conventional fitting

Comparison between our proposed method (SupervisedMLE, Rician) and
conventional fitting (MLE, Rician) highlights additional advantages offered by
our approach. Firstly, Figs 2 and 3 demonstrate qualitatively similar
performance between these methods across the entire parameter space. The fact
that our method, which offers near-instantaneous parameter estimation,
produces similar parameter estimates to well-understood conventional MLE
methods justifies its adoption in and of itself. However, our method not only
mimics but indeed in many cases outperforms (lower bias, variance, and RMSE)
the very same method used to compute those labels. This result not only
motivates its use, but also confirms that DL methods are able to exploit
information shared between training samples beyond what would be possible by
considering each sample in isolation.

###  4.4 A note on RMSE

We note that RMSE is a poor summary measure of network performance. RMSE is
heavily skewed by outliers, and thus favours methods which give parameter
estimates consistently close to mean parameter values. Such estimates, as in
the case of $D_{fast}$, may contain very little information (Fig 4) despite
being associated with low RMSE. Accordingly, we strongly recommend that RMSE
be discontinued as a single summary metric for parameter estimation
performance: it must always be accompanied by bias, variance, and ideally an
analysis of information content.

###  4.5 Justification of parameter marginalisation

The above analysis has been largely based on Figs 2 and 3, which show
parameter estimation performance marginalised over 3 dimensions of $X$. This
choice, made to aid visualisation, was validated against higher dimensional
representations of the same data.

Fig 11 compares SupervisedMLE, Rician and SupervisedGroundtruth performance
across the entire qMRI parameter space. It can be seen that trends observed in
Fig 2 are replicated here; we draw attention to two such examples. Firstly,
Fig 2 suggests SupervisedGroundtruth produces lower $f$ standard deviation
than SupervisedMLE, Rician; Fig 11 confirms this to be the case across all
test data. In contrast, Fig 2 suggests that SupervisedGroundtruth produces
higher $D_{slow}$ bias at low $D_{slow}$ and lower bias at high $D_{slow}$;
Fig 11 confirms a spread of bias differences across the test data: some
favouring one method, and others the other. This effect is explored in Fig 12,
which compares $D_{slow}$ estimation performance as a function of $f$ and
$D_{fast}$ at two specific (non-marginalised) groundtruth $D_{slow}$ values
($0.69,2.71)$. As expected from the marginalised representation in Fig 2, at
low $D_{slow}$ SupervisedGroundtruth produces higher bias across the entire
$f$-$D_{fast}$ parameter space, whereas at high $D_{slow}$ the opposite is
true.

Despite this, it is important to note the limitations of marginalisation. Fig
12 also shows that the relative performance of SupervisedMLE, Rician and
SupervisedGroundtruth varies across all parameter-space dimensions. Consider
$D_{slow}=0.69$, where 2 shows similar marginalised RMSE for these methods. In
fact, by visualising this difference as a function of $f$ and $D_{fast}$, we
reveal two distinct regions: high $f$/low $D_{fast}$ (where SupervisedMLE,
Rician produces lower RMSE), and elsewhere (where it produces higher RMSE).
This highlights (1) the potential pitfalls of producing summary results by
marginalising across entire parameter spaces and (2) the need to choose
parameter-estimation methods appropriate for the specific parameter
combinations relevant to the tissues being investigated [20].

Fig 11: Non-marginalised comparison of parameter estimation performance
between SupervisedMLE, Rician and SupervisedGroundtruth at low SNR (15).
Colour intensity represents density of distribution across all $X$ and all
noise repetitions. Fig 12: Differences in performance (bias, standard
deviation, RMSE) between SupervisedMLE, Rician and SupervisedGroundtruth for
two groundtruth values of $D_{slow}$ at low SNR (15). The outermost columns
(left and right) correspond to $D_{slow}=0.69$ and $D_{slow}=2.71$
respectively, and show mean performance under noise repetition, without
marginalisation. The central column reproduces the corresponding marginalised
representation from Fig 2.

###  4.6 Non-voxelwise approaches

This work has focused on voxelwise DL parameter estimation methods: networks
which map one signal curve to its corresponding parameter estimate. There are,
however, alternatives: convolutional neural network (CNN) methods which map
spatially related clusters (“patches”) of qMRI signals to corresponding
clusters of parameter estimates [21, 22, 23]. Our MLE training label approach
could be incorporated into such methods, and we leave it to future work to
investigate the effect this would have on parameter estimation performance.

##  5 Conclusions

In this work we draw inspiration from state-of-the-art supervised and self-
supervised qMRI parameter estimation methods to propose a novel DNN approach
which combines their respective strengths. In keeping with previous work, we
demonstrate the presence of a bias/variance trade-off between existing
methods; supervised training produces low variance under noise, whereas self-
supervised leads to low bias with respect to groundtruth.

The increased bias of supervised DNNs is counter-intuitive - when labels are
available, these methods have access to more information, and should therefore
outperform, non-labelled alternatives. In light of this, we infer that the
high bias associated with these supervised methods stems from the nature of
the additional information they receive: groundtruth training labels. By
careful adjustment of these labels, we show that the low-bias performance
previously limited to self-supervised approaches can be achieved within a
supervised learning framework.

This framework forms the basis of a novel low-bias supervised learning
approach to qMRI parameter estimation: training on conventionally-derived
maximum likelihood parameter estimates. This method offers four clear
advantages to competing non-supervised low-bias DNN approaches: (1)
flexibility in choosing inter-parameter loss weighting, which enables network
performance to be boosted for qMRI parameters of interest; (2) incorporation
of Rician noise modelling, which improves parameter estimation at low SNR; (3)
separation between signal model and training loss, which enables the
estimation of non-differentiable qMRI signal models;and, crucially, (4)
ability to interface with existing supervised low-variance approaches, to
produce a tunable hybrid parameter estimation method.

This final point underpins the key contribution of this work: unifying low-
bias and low-variance parameter estimation under a single supervised learning
umbrella. When faced with a parameter estimation problem, we no longer need to
choose between extremes of the bias/variance trade-off; we can now tune DNN
parameter estimation performance to the specific needs of the task at hand.
This sets the stage for future work, where this tuning constant is optimised
as part of a computational, task-driven, experimental design framework [20].

## References

  * 1. Cercignani M, Dowell NG, Tofts PS.  Quantitative MRI of the Brain.  CRC Press; 2018.
  * 2. Bishop CM, Roach CM.  Fast curve fitting using neural networks.  Review of Scientific Instruments. 1992;63(10):4450–4456.  doi:10.1063/1.1143696.
  * 3. Golkov V, Dosovitskiy A, Sperl JI, Menzel MI, Czisch M, Sämann P, et al.  q-Space Deep Learning: Twelve-Fold Shorter and Model-Free Diffusion MRI Scans.  IEEE Transactions on Medical Imaging. 2016;35(5):1344–1351.  doi:10.1109/TMI.2016.2551324.
  * 4. Bertleff M, Domsch S, Weingärtner S, Zapp J, O’Brien K, Barth M, et al.  Diffusion parameter mapping with the combined intravoxel incoherent motion and kurtosis model using artificial neural networks at 3 T.  NMR in biomedicine. 2017;30(12).  doi:10.1002/NBM.3833.
  * 5. Liu H, Xiang QS, Tam R, Dvorak AV, MacKay AL, Kolind SH, et al.  Myelin water imaging data analysis in less than one minute.  NeuroImage. 2020;210:116551.  doi:10.1016/J.NEUROIMAGE.2020.116551.
  * 6. Barbieri S, Gurney-Champion OJ, Klaassen R, Thoeny HC.  Deep learning how to fit an intravoxel incoherent motion model to diffusion-weighted MRI.  Magnetic Resonance in Medicine. 2020;83(1):312–321.  doi:10.1002/MRM.27910.
  * 7. Palombo M, Ianus A, Guerreri M, Nunes D, Alexander DC, Shemesh N, et al.  SANDI: A compartment-based model for non-invasive apparent soma and neurite imaging by diffusion MRI.  NeuroImage. 2020;215:116835.  doi:10.1016/J.NEUROIMAGE.2020.116835.
  * 8. Yoon J, Gong E, Chatnuntawech I, Bilgic B, Lee J, Jung W, et al.  Quantitative susceptibility mapping using deep neural network: QSMnet.  NeuroImage. 2018;179:199–206.  doi:10.1016/J.NEUROIMAGE.2018.06.030.
  * 9. Aliotta E, Nourzadeh H, Patel SH.  Extracting diffusion tensor fractional anisotropy and mean diffusivity from 3-direction DWI scans using deep learning.  Magnetic Resonance in Medicine. 2021;85(2):845–854.  doi:10.1002/MRM.28470.
  * 10. Yu T, Canales-Rodríguez EJ, Pizzolato M, Piredda GF, Hilbert T, Fischi-Gomez E, et al.  Model-informed machine learning for multi-component T2 relaxometry.  Medical Image Analysis. 2021;69:101940.  doi:10.1016/J.MEDIA.2020.101940.
  * 11. Gyori NG, Palombo M, Clark CA, Zhang H, Alexander DC.  Training data distribution significantly impacts the estimation of tissue microstructure with machine learning.  Magnetic Resonance in Medicine. 2022;87(2):932–947.  doi:10.1002/MRM.29014.
  * 12. Grussu F, Battiston M, Palombo M, Schneider T, Wheeler-Kingshott CAMG, Alexander DC.  Deep Learning Model Fitting for Diffusion-Relaxometry: A Comparative Study.  Mathematics and Visualization. 2021; p. 159–172.  doi:10.1007/978-3-030-73018-5.
  * 13. Kaandorp MPT, Barbieri S, Klaassen R, van Laarhoven HWM, Crezee H, While PT, et al.  Improved unsupervised physics-informed deep learning for intravoxel incoherent motion modeling and evaluation in pancreatic cancer patients.  Magnetic Resonance in Medicine. 2021;86(4):2250–2265.  doi:10.1002/MRM.28852.
  * 14. Murphy KP.  Probabilistic Machine Learning: An introduction.  MIT Press; 2022.  Available from: probml.ai.
  * 15. Le Bihan D, Breton E, Lallemand D, Grenier P, Cabanis E, Laval-Jeantet M.  MR imaging of intravoxel incoherent motions: application to diffusion and perfusion in neurologic disorders.  Radiology. 1986;161(2):401–407.  doi:10.1148/radiology.161.2.3763909.
  * 16. Mastropietro A, Procissi D, Scalco E, Rizzo G, Bertolino N.  A supervised deep neural network approach with standardized targets for enhanced accuracy of IVIM parameter estimation from multi-SNR images.  NMR in Biomedicine. 2022;35(10):e4774.  doi:https://doi.org/10.1002/nbm.4774.
  * 17. Rozowski M, Palumbo J, Bisen J, Bi C, Bouhrara M, Czaja W, et al.  Input layer regularization for magnetic resonance relaxometry biexponential parameter estimation.  Magnetic Resonance in Chemistry. 2022;doi:https://doi.org/10.1002/mrc.5289.
  * 18. Clevert DA, Unterthiner T, Hochreiter S.  Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs).  4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings. 2015;doi:10.48550/arxiv.1511.07289.
  * 19. Zhao Yh, Li Sl, Liu Zy, Chen X, Zhao Xc, Hu Sy, et al.  Detection of Active Sacroiliitis with Ankylosing Spondylitis through Intravoxel Incoherent Motion Diffusion-Weighted MR Imaging.  European Radiology. 2015;25(9):2754–2763.  doi:10.1007/s00330-015-3634-2.
  * 20. Epstein SC, Bray TJP, Hall-Craggs MA, Zhang H.  Task-driven assessment of experimental designs in diffusion MRI: A computational framework.  PLOS ONE. 2021;16.  doi:10.1371/JOURNAL.PONE.0258442.
  * 21. Fang Z, Chen Y, Lin W, Shen D.  Quantification of relaxation times in MR Fingerprinting using deep learning.  Proceedings of the International Society for Magnetic Resonance in Medicine Scientific Meeting and Exhibition International Society for Magnetic Resonance in Medicine Scientific Meeting and Exhibition. 2017;25.
  * 22. Ulas C, Das D, Thrippleton MJ, Valdés Hernández MdC, Armitage PA, Makin SD, et al.  Convolutional Neural Networks for Direct Inference of Pharmacokinetic Parameters: Application to Stroke Dynamic Contrast-Enhanced MRI.  Frontiers in Neurology. 2019;9(JAN):1147.  doi:10.3389/fneur.2018.01147.
  * 23. Li S, Wu J, Ma L, Cai S, Cai C.  A simultaneous multi-slice T2 mapping framework based on overlapping-echo detachment planar imaging and deep learning reconstruction.  Magnetic Resonance in Medicine. 2022;87(5):2239–2253.  doi:https://doi.org/10.1002/mrm.29128.