Integrated thermal and magnetic susceptibility modeling for air-motion artifact correction in proton resonance frequency shift thermometry

Abstract

Purpose

Hyperthermia treatments are successful adjuvants to conventional cancer therapies in which the tumor is sensitized by heating. To monitor and guide the hyperthermia treatment, measuring the tumor and healthy tissue temperature is important. The typical clinical practice heavily relies on intraluminal probe measurements that are uncomfortable for the patient and only provide spatially sparse temperature information. A solution may be offered through recent advances in magnetic resonance thermometry, which allows for three-dimensional internal temperature measurements. However, these measurements are not widely used in the pelvic region due to a low signal-to-noise ratio and presence of image artifacts.

Methods

To advance the clinical integration of magnetic resonance-guided cancer treatments, we consider the problem of removing air-motion-induced image artifacts. Thereto, we propose a new combined thermal and magnetic susceptibility model-based temperature estimation scheme that uses temperature estimates to improve the removal of air-motion-induced image artifacts. The method is experimentally validated using a dedicated phantom that enables the controlled injection of air-motion artifacts and with in vivo thermometry from a clinical hyperthermia treatment.

Results

We showed, using probe measurements in a heated phantom, that our method reduced the mean absolute error (MAE) by 58% compared to the state-of-the-art near a moving air volume. Moreover, with in vivo thermometry our method obtained a MAE reduction between 17% and 95% compared to the state-of-the-art.

Conclusion

We expect that the combined thermal and magnetic susceptibility modeling used in model-based temperature estimation can significantly improve the monitoring in hyperthermia treatments and enable feedback strategies to further improve MR-guided hyperthermia cancer treatments.

Keywords:

• Noninvasive thermometry
• artifact correction
• air motion
• modelling
• temperature estimation

1. Introduction

Hyperthermia is a successful adjuvant to cancer treatments as it improves the effects of chemo- and radiotherapy [1]. During a hyperthermia treatment, the tumor is heated to temperatures ranging between 39 and 44 degrees Celsius for a duration of 60–90 min, which sensitizes the tumor. Noninvasive heating of a tumor can be accomplished by, e.g., Radio Frequency (RF) antennas or Focused Ultrasound (FUS) [2–4]. A schematic illustration of a typical treatment setup for RF-based hyperthermia is shown in Figure 1. As the treatment outcome is correlated to the thermal dose in the tumor [5], temperature monitoring and control strategies can improve the treatment quality. Moreover, real-time monitoring of the patient’s internal temperature is important given a strong, in-homogeneous, and patient-specific thermoregulatory response.

Figure 1. Schematic illustration of a typical RF hyperthermia treatment setup for the pelvic region. Note that the waterbolus aids the transmission of the electromagnetic fields into the patient, while simultaneously cooling the skin.

Monitoring the temperature in clinical practice is typically achieved through intraluminal probes. While these sensors have excellent accuracy and sampling rate, they are spatially sparse and uncomfortable for the patient. To alleviate these issues, Magnetic Resonance Imaging (MRI) can be used to obtain three-dimensional (3D) temperature measurements [6]. However, the so-called Magnetic Resonance Thermometry (MRT) is plagued by a low signal-to-noise ratio (SNR) and low sampling rates [7,8]. To compensate for the low SNR and sample rates, model-based state observers seem a promising solution [9–11]. However, possibly the most disrupting measurement corruption, which is currently not corrected for in a model-based observer, is due to the magnetic susceptibility change resulting from air-motion [12,13]. In fact, a recent literature study showed that obtaining an accuracy of 1 degree Celsius is not possible in most regions of the patient [14]. As a result, the limited reliability of MRT hampers its clinical adoption as a replacement for the intraluminal probes despite the benefits of MRT.

Recently, Wu et al. [13] proposed a method for correcting such air-motion-induced susceptibility artifacts. Loosely speaking, [13] employed methods conventionally used in Quantitative Susceptibility Mapping (QSM), namely, the Projection onto Dipole Fields (PDF) [15] and the Laplacian Boundary Value problem (LBV) [16], to remove the dipole-shaped artifacts from thermometry measurements. While both the PDF and LBV algorithm in the method by Wu et al. showed great potential, applying techniques from QSM to thermometry is not straightforward. For example, removing air-motion-induced artifacts using QSM techniques can also result in the removal of part of the temperature-induced phase shift. This issue was partly addressed in [13] by temporarily removing spatially constant and linear terms from the phase measurement prior to removing the air-motion-induced susceptibility artifacts. However, this leaves higher-order terms in the temperature-induced phase shift at risk of being (partially) removed. In particular, the proposed solution is not suitable for scenarios where the heated region is close to air-motion.

To overcome the downsides, in this paper, we propose an integrated thermal and magnetic susceptibility model-based state observer (TSM) to improve the correction of air-motion-induced susceptibility artifacts in MRT. Note that state observers are commonly used to estimate the true state of a system from noisy measurements [17, Ch. 6]. We validate our method using a heated phantom experiment equipped with a movable air volume, see Figure 2. Moreover, we perform an in vivo validation of our TSM method on a data set from a clinical hyperthermia cancer treatment. For both the phantom and in vivo validation, we compare our results to the current state-of-the-art [13].

Figure 2. Schematic cross-sections of the phantom. The phantom is filled with Perfax, which has analogous properties to muscle tissue. In the left cross-section, the tube containing the movable air volume is clearly visible. The black crosses denote four Bowman thermal probes, which are used for validation.

