Simultaneous estimation of SAR, thermal diffusivity, and damping using periodic power modulation for MRgFUS quality assurance

Abstract

Purpose: A crucial aspect of quality assurance in thermal therapy is periodic demonstration of the heating performance of the device. Existing methods estimate the specific absorption rate (SAR) from the temperature rise after a short power pulse, which yields a biased estimate as thermal diffusion broadens the apparent SAR pattern. To obtain an unbiased estimate, we propose a robust frequency-domain method that simultaneously identifies the SAR as well as the thermal dynamics.

Methods: We propose a method consisting of periodic modulation of the FUS power while recording the response with MR thermometry (MRT). This approach enables unbiased measurements of spatial Fourier coefficients that encode the thermal response. These coefficients are substituted in a generic thermal model to simultaneously estimate the SAR, diffusivity, and damping. The method was tested using a cylindrical phantom and a 3 T clinical MR-HIFU system. Three scenarios with varying modulation strategies are chosen to challenge the method. The results are compared to the well-known power pulse technique.

Results: The thermal diffusivity is estimated at 0.151 mm²s^–1 with a standard deviation of 0.01 mm²s^–1 between six experiments. The SAR estimates are consistent between all experiments and show an excellent signal-to-noise ratio (SNR) compared to the well established power pulse method. The frequency-domain method proved to be insensitive to B₀-drift and non steady-state initial temperature distributions.

Conclusion: The proposed frequency-domain estimation method shows a high SNR and provided reproducible estimates of the SAR and the corresponding thermal diffusivity. The findings suggest that frequency-domain tools can be highly effective at estimating the SAR from (biased) MRT data acquired during periodic power modulation.

Keywords:

• Specific absorbtion rate
• focused ultrasound
• thermometry
• parameter estimation

I. Introduction

Hyperthermia treatments, where a tumor is heated to temperatures ranging between 39 and 45 degrees Celsius for a duration of 60–90 min, are shown to be an effective adjuvant to conventional cancer therapies [1]. The effect of the hyperthermia treatment depends on, amongst others, the applied thermal dose and the increased perfusion. A priori treatment modeling helps to support the hyperthermia treatment as it enables virtual studies to optimize the outcome [2]. Naturally, these virtual studies require accurate applicator and patient models. Additionally, as applicators are becoming more spatially selective, obtaining accurate models is becoming more relevant. For this reason, methods to verify the accuracy of Radio-Frequency electromagnetic or Focused Ultrasound (FUS) applicator models are essential.

A well-established method, which is incorporated in ESHO QA protocols [3,4], to estimate the Specific Absorption Rate (SAR) supplied by an applicator, utilizes a short power pulse and monitors the temperature rise. For sufficiently short time spans and small spatial SAR gradients, the measured temperature rise indeed approximates the SAR well [5]. However, the approximation deteriorates when diffusion and perfusion effects start (significantly) affecting the temperature distribution. As a result, there is a delicate balance between the signal-to-noise ratio (SNR) of the estimated SAR and the validity of the assumption that measured temperature is proportional to the SAR. When the duration of the power pulse is short, diffusion effects are small compared to the temperature rise. However, the temperature increase itself is also small. Extending the duration of the power pulse increases the temperature further, but also increases the effects from diffusion and perfusion. To complicate matters, diffusion effects are typically dominant in the region of interest, i.e., regions with large SAR gradients such as at the focus of a FUS applicator.

Besides estimating the SAR for quality assurance purposes, estimating the thermal (tissue) parameters is of interest [5–11]. For example, thermal diffusivity (the ratio between conductivity, density, and heat capacity) can be estimated using an invasive ‘self-heated thermistor’ and a temperature sensor [6]. Here, at the first location, a constant temperature is maintained by the thermistor. Then, at the second location, a sensor measures the resulting steady-state temperature increase. The thermal diffusivity is then estimated from the temperature difference. However, this technique requires invasive probes and is only applicable to small volumes. Other diffusivity estimation methods utilize a sinusoidal thermal excitation on the skin and estimates the tissue parameters by extracting the amplitude and phase shifts of the measured surface temperature [9,12]. Besides point and surface measurements, advances in Magnetic Resonance Thermometry [13] (MRT) opened new opportunities to non-invasively estimate the thermal conductivity and perfusion [7]. These existing methods require that the spatial distribution of the SAR is known a priori.

In this paper, we use frequency-domain system identification techniques [14–17] to estimate the SAR distribution and the thermal diffusivity of a phantom using a FUS heating device (the phantom lacks perfusion, and, thus, damping). Our method will explicitly account for the thermal dynamics by estimating the coefficients in the Pennes’ Bioheat Equation [18] (PBE). To achieve this, we periodically modulate the FUS power, which enables the use of advanced frequency-domain tools, such as the Local Rational Method [19] (LRM), to remove the influence of disturbances like B₀-drift in MRT and non steady-state initial temperature distributions. We will verify the proposed method by comparing the estimated SAR and thermal diffusivity between three FUS scenarios where we vary the mean power, modulation frequencies, and experiment duration.

II. Theory

A. Thermal model

We consider a thermal model described by the heat equation with damping and a source [20], of which PBE is a particular realization to model heat transport in patients [18]: (1a) $$\dot{T}(r,t)=\alpha (r){\nabla }^{2}T(r,t)+d(r)T(r,t)+b(r)u(t).$$(1a)

The boundary conditions are (1b) $$T(r,t)=T_{BC}(r,t),\hspace{1em}\text{for}r\text{ on boundary}.$$(1b)

