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 mm2s–1 with a standard deviation of 0.01 mm2s–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 B0-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.


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][6][7][8][9][10][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 noninvasively 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][15][16][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 0 -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.

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]: The boundary conditions are Tðr, tÞ ¼ T BC ðr, tÞ, for r on boundary: Here, T(r,t) 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 fr 1 , :::, r n g ¼ D I � R 3 and boundary voxels fr nþ1 , :::, r nþn BC g ¼ D BC � 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 loca- Similar to sampling the temperature on the discrete grid, we stack the parameters a, d, and b at the discrete locations r i 2 D I in a vector h 2 R 3n , i.e., h ¼ aðr 1 Þ ::: 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 h, we rewrite (1) in state-space form [22] where the matrix functions A : R 3n !R n�ðnþn BC Þ and B : R 3n !R n are given in Appendix A. Note that AðhÞ is nonsquare 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 where j 2 ¼ −1, x denotes the angular frequency, and Z : R !C nþnBC , and U : R !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 h by solving the algebraic equation for h given measurements ZðxÞ and UðxÞ: Opposed to measuring UðxÞ, we will obtain UðxÞ 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ðxÞ 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 2 N: We collect these measurements in a vector yðt k Þ 2 R nþn BC , which is modeled as a combination of the true temperature and noise, Here, gðt k Þ is assumed to be normally distributed with mean lðr i , t k Þ, which denotes a slowly time-varying bias and position dependent variance rðr i Þ: In PRFS thermometry, the slowly time-varying measurement bias typically represents the B 0 -drift.
Given a sequence of measurements at regularly spaced intervals yðt 1 Þ, yðt 2 Þ, :::, yðt N Þ, we introduce the Fourier transformed measurement vector at discrete frequencies x k 2 X ¼ 2p tN f0, 1, :::, N 2 g as Yðx k Þ 2 C n , assuming N to be even.However, Yðx k Þ cannot be directly substituted into (5) as a measurement of Zðx k Þ, as (5) models the thermal dynamics in terms of the steady-state response and Yðx 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 0 -drift.To alleviate this problem, we correct Yðx k Þ using the Local Rational Method [14,19] To apply the LRM, we assume that the measurements Yðx k Þ have the following structure: 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 ðx k Þ accounts for all transient phenomena, such as B 0 -drift, non steady-state initial temperature distributions, and non-periodic applicator power.To isolate Zðx k Þ from Yðx k Þ, we utilize two cornerstone properties of linear timeinvariant 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 X U � X, Here, U k 2 C denotes the complex-valued Fourier coefficient at x k 2 X U and j � j and / denote the magnitude and angle, respectively.The resulting frequency-domain representation (for positive frequencies) of u(t) is then As expected, Uðx k Þ is only non-zero at the modulation frequencies X U .The effect of this particular modulation strategy is that the thermometry in the frequency-domain is given by 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 ðx k Þ at x k 2 X (including at X U ) based on neighboring frequencies that are not excited, see (9).Hereto, we estimate , where l is the order of the rational function and p i ðx k Þ, q i ðx k Þ 2 C nþnBC are the frequencydependent vectors of polynomial coefficients.These coefficients are computed, for each frequency x k 2 X, by solving the optimization problem We solved (10) and computed Y D ðx k Þ using the rkfit toolbox [23].After removing the transient Here, � Y ðx k Þ denotes a noisy measurement of the steadystate response, which can be substituted into (5) for x k 2 X U in order to estimate h.Based on (11), we estimate the noise variance of � Y ðx k Þ using the frequencies that are not excited by the periodic modulation, i.e., where 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 h.

C. Errors in variables identification
In this section, we present the optimization method we use to estimate h from noisy LRM corrected measurements � Y ðx k Þ: Hereto, we start by defining an estimation error in the frequency-domain based on (5), for x k 2 X U : Loosely speaking, eðx k , hÞ represents how well the dynamics, that depend on the parameter vector h, match the measurements.A crucial property is that eðx k , hÞ is linear in h, which is immediate from AðhÞ and BðhÞ being linear in h (see Appendex A).Based on ( 14), we introduce the Maximum Likelihood (ML) estimator for h, where h ?denotes the parameter estimate and C : X � R 3n !R n�n denotes the covariance matrix of eðx k , hÞ: We solve (15) using gradient descend, where the gradient is replaced by the easier-to-compute pseudo gradient [24].

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.
The phantom temperature was measured every 2.5 s using PRFS thermometry, which consisted of seven sagittal multishot EPI accelerated scans with a 112 � 112 reconstruction matrix, an EPI acceleration factor of 9, voxel size of 1:8 � 1:8 � 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).
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.
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 0 -drift, is corrected using the measured phase of the central k-space voxel [25].We only correct the image translation resulting from the B 0 -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 fYðx 1 Þ, :::, Yðx N=2 Þg is compensated for the non-instantaneous (delayed) MRI scans, see Figure 4.The correction is given by fe −jsx 1 Yðx 1 Þ, :::, e −jsx N=2 Yðx N=2 Þg, where s ¼ 0:75t s : The frequency domain input fUðx 1 Þ, :::, Uðx N=2 Þg 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 sinc The LRM correction interpolates Y D ðx k Þ over a window of width 20 (excluding frequencies that are modulated), i.e., X k, 10 ¼ fx k−10 , :::x kþ10 gnX U , with a first-order rational function (l ¼ 1) to obtain, � Y ðx k Þ and varð � Y ðx k ÞÞ for x k 2 X U according to (11) and (12), see Figure 5.
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 2 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., aðr i Þ ¼ aðr k Þ for all r i , r k 2 D I : Last, the optimization problem (15) is solved using the frequencies in X U .
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 0 -drift.
Figure 5 (top) demonstrates the LRM at the focal point on the transducer axis for experiment A1.Distinct 'peaks' at the FUS modulation frequencies x k 2 X U are clearly seen.Additionally, the first-order rational function describes Y D ðx 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 0 -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 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 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.
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 2 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 2 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.

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ðhÞ and BðhÞ: Loosely speaking, AðhÞ 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 s, 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.

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 frequencydomain 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 frequencydomain 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 0 -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 2 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 SAR modulation frequency ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi number of scans voxel noise q :

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

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.

Matrix functions A and B
In this appendix, we present the definition of the matrices AðhÞ and BðhÞ that is well-suited to large-scale systems, such as, the three-dimensional heat equation.Hereto, we directly define AðhÞZðxÞ and BðhÞUðxÞ: First, we start by approximating the Laplacian of the temperature field, i.e., r 2 Tðr i , xÞ for r i 2 D I : We approximate the Laplacian using finite differences, Here, Dr x , Dr y , and Dr z denote the distance to neighboring voxel in the x, y, and z directions, respectively.Note that we require r 2 D BC to estimate the Laplacian at r 2 D I : Given the Laplacian of the temperature, AðhÞZðxÞ is given by The Dirichlet boundary condition is captured by estimating r 2 Tðr i , xÞ on r 2 D I using a finite difference scheme.Indeed, we require all voxels in D I [ D BC to estimate the Laplacian on D I :

Figure 1 .
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 .
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.

Figure 3 .
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 X U and the empty markers denote the discarded frequencies due to insufficient SNR.

Figure 4 .
Figure 4. Illustrating the effective measurement time concept for delayed noninstantaneous measurements.

Figure 7 .
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 .
Figure 8. From left to right: the estimated NPDP in slices 2-6 for scenario A1.

Table 1 .
Experimental design for three scenarios.

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