2. Theory

In this section, we introduce the thermal and the air-motion-induced susceptibility artifact models which we will use in Section 3.

2.1. Thermal model

In this work, we model the thermal dynamics of a patient or phantom using the Pennes’ bioheat equation [18], which is a convection-diffusion equation with a distributed source and internal heat loss given by (1a) $$\begin{array}{c}\mathit{\rho }\mathit{c}\dot{T}(\mathit{r},t)=\nabla ·(k\nabla T(\mathit{r},t))+\rho Q_{\text{SAR}}(\mathit{r},t)\\ +h_{\mathrm{b}}(T(\mathit{r},t)-T_{\text{blood}}(\mathit{r}))\end{array}$$(1a) with convective boundary condition (1b) $$k\nabla T(\mathit{r},t)·\mathit{n}=h(T(\mathit{r},t)-T_{\text{bc}}(\mathit{r}))$$(1b) at the skin. Here, the symbols are defined as: T temperature [K], r position vector [m], $Q_{\text{SAR}}$ specific absorption rate [Wkg^–1], ρ density [kgm^–3], c heat capacity [Jkg^–1K^–1], k thermal conductivity [Wm^–1K^–1], h convective heat transfer [Wm^–2 K^–1], $h_{\mathrm{b}}$ blood heat transfer [WK^–1], n outward facing surface normal [–], $T_{\text{blood}}$ blood temperature [K], and $T_{\text{bc}}$ boundary temperature [K], respectively. Note that the thermal parameters are tissue-specific, possibly even in-homogeneous within a single tissue, and therefore position-dependent. However, the position dependence is omitted in the notation for brevity.

Heating is achieved through the absorption of RF waves emitted by m fixed-frequency antennas surrounding the patient, see, e.g., Figure 1. By controlling the power and phase of each antenna the emitted RF waves constructively and destructively interfere and subsequently deposit a localized heat load [19]. Mathematically, the specific absorption rate is given by (2) $$Q_{\text{SAR}}(\mathit{r},t)=\frac{\sigma }{2\rho }{\Vert \sum _{k=1}^m{\mathit{E}}_k(\mathit{r})p_k(t)\Vert }_2^2.$$(2)

Here, σ [Sm^–1], ${\mathit{E}}_k\in {\mathbb{C}}^{3\times 1}$ [Vm^–1], $p_k\in \mathbb{C},$ and ${\parallel ·\parallel }_2$ denote the electrical conductivity, the complex-valued time-harmonic electric field vector per antenna channel, the complex-valued phasor per antenna channel, and the two-norm, respectively. The time-harmonic electric field generated per antenna for $k=1,\dots ,m$ can be computed using a commercial solver such as Sim4Life [20].

Solving (1) with (2) forward in time is typically performed by spatially discretizing the domain to obtain many coupled ordinary differential equation, which can be solved with well-known numerical integration techniques, see, e.g., [10]. Throughout this paper we define $\Delta T(\mathit{r},t)=T(\mathit{r},t)-T(\mathit{r},t_0)$ as a relative temperature with respect to a baseline at time t₀, e.g., room temperature for the phantom and approximately 37 degrees Celsius for a patient. In this paper, we consider a homogeneous baseline temperature, however, the presented method is not limited to this particular choice. Note that using relative temperatures is intuitive as the thermometry also measures temperature relative to a baseline.

2.2. Air-motion-induced susceptibility artifacts

Hyperthermia cancer treatments are typically monitored using Proton Resonance Frequency Shift (PRFS) thermometry. The PRFS principle exploits the temperature-dependent phase measured by the MRI scanner [7,8]. Besides its temperature sensitivity, the MR signal phase is also dependent on the macroscopic magnetic field [21–23], i.e., (3a) $$\Delta \phi (\mathit{r},t)=\phi (\mathit{r},t)-\phi (\mathit{r},t_0),$$(3a) (3b) $$=\gamma \text{TE}(B_{\text{nuc}}(\mathit{r},t)-B_{\text{nuc}}(\mathit{r},t_0)),$$(3b) (3c) $$\approx \gamma \text{TE}\alpha B_0(T(\mathit{r},t)-T(\mathit{r},t_0))$$(3c) $$\begin{array}{c}-\gamma \text{TE}{B}_0\frac{2}{3}(\chi (\mathit{r},t)-\chi (\mathit{r},t_0))\\ +\gamma \text{TE}(B_{\text{mac}}(\mathit{r},t)-B_{\text{mac}}(\mathit{r},t_0)).\end{array}$$

Here, the symbols are defined as: φ phase of the MR-signal [rad], α PRFS sensitivity [K^–1], γ gyromagnetic ratio [rads^–1K^–1], TE echo time [s], B₀ nominal magnetic field [T], $B_{\text{nuc}}$ magnetic field at the nucleus [T], $B_{\text{mac}}$ macroscopic magnetic field [T], and χ magnetic susceptibility distribution [–], respectively. Typically, using (3), the measured temperature is computed as (4) $$\Delta T_{\text{measured}}(\mathit{r},t):=\frac{\Delta \phi (\mathit{r},t)}{\gamma \text{TE}\alpha B_0}.$$(4)