Here, T(r,t) [°C], $T_{BC}(r,t)$ [°C], $r\in {\mathbb{R}}^3$ [m], t [s], $\alpha (r)$ [m²s^–1], d(r) [s^–1], b(r) [°CJ^–1], u(t) [W] denote the temperature, prescribed boundary temperature, position vector, time, heterogeneous thermal diffusivity, heterogeneous damping (from perfusion), normalized power deposition profile (NPDP), and scalar input power, respectively. In contrast to the PBE, the left-hand side of (1a) is divided by the density and specific heat capacity to obtain a unique parameterization of the partial differential equation. Additionally, T(r,t) denotes a temperature relative to a steady-state, e.g., room temperature for a phantom, as this is a natural choice when using relative temperature measurements. Through (1b), we specify Dirichlet boundary conditions, where the boundary temperature is described with $T_{BC}(r,t).$ While this may seem strange at first, we will be using thermometry to prescribe the boundary temperature. This isolates our model from the complicated and unknown true boundary conditions of the patient or phantom [15].

When considering temperature measurements using Proton Resonance Frequency Shift (PRFS) thermometry [21], it is natural to discretize (1a) and (1b) on a cartesian grid in space. Hereto, we define a space of interior voxels at locations $\{r_1,\dots ,r_n\}={\mathbb{D}}_I\subset {\mathbb{R}}^3$ and boundary voxels $\{r_{n+1},\dots ,r_{n+n_{BC}}\}={\mathbb{D}}_{BC}\subset {\mathbb{R}}^3.$ We make the distinction between interior and boundary voxels in order to incorporate the boundary conditions in our modeling framework. To translate the spatially discrete thermal dynamics to a state-space framework, we stack the temperature at the discrete spatial locations in a vector, i.e., $z(t)=\left[\begin{array}{c}{z}_I(t)\\ z_{BC}(t)\end{array}\right],$ with (2) $$z_I(t)=\left[\begin{array}{c}T(r_1,t)\\ ⋮\\ T(r_n,t)\end{array}\right],z_{BC}(t)=\left[\begin{array}{c}T(r_{n+1},t)\\ ⋮\\ T(r_{n+n_{BC}},t)\end{array}\right].$$(2)

Similar to sampling the temperature on the discrete grid, we stack the parameters α, d, and b at the discrete locations $r_i\in {\mathbb{D}}_I$ in a vector $\theta \in {\mathbb{R}}^{3n},$ i.e., (3) $$\theta ={[\alpha (r_1)\dots \alpha (r_n)d(r_1)\dots d(r_n)b(r_1)\dots b(r_n)]}^{\top }.$$(3)

We only estimate the coefficients at the interior voxels, as the boundary voxels are used to provide boundary conditions to our model. Using the newly introduced vectors z and θ, we rewrite (1) in state-space form [22] (4) $${\dot{z}}_I(t)=A(\theta )z(t)+B(\theta )u(t),$$(4)

where the matrix functions $A:{\mathbb{R}}^{3n}\to {\mathbb{R}}^{n\times (n+n_{BC})}$ and $B:{\mathbb{R}}^{3n}\to {\mathbb{R}}^n$ are given in Appendix A. Note that $A(\theta )$ is non-square as it describes the evolution of the interior voxels using both the interior and boundary voxels. While this definition deviates from the typical definition in literature, it is a natural choice for our problem, as we will see below.

We apply the Fourier transform to (4), yielding (5) $$j\omega Z_I(\omega )=A(\theta )Z(\omega )+B(\theta )U(\omega ),$$(5) where $j^2=-1,$ ω denotes the angular frequency, and $Z:\mathbb{R}\to {\mathbb{C}}^{n+n_{BC}},$ and $U:\mathbb{R}\to \mathbb{C}$ denote the state vector and input in the frequency domain, respectively. By applying the Fourier transform to (4), we transformed a differential equation to a complex-valued algebraic equation. As a result, we can estimate θ by solving the algebraic equation for θ given measurements $Z(\omega )$ and $U(\omega ).$ Opposed to measuring $U(\omega ),$ we will obtain $U(\omega )$ by applying a Fast Fourier Transform (FFT) to the prescribed FUS modulation power.

Remark: The SAR, measured in Wkg^–1, is commonly used to quantify the heat load applied to a patient. However, our proposed method estimates, what we call, the normalized power deposition profile (NPDP), measured in °CJ^–1. This choice is motivated by the observation that the SAR is not uniquely identifiable without an additional measurement of the specific heat capacity and FUS power. Nevertheless, the SAR can be derived from the NPDP by multiplying it with the specific heat capacity and the FUS power.

B. Measuring the frequency-domain

In this section, we will detail how $Z(\omega )$ is measured from thermometry data. We will start by specifying requirements for the thermometry, which is followed by the design of the periodic power modulation and the post processing.

We consider temperature measurements using MRT at regularly spaced intervals with sample time t_s, i.e., we measure at times t_k = kt_s with $k\in \mathbb{N}.$ We collect these measurements in a vector $y(t_k)\in {\mathbb{R}}^{n+n_{BC}},$ which is modeled as a combination of the true temperature and noise, (6) $$y(t_k)=z(t_k)+\eta (t_k),{\eta }_i(t_k)\sim \mathcal{N}(\mu (r_i,t_k),\sigma (r_i)).$$(6)