Note that (4) is based on the assumption that both the macroscopic magnetic field and the susceptibility distribution in (3) do not change over time. However, due to, e.g., magnetic field drift, we need to account for changes in the macroscopic field, denoted by $\Delta B_{\text{mac}}(\mathit{r},t)=B_{\text{mac}}(\mathit{r},t)-B_{\text{mac}}(\mathit{r},t_0).$ Additionally, note that the phase shift also depends on the local change in magnetic susceptibility. In this paper, we will neglect the smaller direct influence of the magnetic susceptibility on the phase difference, as we are interested in the more dominant far reaching disturbances that occur with air-motion. To this end, we model the effect of air-motion on the macroscopic magnetic field similar to [22,23] as (5a) $$\Delta B_{\text{mac}}(\mathit{r},t)\approx B_0(\Delta \chi (\mathit{r},t)*D(\mathit{r}))+\Delta B_{\text{drift}}(\mathit{r},t),$$(5a) (5b) $$D(\mathit{r}):=\frac{1}{4\pi }\frac{3{({\mathit{e}}_z^{\top }\mathit{r})}^2{\Vert \mathit{r}\Vert }_2^2}{{\Vert \mathit{r}\Vert }_2^5}.$$(5b)

Here, $\Delta B_{\text{drift}}$ [T], $\Delta \chi$ [–], D [–], ${\mathit{e}}_z$ [–], and $*$ denote the magnetic field drift, the difference in susceptibility distribution, the dipole kernel, a basis vector pointing in the direction of the magnetic field, and the convolution operator, respectively.

2.3. Background field removal

In this paper, we will focus on a particular method of QSM to remove the air-motion-induced susceptibility artifacts, namely the PDF method [15]. In this section, we will briefly summarize the underlying principle and assumptions of the PDF method. For a more detailed description, we refer the interested reader to [15] and the references therein.

A key assumption in the PDF method is that the measured temperature field (or equivalently, phase field) has the following decomposition: (6) $$\begin{array}{l}\Delta T_{\text{measured}}(\mathit{r},t_k)=\Delta T(\mathit{r},t_k)+\Delta T_{B_0}(\mathit{r},t_k)+\Delta T_{\chi }(\mathit{r},t_k).\end{array}$$(6)

Here, $\Delta T(\mathit{r},t_k)$ [K] denotes the relative temperature difference between the baseline at time t₀ and current time t_k, $\Delta T_{B_0}=\frac{\Delta B_{\text{drift}}}{\alpha B_0}$ denotes the apparent temperature resulting from B₀-drift, and $\Delta T_{\chi }=\frac{\Delta \chi *D}{\alpha }$ denotes the dipole-shaped apparent temperature resulting from air-motion, which we will refer to as susceptibility artifacts. Note that methods to estimate $\Delta T_{B_0}$ are well-known [7,8], therefore, we assume that an estimate of $\Delta T_{B_0}$ is available. Next, the PDF method uses an optimization-based strategy to place susceptibility sources $\Delta {\chi }_{\text{PDF}}$ in a pre-defined volume to approximate $\Delta T_{\chi }.$ Given an approximation of $\Delta T_{\chi }\approx \frac{\Delta {\chi }_{\text{PDF}}*D}{\alpha },$ we can correct the thermometry using (7) $$\begin{array}{ll}\Delta T_{\text{corrected}}(\mathit{r},t_k)=& \Delta T_{\text{measured}}(\mathit{r},t_k)-\Delta T_{B_0}(\mathit{r},t_k)\\ & -{\alpha }^{-1}\Delta {\chi }_{\text{PDF}}(\mathit{r},t_k)*D(\mathit{r}).\end{array}$$(7)

Computing the susceptibility sources $\Delta {\chi }_{\text{PDF}}$ is performed by solving the optimization problem (8) $$\begin{array}{ll}\Delta {\chi }_{\text{PDF}}(\mathit{r},t_k):=& \underset{\Delta \chi \in {\Omega }_{\mathrm{S}}}{\arg \hspace{0.17em}\min }{\int }_{{\Omega }_{\text{ROI}}}(\Delta T_{\text{measured}}(\mathit{r},t_k)\\ & -\Delta T_{B_0}(\mathit{r},t_k)-{\alpha }^{-1}\Delta \chi (\mathit{r})*D(\mathit{r}){)}^2\mathrm{d}V.\end{array}$$(8)

Here, $\Delta {\chi }_{\text{PDF}}$ [–], ${\Omega }_{\mathrm{S}},$ ${\Omega }_{\text{ROI}}$ denote the identified susceptibility difference distribution, a mask selecting regions with susceptibility artifact sources, and a mask selecting a region of interest in which the clinician would like to measure the temperature, i.e. the patient. The optimization problem (8) can be efficiently solved with the MEDI toolbox [24].

2.4. Mask selection

Approximating the background field containing the air-motion susceptibility artifacts using (8) should be done carefully. As mentioned in [15], the PDF method can partly remove the temperature-induced phase shift, i.e. $\Delta T(\mathit{r},t_k)$ in (6). Loosely speaking, when ${\Omega }_{\mathrm{S}}$ is chosen too large and includes more than just the susceptibility artifact sources, more of the temperature-induced phase shift is removed. With the same logic, choosing ${\Omega }_{\mathrm{S}}$ too small limits the air-motion artifact removal. Clearly, choosing ${\Omega }_{\mathrm{S}},$ i.e. the moving air volume, is a crucial aspect in removing air-motion artifacts. To this end, we construct ${\Omega }_{\mathrm{S}}$ by comparing the anatomical scan at the baseline t₀, and at the current time t_k. Any sufficiently large change in intensity between anatomical scans is assumed to indicate air volume movement, and are thus included in ${\Omega }_{\mathrm{S}}.$ Thereto, we compute ${\Omega }_{\mathrm{S}}$ by applying a threshold to the absolute difference between anatomical scans followed by a morphological closing and opening to close holes and remove small objects. More formally, ${\Omega }_{\mathrm{S}}$ is defined as (9) $${\Omega }_{\mathrm{S}}:=\left(({\Omega }_{\mathrm{S}}^0\oplus \mathcal{S}(d_1))\ominus \mathcal{S}(d_1)\ominus \mathcal{S}(d_2)\right)\oplus \mathcal{S}(d_2),$$(9) where ${\Omega }_{\mathrm{S}}^0:=\left\{\mathit{r}\right||s(\mathit{r},t_0)-s(\mathit{r},t_k)|\ge \text{threshold}\}$ denotes the unprocessed mask, s denotes the signal magnitude corresponding the anatomic MR scan, $\mathcal{S}(d):=\{\mathit{r}{\parallel \mathit{r}\parallel }_2\le d\}$ denotes a structuring element used to for morphological opening and closing operations, d₁ and d₂ denote the radius of the spherical structuring element, and ⊕, ⊖ denote the Minkowski sum and difference, respectively. Note that (9) is easily implemented with binary dilation and erosion algorithms found in image processing toolboxes.

The concept of using the difference between two anatomic scans is illustrated in Figure 3 using a phantom equipped with a movable air volume.

Figure 3. Illustrating the susceptibility source mask selection ${\Omega }_{\mathrm{S}}$ by detecting air movement using anatomic MR scans. In subfigures A and B, the anatomic baseline scan $s(\mathit{r},t_0)$ and current scan $s(\mathit{r},t_k)$ of the phantom in Figure 2 are shown. Subfigure C shows the absolute difference between A and B. Observe that the moving air is easily distinguished in C. Computing the mask ${\Omega }_{\mathrm{S}}$ is performed by applying an appropriate threshold to C such that all moving air is contained in ${\Omega }_{\mathrm{S}}.$

3. Methods

In this section, we start by defining the thermometry acquisition. Hereafter, our TSM air-motion artifact correction scheme is introduced. Next, a heated phantom experiment is detailed, which is used to validate our method against temperature probes. Last, an in vivo validation using data from a hyperthermia cancer treatment is described.

3.1. Thermometry acquisition and B₀-drift correction

The phantom and in vivo temperature is monitored with PRFS-based thermometry acquired by a GE 1.5 T MRI scanner. The MR sequence is a double-echo gradient recalled echo with the first echo at 4.8 ms and the second echo at 19.1 ms. A double-echo scan is used to correct for the electrical conductivity bias [25]. The remaining parameters are: repetition time 620 ms, slice thickness 1 cm, number of slices 25, field of view 50 cm by 50 cm, acquisition matrix 128 × 128, flip angle ${40}^{°},$ and scan time 82 s. The thermometry acquisition interval is 93 s for the phantom and approximately 10 min in vivo.

The B₀-drift is compensated by fitting a 3D polynomial, that is linear in all three dimensions, through four oil-filled reference tubes on the exterior of the phantom and through the reference tubes and internal fat in vivo. Both our TSM method and our implementation of [13] will use the same B₀-drift correction.

3.2. TSM air-motion correction scheme

As mentioned in Section 2.3, computing the mask that includes possible susceptibility sources, is a crucial step in the removal of the air-motion artifacts. Our TSM method aims to reduce the removal of the temperature-induced phase shift from the thermometry, when ${\Omega }_{\mathrm{S}}$ is invariably imperfect.

The TSM data processing pipeline is shown in Figure 4. Here, $\Delta T_{\text{measured}}$ and $\Delta T_{\text{corrected}}$ denote the measured temperature and the air-motion artifact corrected thermometry, respectively. Additionally, $\Delta T_{\text{pred}}$ and $\Delta T_{\text{est}}$ denote the model-based temperature prediction and the estimated temperature field that fuses both the thermometry and the model-based prediction, respectively. To summarize the method in Figure 4, we propose to remove a large part of the temperature-induced phase shift from the thermometry measurement $\Delta T_{\text{measured}}$ using model-based temperature prediction $\Delta T_{\text{pred}}$ prior to applying the background field removal techniques. Loosely speaking, by subtracting the predicted temperature field from the measurement, there is less temperature-induced phase shift in the measurement for the background field removal techniques to inadvertently remove. After the background field removal, the predicted temperature is added back to obtain the corrected measurement $\Delta T_{\text{corrected}}.$ For the remainder of this section, the details of each block in Figure 4 are presented.

Figure 4. Schematic overview of the TSM method. The key contribution is the subtraction of the predicted temperature field prior to the background field removal. To obtain the corrected thermometry the predicted temperate is added back to the result of the background field removal.

3.2.1. ${\mathit{B}}_{\mathit{0}}$-drift correction

The B₀-drift in the measured thermometry is corrected by fitting a linear polynomial through reference objects to obtain $\Delta T_{B_0},$ see Section 3.1 and [7,8].

3.2.2. Prediction