Here, $\eta (t_k)$ is assumed to be normally distributed with mean $\mu (r_i,t_k),$ which denotes a slowly time-varying bias and position dependent variance $\sigma (r_i).$ In PRFS thermometry, the slowly time-varying measurement bias typically represents the B₀-drift.

Given a sequence of measurements at regularly spaced intervals $y(t_1),y(t_2),\dots ,y(t_N),$ we introduce the Fourier transformed measurement vector at discrete frequencies ${\omega }_k\in \Omega =\frac{2\pi }{t_N}\{0,1,\dots ,\frac{N}{2}\}$ as $Y({\omega }_k)\in {\mathbb{C}}^n,$ assuming N to be even. However, $Y({\omega }_k)$ cannot be directly substituted into (5) as a measurement of $Z({\omega }_k),$ as (5) models the thermal dynamics in terms of the steady-state response and $Y({\omega }_k)$ contains a combination of the transient response and the steady-state response. The steady-state response is obtained when the transient dynamics have decayed [10,14] and when there is no B₀-drift. To alleviate this problem, we correct $Y({\omega }_k)$ using the Local Rational Method [14,19] (LRM) to obtain $\overline{Y}({\omega }_k),$ which is a non-biased measurement of $Z({\omega }_k),$ as explained next.

To apply the LRM, we assume that the measurements $Y({\omega }_k)$ have the following structure: $Y({\omega }_k)=Z({\omega }_k)+Y_D({\omega }_k)+Y_N({\omega }_k),$ where Y_D and Y_N denote the transient disturbances and the zero-mean circular complex normally distributed noise. By decomposing the measurements like this, $Y_D({\omega }_k)$ accounts for all transient phenomena, such as B₀-drift, non steady-state initial temperature distributions, and non-periodic applicator power. To isolate $Z({\omega }_k)$ from $Y({\omega }_k),$ we utilize two cornerstone properties of linear time-invariant systems. First, when a linear time-invariant system is excited with a periodic input at a certain frequency, the steady-state response will be at the same frequency. Second, the response to a superposition of inputs is a superposition of the respective individual responses. Hereto, we choose u(t) (the FUS power) as a random phase multi-sine at a limited number of discrete frequencies ${\Omega }_U\subset \Omega ,$ (7) $$u(t)=\sum _{{\omega }_k\in {\Omega }_U}\mathfrak{R}(|U_k|e^{j({\omega }_{k}t+\angle U_k)}).$$(7) Here, $U_k\in \mathbb{C}$ denotes the complex-valued Fourier coefficient at ${\omega }_k\in {\Omega }_U$ and $|·|$ and $\angle$ denote the magnitude and angle, respectively. The resulting frequency-domain representation (for positive frequencies) of u(t) is then (8) $$U({\omega }_k)=\{\begin{array}{cc}{U}_k,\hfill & \text{if}{\omega }_k\in {\Omega }_U,\hfill \\ 0,\hfill & \text{otherwise}.\hfill \end{array}$$(8)

As expected, $U({\omega }_k)$ is only non-zero at the modulation frequencies Ω_U. The effect of this particular modulation strategy is that the thermometry in the frequency-domain is given by (9) $$Y({\omega }_k)=\{\begin{array}{cc}Z({\omega }_k)+Y_D({\omega }_k)+Y_N({\omega }_k),\hfill & \text{if}{\omega }_k\in {\Omega }_U,\hfill \\ Y_D({\omega }_k)+Y_N({\omega }_k),\hfill & \text{otherwise}.\hfill \end{array}$$(9)

Indeed, we expect to see the steady-state response only at the modulation frequencies, while measurement noise and transient disturbances are spread over all frequencies. The key insight exploited by the LRM is that we can estimate $Y_D({\omega }_k)$ at ${\omega }_k\in \Omega$ (including at Ω_U) based on neighboring frequencies that are not excited, see (9). Hereto, we estimate $Y_D({\omega }_k)$ at ${\omega }_k\in \Omega$ by fitting a low order rational function through $Y({\omega }_k)$ in the window ${\Omega }_{k,m}=\{{\omega }_{k-m},\dots ,{\omega }_{k+m}\}\backslash {\Omega }_U,$ i.e., $Y_D({\omega }_k)\approx \frac{{\sum }_{i=0}^l{p}_i({\omega }_k){\omega }_k^i}{{\sum }_{i=0}^l{q}_i({\omega }_k){\omega }_k^i},$ where l is the order of the rational function and $p_i({\omega }_k),q_i({\omega }_k)\in {\mathbb{C}}^{n+n_{BC}}$ are the frequency-dependent vectors of polynomial coefficients. These coefficients are computed, for each frequency ${\omega }_k\in \Omega ,$ by solving the optimization problem (10) $$\underset{p_0,q_0\cdots p_l,q_l}{\arg \min }\sum _{{\omega }_h\in {\Omega }_{k,m}}{\left(Y({\omega }_h)-\frac{\sum _{i=0}^l{p}_i({\omega }_k){\omega }_h^i}{\sum _{i=0}^l{q}_i({\omega }_k){\omega }_h^i}\right)}^2.$$(10)

We solved (10) and computed $Y_D({\omega }_k)$ using the rkfit toolbox [23]. After removing the transient $Y_D({\omega }_k)$ from $Y({\omega }_k),$ we obtain (11) $$\overline{Y}({\omega }_k)=\{\begin{array}{cc}Z({\omega }_k)+Y_N({\omega }_k),\hfill & \text{if}{\omega }_k\in {\Omega }_U,\hfill \\ Y_N({\omega }_k),\hfill & \text{otherwise}.\hfill \end{array}$$(11)

Here, $\overline{Y}({\omega }_k)$ denotes a noisy measurement of the steady-state response, which can be substituted into (5)(5) $$j\omega Z_I(\omega )=A(\theta )Z(\omega )+B(\theta )U(\omega ),$$(5) for ${\omega }_k\in {\Omega }_U$ in order to estimate θ. Based on (11)(11) $$\overline{Y}({\omega }_k)=\{\begin{array}{cc}Z({\omega }_k)+Y_N({\omega }_k),\hfill & \text{if}{\omega }_k\in {\Omega }_U,\hfill \\ Y_N({\omega }_k),\hfill & \text{otherwise}.\hfill \end{array}$$(11) , we estimate the noise variance of $\overline{Y}({\omega }_k)$ using the frequencies that are not excited by the periodic modulation, i.e., (12) $$\text{diag}(\text{var}(\overline{Y}({\omega }_k)))\approx \frac{1}{s_k}\sum _{{\omega }_h\in {\Omega }_{k,m}}\overline{Y}({\omega }_h),$$(12) (13) $$s_k=\text{card}({\Omega }_{k,m})-2(l+1)-1,$$(13) where $\text{card}(·)$ denotes the cardinality of a set. As we will see in the next section, we use the measured forced steady-state response and estimated noise variance to formulate the maximum likelihood estimator for θ.

C. Errors in variables identification

In this section, we present the optimization method we use to estimate θ from noisy LRM corrected measurements $\overline{Y}({\omega }_k).$ Hereto, we start by defining an estimation error in the frequency-domain based on (5)(5) $$j\omega Z_I(\omega )=A(\theta )Z(\omega )+B(\theta )U(\omega ),$$(5) , (14) $$e({\omega }_k,\theta )=A(\theta )\overline{Y}({\omega }_k)+B(\theta )U({\omega }_k)-j{\omega }_k{\overline{Y}}_I({\omega }_k),$$(14) for ${\omega }_k\in {\Omega }_U.$ Loosely speaking, $e({\omega }_k,\theta )$ represents how well the dynamics, that depend on the parameter vector θ, match the measurements. A crucial property is that $e({\omega }_k,\theta )$ is linear in θ, which is immediate from $A(\theta )$ and $B(\theta )$ being linear in θ (see Appendex A). Based on (14)(14) $$e({\omega }_k,\theta )=A(\theta )\overline{Y}({\omega }_k)+B(\theta )U({\omega }_k)-j{\omega }_k{\overline{Y}}_I({\omega }_k),$$(14) , we introduce the Maximum Likelihood (ML) estimator for θ, (15) $${\theta }^{\star }=\underset{\theta \in {\mathbb{R}}^{3n}}{\arg \min }\hspace{1em}\sum _{{\omega }_k\in {\Omega }_U}{e}^{\mathsf{H}}({\omega }_k,\theta )C^{-1}({\omega }_k,\theta )e({\omega }_k,\theta ),$$(15) where ${\theta }^{\star }$ denotes the parameter estimate and $C:\Omega \times {\mathbb{R}}^{3n}\to {\mathbb{R}}^{n\times n}$ denotes the covariance matrix of $e({\omega }_k,\theta ).$ The matrix C is θ-dependent, as noise in $\overline{Y}({\omega }_k)$ is multiplied by $A(\theta )$ and is given by (16) $$C({\omega }_k,\theta )=(A(\theta )-j{\omega }_{k}I)\text{var}(\overline{Y}({\omega }_k)){(A(\theta )-j{\omega }_{k}I)}^{\mathsf{H}}.$$(16)

We solve (15)(15) $${\theta }^{\star }=\underset{\theta \in {\mathbb{R}}^{3n}}{\arg \min }\hspace{1em}\sum _{{\omega }_k\in {\Omega }_U}{e}^{\mathsf{H}}({\omega }_k,\theta )C^{-1}({\omega }_k,\theta )e({\omega }_k,\theta ),$$(15) using gradient descend, where the gradient is replaced by the easier-to-compute pseudo gradient [24].

III. Methods

A. Setup

The parameter estimation method was assessed using a cylindrical phantom in a clinical MR-HIFU system (3 T Achieva, Philips Healthcare, Best, The Netherlands, and Sonalleve V2 HIFU, Profound Medical, Mississauga, Canada), see Figures 1 and 2 for a scan and schematic illustration of the setup.

Figure 1. MRI scan showing a sagittal slice of the phantom placed on top of the FUS applicator. The position of the seven slices that are acquired every 2.5 s are indicated by the yellow rectangles. The approximate acoustic beam path is indicated by the red circular sector.

Figure 2. Schematic overview of the experiment setup. Bottom left: the modulated FUS power, which serves as the input to the system. Right: a sagittal section of the phantom placed on top of the the FUS applicator. Top left: the PRFS thermometry, which serves as the output of the system.

The phantom temperature was measured every 2.5 s using PRFS thermometry, which consisted of seven sagittal multi-shot EPI accelerated scans with a 112 × 112 reconstruction matrix, an EPI acceleration factor of 9, voxel size of $1.8\times 1.8\times 6$ mm, and an effective echo time of 16 ms. Six power modulation experiments, split into three scenarios, were performed to challenge the proposed method across a range of modulation frequencies, experiment duration, and ultrasound powers, see Table 1. Scenario A and C contain two and three experiments, respectively, to demonstrate the consistency under the same experimental conditions. All experiments were performed in succession, allowing the phantom to only partially cool down (approximately five to ten minutes between experiments).

Table 1. Experimental design for three scenarios.

Scenario	A	B	C
Duration [min]	48	48	8
ΩU [rad^–1]	$\frac{2\pi }{720}\{1,3\}$	$\frac{2\pi }{720}\{10,14\}$	$\frac{2\pi }{120}$
ΩU [mHz]	$\{1.39,4.17\}$	$\{13.9,19.4\}$	8.3
Mean power [W]	2	2	3.5
Quantity	2	1	3

The proposed frequency domain method is compared to the power pulse method. For this method, we start with a baseline scan when the phantom is at room temperature. Then, a 25 s pulse of 4 W is applied to the phantom. At 25 s, a second scan is acquired to compute the temperature with respect to the baseline. The NPDP is then estimated by dividing the temperature rise over the time span (25 s) and the applied power (4 W). Both scans are the same multi-shot EPI scans as used for the frequency domain method.

The FUS power modulation signals u(t), and their respective spectra, are seen in Figure 3. Not all frequencies in the u(t) spectrum are sufficiently excited to yield a good SNR, see Figure 3. We discarded the excitation frequencies (empty markers) for which the heating was not clearly visible above the noise. Table 1 lists the excitation frequencies with sufficient SNR.

Figure 3. Time and frequency-domain representation of u(t) for scenarios A (–––), B (–––), C (–––). for scenarios A and B, only part of the input signal is shown. The filled markers in the frequency-domain indicate the frequencies in Ω_U and the empty markers denote the discarded frequencies due to insufficient SNR.

Remark: The SNR can be improved by designing u(t) such that corresponding spectrum concentrates all energy at a few frequencies. This avoids discarding frequencies with insufficient SNR.

B. Data post-processing

Image translation in the frequency encode direction, resulting from a B₀-drift, is corrected using the measured phase of the central k-space voxel [25]. We only correct the image translation resulting from the B₀-drift, the bias that is introduced in the thermometry is inherently corrected by the LRM method, see Section II B. Hereafter, in image-space, the phase is unwrapped and scaled accordingly to obtain the relative temperature measurements from PRFS thermometry [21]. The combined domain of interior and boundary voxels is chosen as the set of voxels with a standard deviation less than 0.3°C. The threshold was motivated by highly spatially correlated noise around the edges of the phantom. The boundary voxels are subsequently defined along the boundary of the previously defined set (naturally, the boundary voxels include the first and last slice of the thermometry). Based on the interior and boundary voxels, a measurement vector $y(t_k)$ according to (6) is constructed for all time-steps.

The frequency domain temperature measurement $\{Y({\omega }_1),\dots ,Y({\omega }_{N/2})\}$ is compensated for the non-instantaneous (delayed) MRI scans, see Figure 4. The correction is given by $\{e^{-j\tau {\omega }_1}Y({\omega }_1),\dots ,e^{-j\tau {\omega }_{N/2}}Y({\omega }_{N/2})\},$ where $\tau =0.75t_s.$ The frequency domain input $\{U({\omega }_1),\dots ,U({\omega }_{N/2})\}$ is multiplied by a zero-order hold as the signal is saved as opposed to being measured in a band limited setting [14]. The correction is given by $\text{sinc}\left(\frac{{\omega }_k}{2\pi }{t}_s\right)e^{-j\frac{{\omega }_k}{2}{t}_s}U({\omega }_k)$ for ${\omega }_k\in {\Omega }_U.$

Figure 4. Illustrating the effective measurement time concept for delayed non-instantaneous measurements.

The LRM correction interpolates $Y_D({\omega }_k)$ over a window of width 20 (excluding frequencies that are modulated), i.e., ${\Omega }_{k,10}=\{{\omega }_{k-10},\dots {\omega }_{k+10}\}\backslash {\Omega }_U,$ with a first-order rational function (l = 1) to obtain, $\overline{Y}({\omega }_k)$ and $\text{var}(\overline{Y}({\omega }_k))$ for ${\omega }_k\in {\Omega }_U$ according to (11)(11) $$\overline{Y}({\omega }_k)=\{\begin{array}{cc}Z({\omega }_k)+Y_N({\omega }_k),\hfill & \text{if}{\omega }_k\in {\Omega }_U,\hfill \\ Y_N({\omega }_k),\hfill & \text{otherwise}.\hfill \end{array}$$(11) and (12)(12) $$\text{diag}(\text{var}(\overline{Y}({\omega }_k)))\approx \frac{1}{s_k}\sum _{{\omega }_h\in {\Omega }_{k,m}}\overline{Y}({\omega }_h),$$(12) , see Figure 5.

Figure 5. Visual demonstration of the LRM method for scenario A1. Top: frequency-domain representation of the measurement, with $Y({\omega }_k)$ (–––) and the estimated transient $Y_D({\omega }_k)$ (–––). Bottom: uncorrected measurement $y(t_k)$ (–––) and corrected measurement ${\mathcal{F}}^{-1}(\overline{Y}(\omega ))$ (–––), where ${\mathcal{F}}^{-1}$ denotes the inverse Fourier transform. As expected, the LRM corrected measurement automatically compensates for the strong drift by focusing on the periodic excitation frequencies.

The following constraints are added to the optimization problem. The damping is set to zero, i.e., $d(r_i)=0$ for all $r_i\in {\mathbb{D}}_I,$ as the phantom lacks perfusion. The thermal diffusivity cannot be (reliably) estimated in regions without heating, for this reason, we constrain the diffusivity to be constant over the interior voxels, i.e., $\alpha (r_i)=\alpha (r_k)$ for all $r_i,r_k\in {\mathbb{D}}_I.$ Last, the optimization problem (15)(15) $${\theta }^{\star }=\underset{\theta \in {\mathbb{R}}^{3n}}{\arg \min }\hspace{1em}\sum _{{\omega }_k\in {\Omega }_U}{e}^{\mathsf{H}}({\omega }_k,\theta )C^{-1}({\omega }_k,\theta )e({\omega }_k,\theta ),$$(15) is solved using the frequencies in Ω_U.

IV. Results

Figure 3 (top) shows the FUS power modulation, u(t), for each scenario. Figure 3 (bottom), shows the FUS power modulation in the frequency-domain, $U({\omega }_k),$ for each scenario. The filled markers denote the frequencies in Ω_U, which are used in the estimation algorithm. The empty markers are discarded frequencies due to insufficient SNR. Figure 6 shows the PRFS thermometry at the focal point along the transducer axis over time, for each experiment. The periodic steady-state response is visible in the time series data, but it is dominated by thermal transient dynamics and B₀-drift.

Figure 6. Temperature at the focal point along the transducer axis over time (see blue arrow Figure 7). The different graphs correspond to experiments A1 (–––), A2 (– – –), B1 (–––), C1 (–––), C2 (– – –), and C3 (–- –- –-), see Table 1 for the definitions. The negative trend in the top plot is predominantly caused by B₀-drift. This drift is not explicitly corrected as the LRM method is insensitive to slowly time-varying measurement biases.

Figure 5 (top) demonstrates the LRM at the focal point on the transducer axis for experiment A1. Distinct ‘peaks’ at the FUS modulation frequencies ${\omega }_k\in {\Omega }_U$ are clearly seen. Additionally, the first-order rational function describes $Y_D({\omega }_k)$ well at the frequencies not excited by the FUS power modulation. Figure 5 (bottom) illustrates the LRM method in the time-domain. The time-domain response suggests that the LRM successfully suppressed transient phenomena (such as the B₀-drift) and extracts the steady-state response from the unprocessed thermometry.

Figure 7 illustrates the LRM and how the results can be interpreted to estimate the thermal diffusivity and NPDP. Figure 7 (left) shows the measured temperature on the central slice with three colored arrows at different distances to the transducer axis. The arrow colors correspond to the colors in the time series. Figure 7 (middle) shows the measured temperature over time at different distances from the transducer axis. The periodic response is clearly visible, as well as a decreasing amplitude and delayed response with an increasing distance from the transducer axis. Figure 7 (right) shows the LRM corrected measurements where transient phenomena are suppressed and noise is reduced. The change in amplitude and phase shifts are used to separate diffusion effects from direct heating.

Figure 7. Magnified temperature map of the focal region for experiment C1 (left). The measured temperature at the locations indicated by the corresponding colored arrows (middle). The LRM corrected temperature transformed back into the time-domain for the corresponding locations (right).

Figure 8 shows the estimated NPDP (which is proportional to the SAR) at slices two until six for experiment A1. Power deposition is detected in all five slices. Additionally, in the middle slice, a clear focus, including near and far field heating, is observed. Recall that slices one and seven are used as a boundary condition and no power deposition is estimated in these slices. Figure 9 compares the estimated NPDP at the central slice for scenarios A1, B1, C1, and the power pulse method, respectively. The NPDP is consistent between all scenarios. Additionally, the overall shape of the estimated NPDP aligns well with the estimate resulting from the power pulse method.

Figure 8. From left to right: the estimated NPDP in slices 2–6 for scenario A1.

Figure 9. Estimated NPDP, from left to right: scenario A1 (duration: long; frequency: low; power: low), scenario B1 (duration: long; frequency: high; power: low), scenario C1 (duration: short; frequency: medium; power: high), and the power pulse method.

Table 2 lists the estimated thermal diffusivity and several properties of the estimated NPDP. The standard deviation of the estimated thermal diffusivity is 0.0088 [mm²s^–1] between all experiments from scenarios A, B, and C. When focusing on the experiments from scenario C, the standard deviation is 0.0012 [mm²s^–1]. This drop in standard deviation could suggest that residual biases are a significant factor in the quality of the resulting estimate. The peak power deposition for all experiments is approximately 20% higher compared to peak power deposition as estimated by the power pulse method. This is expected as diffusion broadens the apparent power deposition profile. The broadening of the focus is also observed when comparing the TC25 and TC50 (number of voxels with at least 25%, or 50%, of the peak value) between the proposed method and the power pulse method. Finally, the spatial standard deviation of the NPDP is between 5 and 10 times lower for the proposed method, compared to the power pulse method.

Table 2. Estimated thermal diffusivity for each scenario and noise standard deviation of the estimated NPDP.

Scenario	A1	A2	B1
α [mm^2s^–1]	0.138	0.142	0.153
max(b) [°CJ^–1]	$801·{10}^{-4}$	$888·{10}^{-4}$	$833·{10}^{-4}$
std(b) [°CJ^–1]	$5.02·{10}^{-4}$	$4.67·{10}^{-4}$	$8.13·{10}^{-4}$
TC25 [voxel]	47	41	40
TC50 [voxel]	17	16	16
Scenario	C1	C2	C3
α [mm^2s^–1]	0.157	0.159	0.157
max(b) [°CJ^–1]	$972·{10}^{-4}$	$900·{10}^{-4}$	$855·{10}^{-4}$
std(b) [°CJ^–1]	$9.62·{10}^{-4}$	$9.83·{10}^{-4}$	$9.69·{10}^{-4}$
TC25 [voxel]	79	91	90
TC50 [voxel]	20	20	20
Scenario	Power pulse
max(b) [°CJ^– 1]	$719·{10}^{-4}$
std(b) [°CJ^– 1]	$48.1·{10}^{-4}$
TC25 [voxel]	142
TC50 [voxel]	28

V. Discussion

In this section, we will discuss the accuracy and precision of our estimated thermal diffusivity and NPDP. This is followed by a brief discussion on experiment design and future applications.

A. Notes on accuracy and precision

The discrepancy between the standard deviation when considering all experiments or just the experiments in scenario C, could suggest that biases contribute more to the error compared to noise. We will discuss potential mechanisms that can introduce a bias in our methodology.

LRM interpolation

A crucial step in the method is the LRM correction, which suppresses non-periodic transients from the thermometry. In doing so, it is assumed that the transient disturbances are smooth in the frequency-domain and can be locally described by a low order rational function (first-order in our setting). Figure 5 shows that this is a reasonable assumption, however, any interpolation error is directly propagated through the estimation algorithm. Hence, this process must be performed with care. Based on our experience, interpolating the transient disturbances in the frequency-domain is more difficult at lower frequencies due to the larger curvature, see, e.g., Figure 5 (top). This effectively limits the lowest possible power modulation frequency (at least three or four periods must be measured for the LRM to be successful).

Estimating the Laplacian

As seen in (14), the error measure encodes the thermal dynamics through the matrices $A(\theta )$ and $B(\theta ).$ Loosely speaking, $A(\theta )$ is the linear map that computes the Laplacian of the temperature field, which is needed to estimate the amount of diffusion at all discrete locations. As the feature size of the FUS applicator is similar in size to the voxels (especially in the out of plane direction), computing the Laplacian through finite differences yields a biased result. Lower power modulation frequencies or smaller voxel sizes are possible measures to reduce this effect. Another interesting approach would be to infer a spatially continuous temperature field from the discrete measurements, using, e.g., Gaussian processes [15]. Such methods derive the Laplacian analytically and do not suffer from the aforementioned limitations, provided that an underlying temperature field can be reconstructed from the thermometry.

Measurement timing uncertainty

Phase shifts in the periodic thermal response are used to differentiate diffusion effects from the NPDP (see Figure 7). For this reason, unaccounted phase shifts originating from timing uncertainties can bias the results. To minimize this problem, we used short scan times and estimated the effective measurement time τ, see Figure 4. However, the exact delay is difficult to estimate as it varies over the volume and is a combination of MR acquisition time and the delay between the commanded FUS power and applied FUS power. Accurately synchronizing the FUS power update intervals with the MR scans and time stamping individual slices will reduce the uncertainty in the (slice dependent) effective measurement time.

Band limited assumption

Band limited experimental conditions are often assumed in the system identification literature, i.e., it is assumed that there is no spectral content above the Nyquist frequency [14]. Loosely speaking, the band limited assumption guarantees that no high frequency spectral content can alias onto the lower frequencies. This condition is typically satisfied through the use of anti-alias filters. However, it is not straightforward to implement such filters for MRI scanners. To minimize the effect of aliasing, it is preferred to externally measuring the FUS power, instead of saving the commanded value. Besides measuring the applied FUS power, fast MR scans are crucial to reduce aliasing effects in the thermometry. In general, thermal dynamics are slow and naturally suppress high frequency excitations. Hence, when the Nyquist frequency is sufficiently high, the thermal dynamics themselves serve as an anti-alias filter. Reducing aliasing effects is therefore the second reason, besides reducing timing uncertainties, to use fast MR scans.

B. Experiment design

Expected SNR

Experiment duration, FUS power amplitude, and excitation frequencies influence the noise standard deviation of the estimated NPDP, as shown in Figure 9 and Table 2. Generally speaking, increasing the experiment duration boosts the SNR at the modulation frequencies, which is an important benefit of the frequency-domain approach. This allows the frequency-domain estimation method to work, even in settings with a bad SNR, given a sufficiently long experiment duration. This is in contrast to the well-established power pulse method, for which an increased experiment duration results in a large bias. Another method to increase the SNR is to use higher FUS power at the modulation frequencies. As briefly mentioned, there are benefits in reducing the acquisition time of a single MRI scan. However, faster scans typically reduce the SNR for the respective scan. Nevertheless, when considering a fixed experiment duration, faster scans result in more measurements. As the expected noise covariance in the frequency-domain is proportional to the reciprocal of the number of scans, the lower SNR from a faster scan is approximately compensated by the increase in the number of scans. This relationship holds provided the increasing noise variance for a single scan is approximately inversely proportional to the scan time. Hence, faster scans have the potential to minimally affect the expected SNR for a fixed experiment duration while increasing the Nyquist frequency and reducing biases.

FUS modulation frequency

An important aspect for the estimation quality is the design of the input signal, i.e., choosing the excitation frequencies and their respective amplitudes. It has been shown that each parameter (e.g., the NPDP and thermal diffusivity) has a varying sensitivity at different frequencies [17]. As a result, choosing the FUS modulation frequencies can have an impact on the quality of the parameter estimate. An important tradeoff is that higher modulation frequencies typically yield a lower SNR as thermal systems attenuate high frequencies, see, e.g., the differences between scenarios A and B in Figure 6. On the other hand, low modulation frequencies can result in longer experiment duration, as at least three or four periods must be measured. Hence, the desired experiment duration implicitly limits the lowest modulation frequency. Clearly, optimal experiment design to obtain the desired estimation accuracy is an interesting topic for further research.

Future applications

Besides quality assurance, in vivo applications are of interest for future research. In vivo estimates of the power deposition profile, thermal diffusivity and damping could, for example, enable personalized treatments. However, extending the method to an in vivo setting introduces additional challenges. For example, motion during identification violates the assumptions of the linear time-invariant model (1). It is expected that sensitivity to motion is particularly high for FUS applicators due to the small size of the focal region. Applying the method in a setting with heterogeneous diffusivity and damping parameters is expected to increase the variance of the parameter estimates. This motivates the need for optimized periodic excitations and novel post-processioning techniques to obtain the desired estimation accuracy within a reasonable measurement time.

VI. Conclusion

In this paper, we presented a frequency-domain system identification approach to identify thermal parameters such as the thermal diffusivity, damping, and spatial power deposition profile. In particular, by transitioning from time-domain to frequency-domain identification, in combination with a periodic FUS modulation, we concentrate the signal power at a few frequencies, resulting in a good SNR. Moreover, by using advanced frequency-domain processing techniques, e.g., the LRM, our method is insensitive to B₀-drift and non steady-state initial temperature distributions. Moreover, the presented method explicitly accounts for measurement noise and thermal dynamics including diffusion and perfusion.

We demonstrated the feasibility of the method by estimating the thermal diffusivity and spatial power deposition in a phantom using an MR-HIFU setup. We performed six experiments based on three different scenarios, each with different input signals and experiment duration. The comparison to the power pulse method showed that our method has a significant increase in SNR and an apparent improvement in the estimated power deposition profile at the focal point of the applicator. The thermal diffusivity was estimated with a standard deviation of 0.01 mm²s^–1 between the six experiments. The spatial power deposition clearly showed the near- and far-field heating and a well-defined focus. Moreover, the spatial shape and magnitude of the estimated power deposition profile is consistent between all experiments, which further supports the expected accuracy of the method.

Finally, we provide several rules of thumb that summarize important aspects of experimental design:

• Number of periods: at least three or four.

• An integer number of periods must be measured.

• Temperature must be measured at equidistant time intervals.

• Scans should be as fast as possible, provided the noise increase is approximately proportional to the decrease in scan time.

• SNR of the estimated NPDP is (approximately) proportional to $\frac{\text{SAR}}{\text{modulation frequency}}\sqrt{\frac{\text{number of scans}}{\text{voxel noise}}}.$

Disclosure statement

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

Data availability statement

The data that support the findings of this study are openly available in “4TU” at http://doi.org/10.4121/94fca849-349e-48d6-a971-e40e165ca684

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

Matrix functions A and B

In this appendix, we present the definition of the matrices $A(\theta )$ and $B(\theta )$ that is well-suited to large-scale systems, such as, the three-dimensional heat equation. Hereto, we directly define $A(\theta )Z(\omega )$ and $B(\theta )U(\omega ).$

First, we start by approximating the Laplacian of the temperature field, i.e., ${\nabla }^{2}T(r_i,\omega )$ for $r_i\in {\mathbb{D}}_I.$ We approximate the Laplacian using finite differences, (17) $$\begin{array}{c}{\nabla }^{2}T(r_i)\approx \frac{T(r_i+\Delta r_x)-2T(r_i)+T(r_i-\Delta r_x)}{\Delta x^2}\\ +\frac{T(r_i+\Delta r_y)-2T(r_i)+T(r_i-\Delta r_y)}{\Delta y^2}\\ +\frac{T(r_i+\Delta r_z)-2T(r_i)+T(r_i-\Delta r_z)}{\Delta z^2}.\end{array}$$(17)

Here, $\Delta r_x,\Delta r_y,$ and $\Delta r_z$ denote the distance to neighboring voxel in the x, y, and z directions, respectively. Note that we require $r\in {\mathbb{D}}_{BC}$ to estimate the Laplacian at $r\in {\mathbb{D}}_I.$ Given the Laplacian of the temperature, $A(\theta )Z(\omega )$ is given by (18) $$\begin{array}{c}A(\theta )Z(\omega )=\text{diag}\left(\left[\begin{array}{c}{\nabla }^{2}T(r_1,\omega )\\ ⋮\\ {\nabla }^{2}T(r_n,\omega )\end{array}\right]\right)\left[\begin{array}{c}\alpha (r_1)\\ ⋮\\ \alpha (r_n)\end{array}\right]\\ +\text{diag}\left(\left[\begin{array}{c}T(r_1,\omega )\\ ⋮\\ T(r_n,\omega )\end{array}\right]\right)\left[\begin{array}{c}d(r_1)\\ ⋮\\ d(r_n)\end{array}\right].\end{array}$$(18)

Second, $B(\theta )U(\omega )$ is given by (19) $$B(\theta )U(\omega )=\left[\begin{array}{c}b(r_1)\\ ⋮\\ b(r_n)\end{array}\right]U(\omega ).$$(19)

The Dirichlet boundary condition is captured by estimating ${\nabla }^{2}T(r_i,\omega )$ on $r\in {\mathbb{D}}_I$ using a finite difference scheme. Indeed, we require all voxels in ${\mathbb{D}}_I\cup {\mathbb{D}}_{BC}$ to estimate the Laplacian on ${\mathbb{D}}_I.$