The first step in the model-based observer is the prediction step. Here, the temperature estimate from the previous sample time $t_{k-1}$ is forward simulated using the thermal model to obtain a predicted temperature at time t_k. More specifically, we integrate (1) with heat load $Q_{\text{sar}}(\mathit{r},t_{k-1})$ in (2) over the interval $t_k-t_{k-1}$ starting from the initial state $\Delta T_{\text{est}}(\mathit{r},t_{k-1}).$ We denote this forward simulation by (10) $$\Delta T_{\text{pred}}(\mathit{r},t_k)=f(\mathit{r},\Delta T_{\text{est}}(t_{k-1}),Q_{\text{sar}}(t_{k-1}),t_k-t_{k-1}),$$(10) where $f(·)$ denotes a numerical integration scheme. In our work, we implement $f(·)$ based on a continuous-time state-space model (with a finite-dimensional state) obtained by discretizing the spatial domain [10] and subsequently approximate the matrix exponential, which is central in the solution to a set of coupled linear time-invariant differential equations [26]. Note that the particular method of numerical integration is not an intrinsic part of the presented method and alternative integration schemes can be used.

3.2.3. Background field removal

The background field removal is implemented by solving (8) and substituting $\Delta T_{\text{measured}}(\mathit{r},t_k)$ by $\Delta T_{\text{measured}}(\mathit{r},t_k)-\Delta T_{\text{pred}}(\mathit{r},t_k).$ After solving (8), the resulting susceptibility distribution is subtracted from the measurement to obtain $\Delta T_{\text{corrected}}(\mathit{r},t_k)$ using (7). Note that (7) implicitly accounts for the addition of $\Delta T_{\text{pred}}$ as seen in Figure 4. In Figure 5, the intermediate temperature $\Delta T_{\text{measured}}(\mathit{r},t_k)-\Delta T_{B_0}(\mathit{r},t_k)-\Delta T_{\text{pred}}(\mathit{r},t_k),$ susceptibility field $\Delta {\chi }_{\text{PDF}},$ and identified air-motion-induced artifact $\Delta T_{\chi }$ are shown to provide better insight into the background field removal method.

Figure 5. Transverse slice of the internal temperature and susceptibility fields used in the TSM background field removal method. Subfigure A shows the difference between the B₀-drift corrected thermometry and the predicted temperature, i.e. $\Delta T_{\text{measured}}(\mathit{r},t_k)-\Delta T_{B_0}(\mathit{r},t_k)-\Delta T_{\text{pred}}(\mathit{r},t_k).$ This difference is used to compute the PDF susceptibility distribution. Note that by limiting the PDF method to the difference in measured and predicted temperature, there is less temperature-induced phase shift to inadvertently remove. In subfigure B the identified susceptibility distribution $\Delta {\chi }_{\text{PDF}}(\mathit{r},t_k)$ is shown and in subfigure C the identified susceptibility artifact is shown, i.e. $\Delta T_{\chi }(\mathit{r},t_k)={\alpha }^{-1}\Delta {\chi }_{\text{PDF}}(\mathit{r},t_k)*D(\mathit{r}).$ Observe that that $\Delta {\chi }_{\text{PDF}}$ is confined to a small region defined by ${\Omega }_{\mathrm{S}}$ and note that due to the convolution between the susceptibility distribution and the dipole the impact is observed in larger region than only ${\Omega }_{\mathrm{S}}.$

3.2.4. Fusion

The last step in our TSM observer is fusing the corrected thermometry with the predicted temperature from the thermal model. A simple policy to combine the predicted and measured temperature is (11) $$\begin{array}{ll}\Delta T_{\text{est}}(\mathit{r},t_k)=& (1-w(\mathit{r}))\Delta T_{\text{pred}}(\mathit{r},t_k)\\ & +w(\mathit{r})\Delta T_{\text{corrected}}(\mathit{r},t_k),\end{array}$$(11) with $0\le w(\mathit{r})\le 1.$ Intuitively, $w(\mathit{r})$ reflects the trust in either the model-based prediction or the measurement, e.g., $w(\mathit{r})\approx 0$ in regions with unreliable thermometry and $w(\mathit{r})\approx 1$ in regions with reliable thermometry. For more details on (11), see Appendix B. Last, note that the estimated temperature $\Delta T_{\text{est}}(\mathit{r},t_k)$ will serve as the initial condition for the predicted temperature at time $t_{k+1}.$

3.3. Phantom validation

The TSM air-motion artifact correction scheme is validated using a heated phantom experiment. The phantom is equipped with a movable air volume that can be used to inject air-motion artifacts in a controlled manner, see Figure 2 for a schematic illustration. The movable air volume consists of a 10 mm hollow plastic ball attached to a string and inserted inside a tube filled with distilled water. The air volume is moved in between MR scans to introduce air-motion-induced susceptibility artifacts.

The resulting estimated temperature field $\Delta T_{\text{est}},$ is validated by comparison to four Bowman temperature probes, as seen in Figure 2. The comparison is quantified using the Mean Absolute Error (MAE) and Standard Deviation (STD) with respect to the probe data for $j=1,\dots ,4.$ (12a) $${\text{MAE}}_j:=\frac{1}{n}\sum _{k=1}^n|\Delta T_{\text{probe},j}(t_k)-\Delta T_{\text{est}}({\mathit{r}}_j,t_k)|,$$(12a) (12b) $${\text{STD}}_j:=\sqrt{\frac{1}{n-1}\sum _{k=1}^n(\Delta T_{\text{probe},j}(t_k)-\Delta T_{\text{est}}({\mathit{r}}_j,t_k)-e_j{)}^2},$$(12b) with (12c) $$e_j:=\frac{1}{n}\sum _{k=1}^n\Delta T_{\text{probe},j}(t_k)-\Delta T_{\text{est}}({\mathit{r}}_j,t_k).$$(12c)

Here, $\Delta T_{\text{probe},j}(t_k)$ denotes the relative temperature with respect to time t₀ measured by the j-th probe and ${\mathit{r}}_j$ denotes the position of the j-th probe. Additionally, we compare our method with the current state-of-the-art [13].

Our implementation of the method in [13] will use the same B₀-drift correction. The background field removal is similar, but instead of subtracting the predicted temperature field, the spatially constant and linear terms are subtracted from the B₀-drift corrected thermometry. Additionally, the mask ${\Omega }_{\mathrm{S}}$ is chosen as all voxels in the anatomical scan whose magnitude is below a given threshold, as proposed in [13].

3.4. In vivo validation

Besides the phantom validation, the TSM method is also validated in vivo using thermometry from a hyperthermia cancer treatment. The patient included in the study is a 43-year-old female diagnosed with primary cervical cancer and treated with curative intent using thermoradiotherapy at Erasmus MC. The patient had a histological confirmed cervical carcinoma with FIGO stage IIIb. The research protocol for this investigation was approved by the medical ethics committee of Erasmus MC, University Medical Center Rotterdam. Note that the approval was obtained prior to the start of the study and the ethics code is MEC 2015-108. Additionally, we have received the informed consent from the patient included in this study.

Similar to the phantom, we placed temperature probes in three intraluminal cavities, namely, the rectum, bladder, and vagina. These probes are used to quantify the performance as detailed in (12). The probe locations with respect to the patients anatomy are shown in Figure 6. Note, that the probes in the patient measure the average temperature along the depth, as seen in right hand side from Figure 6. This behavior is taken into account in the validation. Similar to the phantom validation, we also compare our method to the B₀-drift corrected thermometry and the method in [13].

Figure 6. Left, patient anatomy from the transverse plane showing the probe locations. Right, patient anatomy from the sagittal plane through the probes, showing the depth of the probes. The green dashed line indicated the position of the transverse plane. The probes, in numerical order, are located in the vagina, bladder, and rectum, respectively. Note that an average temperature is measured along the length of the probe.

The thermal model of the patient is computed using the segmentation in Figure 7. Here, the patient is segmented with four dominant tissues based on MR-scans at the start of the treatment [27]. The segmentation used by the thermal model is constant throughout the treatment. Note that a static thermal model can lead to model mismatch as the patients anatomy changes with air motion.

Figure 7. Left, patient anatomy in the transverse plane. Right, corresponding segmentation used in the thermal model. Note that the segmentation is based on the anatomy at the start of the treatment.

Figure 8. Transverse slice located at the center of the phantom at 70 min. The subfigures A, B, C, and D denote anatomic scan of the phantom, the B₀-drift corrected thermometry, the air-motion-corrected thermometry using the TSM method, and the estimated temperature field from the observer. Gray areas indicate that no temperature information is available. Note that in B, several air-motion-induced susceptibility artifacts are present (green arrows) that are removed in C and D. The dashed green arrow in A defines the slice along which temperature is sampled for Figure 10.

4. Results

4.1. Phantom validation

The probe measurements together with the results of our method, the method in [13], and the B₀-drift corrected thermometry are shown in Figure 9. A spatial visualization of the thermometry for a single time-step is shown in Figure 8. The performance criteria from (12) applied to the results in Figure 9 are shown in Table 1.

Figure 9. Temperature probe correlation for the phantom experiment with B₀-drift corrected PRFS thermometry, TSM model-based temperature observer, current state-of-the-art [13]. The light gray area denotes the time-window in which heating was applied. Note that probe 4 is used to calibrate $\Delta T_{\text{measured}}-\Delta T_{B_0},$ see Appendix A.

Table 1. Mean absolute error and standard deviation with respect to the probes for the phantom experiment. For details on probe 4, see Appendix A.

MAE (STD) [K]	PRFS	TSM	Wu et al.
probe 1	1.60 (3.37)	0.15 (0.13)	0.36 (0.37)
probe 2	0.14 (0.16)	0.10 (0.09)	0.52 (0.24)
probe 3	0.25 (0.24)	0.07 (0.08)	0.49 (0.24)
probe 4${}^{*}$	0.20 (0.23)	0.08 (0.10)	0.86 (0.25)

The air-motion-induced susceptibility artifacts are observed in the B₀-drift corrected thermometry at probe 1 at approximately 30 and 70 min, see Figure 9. Both our TSM method and the method in [13] mitigate the effect of the susceptibility artifact significantly, with improvements in MAE at probe 1 relative to the B₀-drift corrected thermometry of ten and four times, respectively (Table 1). For the remaining probes, the B₀-drift corrected thermometry is not corrupted by the air-motion-induced susceptibility artifacts, as seen in Figure 9. We observe that our TSM method is consistent with probes 2, 3, and 4 resulting in a MAE of $\le 0.1$ K. The current state-of-the-art yields significantly higher MAE of approximately 0.5 K, as seen by the offset at probes 2, 3, and 4 in Figure 9. Clearly, the method in [13] is susceptible to inadvertently removing the temperature-induced phase shift from the thermometry by the background field removal techniques. Besides the improvements in MAE by the TSM method, we observe a significant reduction in standard deviation, with respect to the B₀-drift corrected thermometry and the method in [13]. This improvement in standard deviation is the result of the spatial and temporal smoothing from the integrated thermal model.

In addition to the spatially sparse probe comparison, the transverse slice in Figure 8 provides a more visual result of our method. In Figure 8B, observe that the thermometry is corrupted by two dominant susceptibility artifacts. In Figure 8C, the susceptibility artifacts are corrected, while minimally affecting the temperature-induced phase shift. In Figure 8D, the observer-based temperature estimate is shown. Note that the observer can estimate the temperature in regions where no thermometry is available and simultaneously improve the effective SNR by spatially and temporally smoothing the measurement.

For a spatial comparison between our TSM method and the method [13], the temperature along the slice indicated in Figure 8A is shown in Figure 10. From Figure 10, it is clear that the method from [13], partly removed the temperature-induced phase shift from the thermometry. In contrast, our TSM method is consistent with the thermometry when sufficiently far from the susceptibility artifact (Figure 10).

Figure 10. The temperature along the dashed arrow from Figure 8A. The y-axis scale is chosen to demonstrate the severity of the air-motion-induced susceptibility artifact (centered at pixel 95). The gray areas denote regions in which no thermometry is available due to, e.g. air or fat voxels. Observe that the TSM method is successful at removing the susceptibility artifact as the estimated temperature field is consistent with the probes and the B₀-drift corrected thermometry in regions not corrupted by the susceptibility artifact.

Figure 11. Transverse slice centrally located in the patient at 60 min. The subfigures B, C, and D denote the B₀-drift correct thermometry, the air-motion corrected thermometry, and the estimated temperature field from the observer. Gray areas indicate that no temperature information is available. Note that in B, a significant air-motion-induced susceptibility artifact is corrupting a large part of the thermometry (green arrow).

4.2. In vivo validation

The probe measurements overlaid with the thermometry are shown in Figure 12. Similar to the phantom experiment, a transverse slice for a single time-step is shown in Figure 11. Last, the MAE and STD are provided in Table 2. Observe that the probes are not located in the transverse slice shown in Figure 11, the particular slice location is chosen as it best demonstrated the susceptibility artifact and the corresponding correction.

Figure 12. Temperature probe correlation for the in vivo validation with B₀-drift corrected PRFS thermometry, TSM model-based temperature observer, current state-of-the-art [13].

Table 2. Mean absolute error and standard deviation with respect to the probes for the patient data.

MAE (STD) [K]	PRFS	TSM	Wu et al.
probe 1, vagina	7.90 (4.01)	0.50 (0.66)	3.12 (1.88)
probe 2, rectum	13.2 (7.31)	0.26 (0.31)	7.23 (2.91)
probe 3, bladder	6.23 (3.58)	1.99 (1.24)	2.40 (1.47)

Notice that, in contrast to the phantom, the in vivo B₀-drift corrected thermometry is continuously corrupted with air-motion-induced susceptibility artifacts as, for example, the intestinal gas is continuously in motion, see, e.g., Figures 11 and 12. As a result, a large region of the thermometry is corrupted by susceptibility artifacts, especially in the vicinity of the intestines. Additionally, due to the longer inter-scan time (approximate 10 min) there is less correlation between subsequent scans, compared to the phantom thermometry. As a result, the temporal smoothing effect of the model-based temperature estimation scheme is reduced.

Based on the quantitative results from Table 2, it is clear that the TSM method outperforms the current state-of-the-art and the B₀-drift corrected thermometry, both in terms of MAE and STD. Most notably, probes 1 and 2 show a good correspondence with a MAE $\le 0.5$ K. The larger discrepancy at probe 3 (MAE $\approx 2$ K) can be the result of model imperfections in combination with unreliable thermometry. Alternatively, bulk motion or imperfect air-motion correction are also possible explanations. Nevertheless, our method outperforms the B₀-drift corrected thermometry at probe 3.

5. Discussion

The phantom with movable air volume provides a controlled experimental setup to verify air-motion artifact correction schemes. This controlled setup clearly showed that the current state-of-the-art potentially removes temperature-induced phase shift from the thermometry. This behavior is also suggested by the underlying theory. From Figure 9, it is clear that augmenting the background field removal techniques with temperature field estimates yields a more robust method that is less susceptible to inadvertent removal of part of the temperature-induced phase shift.

The performance of our TSM method correlates well with the probe data. Note, however, that this comparison is limited to spatially sparse locations. Nevertheless, we believe that the estimated temperature $\Delta T_{\text{est}},$ using our TSM method is accurate over the whole phantom. The main argument for this, is that the estimated temperature field is consistent with the low spatial frequencies of the thermometry when sufficiently far from susceptibility artifacts, see, e.g., Figure 10 between pixels 20 and 70. This suggest that the background field removal techniques in combination with our model-based state observer prevents removing the temperature-induced phase shift from the thermometry measurements.

The TSM method was also applied to in vivo thermometry from a hyperthermia cancer treatment. Here, we obtained promising results that suggest more of the temperature-induced phase shift is preserved, while simultaneously removing the susceptibility artifacts when compared to the method in [13]. Nevertheless, there are several aspects that could further improve the quality of the estimated temperature field. First, the selection of ${\Omega }_{\mathrm{S}},$ i.e. the sources of the susceptibility artifacts is performed by comparing the baseline and current anatomical scans. While this works well for the phantom, its application in vivo is less robust. For example, bulk motion could be misidentified as a source for air-motion-induced susceptibility artifacts. Future work on directly detecting the characteristic shape of the susceptibility artifact can improve the robustness of the method. Second, reducing the time in between thermometry scans reduces the impact of limited model quality by increasing the temporal correlation between scans. Third, more sophisticated state observer designs, as seen in [9,10], could result in better temperature estimates from the available measurements. Fourth, more accurate thermal models with, for instance, dynamic segmentation based on the location of the air volume, would further improve the quality of the temperature estimation scheme. Last, in addition to air-motion, in patients with significant bulk movement, we observed significant temperature errors (especially on tissue boundaries), thereby limiting the thermometry quality. Correcting for this type of motion is not included in the TSM model-based temperature estimation scheme. Hence, methods that can compensate for bulk motion can further remove image artifacts from the thermometry. Nevertheless, the in vivo results show that utilizing model-based observers for monitoring hyperthermia treatments is a promising research direction for simultaneous artifact correction and temperature estimation.

6. Conclusion

In this paper, we presented an integrated thermal and susceptibility model to effectively remove the air-motion-induced susceptible artifacts in PRFS thermometry using an observer. This TSM method is experimentally validated with a dedicated phantom that allowed for a controlled injection of susceptibility artifacts and in vivo with thermometry data from a hyperthermia cancer treatment. The results obtained in the phantom agreed well with four temperature probes that serve as a ground truth. Additionally, the TSM method substantially improved the standard deviation of the temperature estimate through spatial and temporal smoothing using the integrated thermal model. Similarly, the in vivo results show that the TSM method improves the MAE with respect to the B₀-drift corrected thermometry and the current state-of-the-art. These results show that, the integrated thermal and magnetic susceptibility modeling using state observers and background field removal methods is a promising research direction that strives to estimate temperature fields from noisy measurements corrupted with motion artifacts.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research is supported by KWF Kankerbestrijding and NWO Domain AES, as part of their joint strategic research programme: Technology for Oncology II. The collaboration project is co-funded by the PPP Allowance made available by Health ∼ Holland, Top Sector Life Sciences & Health, to stimulate public-private partnerships.

Appendix A

PRFS-coefficient calibration with probe 4

As there is significant uncertainty in the temperature sensitivity of Perfax, our soft tissue analog, used in the phantom, we use probe 4 to calibrate α from (3) using the B₀-drift corrected thermometry, i.e. $\Delta T_{\text{measured}}-\Delta T_{B_0}.$ Probe 4 is selected as it is furthest from any air-motion artifact. Note that the calibration is identical between our TSM method and the method in [13] since the B₀-drift correction is identical between them. For the remainder of this paper, the fourth probe is marked with ‘${}^{*}$’ to denote it is used for calibration.

Appendix B

Fusing model-based temperature predictions and thermometry

This appendix will provide more details on fusing the model-based temperature prediction and the susceptibility artifact corrected thermometry. Thereto, our goal is to estimate $\Delta T_{\text{est}}$ such that its uncertainty is lower than the individual uncertainty of both the model-based prediction and the noisy measurement. In this work, we propose a simple fusing scheme (11) to achieve this. As mentioned briefly, we need design $w(\mathit{r},t_k)$ to reflect our goal of reducing the uncertainty. Mathematically, balancing the uncertainties is achieved by (13) $$w(\mathit{r},t_k)=\frac{\text{var}(\Delta T_{\text{pred}}(\mathit{r}))}{\text{var}(\Delta T_{\text{pred}}(\mathit{r}))+\text{var}(\Delta T_{\text{corrected}}(\mathit{r},t_k))},$$(13) for r in water-based tissue/material. In (13), $\text{var}(·)$ denotes the variance. For materials/tissues that have no temperature sensitivity $w(\mathit{r})=0.$ First, observe that (13) ensures $w(\mathit{r})\in [0,1].$ Second, (13) is analogous to the one-dimensional Kalman filter. Loosely speaking, the Kalman filter is an optimal state estimator when the model-based prediction and measurement both have a Gaussian uncertainty distribution.

The variance of $\Delta T_{\text{corrected}}(\mathit{r})$ is estimated from the signal magnitude using (14) $$\begin{array}{ll}\text{var}(\Delta T_{\text{corrected}}(\mathit{r},t_k))=& {\left(\frac{{\sigma }_0}{\gamma \text{TE}\alpha B_0}\right)}^2\\ & \left(\frac{1}{s^2(\mathit{r},t_k)}+\frac{1}{s^2(\mathit{r},t_0)}\right).\end{array}$$(14)

Here, σ₀ denotes the standard deviation of the complex-valued MR signal. Note that we estimate the variance of the thermometry using the absolute MR-signal at both t₀ and t_k. This approximation of the variance is derived from the signal-to-noise ratio model [7].

The variance of the model-based prediction $\Delta T_{\text{pred}}(\mathit{r})$ is chosen as a constant and used to tune the state observer.