Robust sound speed estimation for ultrasound-based hepatic steatosis assessment

Hepatic steatosis, due to fat accumulation in the liver is the most common cause of chronic liver disease and may lead to severe liver conditions (Wieckowska and Feldstein 2008). Biopsy and MRI are gold standard techniques to diagnose hepatic steatosis with a percentage of fat in the liver. These techniques are not without limitations. First, liver biopsy suffers from sampling problems: liver biopsies sample as little as 1/50 000 of the total mass of the liver, often resulting in insufficient information for a definitive diagnosis (Janiec et al 2005, Ratziu et al 2005). Second, it is an invasive method involving certain risks and added stress and expense. Finally, the histologic evaluation is subjective and dependent on the experience of the pathologist. On the other hand, MRI is currently the most accurate noninvasive imaging method for diagnosing and quantifying hepatic steatosis, which estimates the proton-density fat fraction (PDFF) as a measure of fractional fat content (Leporq et al 2014). However MRI also has several limitations including cost, contra-indications and poor availability that could be overcome by using ultrasound. Although conventional ultrasound is widely used by radiologists for hepatic steatosis assessment, it remains qualitative. Therefore, there is a medical need to develop noninvasive techniques that can robustly quantify the degree of hepatic steatosis.

Ultrasonic imaging has been explored for many years for its ability to detect and characterize liver disease and is highly accurate to diagnose liver cirrhosis (Tchelepi et al 2002, Mishra and Younossi 2007, Deffieux et al 2015). Yet, current conventional ultrasonic techniques do not allow for quantification of the degree of fatty liver. Presence of liver steatosis is depicted on conventional ultrasound when the liver appears bright but this finding has poor sensitivity and requires 30% of hepatic steatosis (Mehta et al 2008, Dasarathy et al 2009). Moreover, this qualitative assessment is highly subjective and depends on the expertise and experience of the operator (Zwiebel 1995). Researchers have investigated the liver-kidney contrast to quantify liver fat content using the ultrasound hepatic/renal ratio and the hepatic attenuation rate from ultrasound hepatic and right kidney images (Xia et al 2012). However, this technique still needs standardization, further testing in a clinical setting and will always suffer from the impact of variable acoustic window. Lastly, the effectiveness of ultrasound for diagnosing hepatic steatosis is reduced in patients who are morbidly obese (de Moura Almeida et al 2008), mainly due to reverberation in abdominal tissue (Lediju et al 2009) and wave front distortion induced by the subcutaneous fat (Browne et al 2005).

Suzuki and coworkers observed that the ultrasonic attenuation depends on fatty infiltration of the liver and to a lesser extent on fibrosis (Suzuki et al 1992). This theory was also used to develop an attenuation parameter based on the ultrasonic properties of the radiofrequency backpropagated signal called Controlled Attenuation Parameter (CAP) (Sasso et al 2010). The CAP has the ability to detect minimal hepatic steatosis such as 10% steatosis assessed on pathology. This CAP technique is performed in conjunction with liver fibrosis assessment based on 1D Transient Elastography (Foucher et al 2006b, Sandrin et al 2002). Other recent studies have used ultrasound to detect steatosis (Lin et al 2015, Son et al 2015) but they are based on a single cutoff value. Beyond such binary below/above diagnosis, there is a real clinical unmet need for a non-expensive, widely available, and highly reliable technique to precisely grade hepatic steatosis as MRI does.

The objective of this study was to assess the accuracy for steatosis detection and staging of a specific ultrasonic sequence for sound speed estimation (SSE). This work presents a method able to precisely calculate the sound speed in the liver. As it is well known that speed of sound within soft tissues varies slightly with fat content, a relationship between sound speed and percentage of fat in the liver can be found. Indeed, an increase in fat content leads to a decrease in wave speed (Duck 1990). Bamber and Hill reported higher mean sound speed in excised normal liver than in fatty human livers (Bamber and Hill 1981). In another in vivo study, researchers reported higher sound speed in normal liver than in fatty liver without fibrosis from humans (Chen et al 1987). However, these studies were performed in excised organs. For non-invasive sound speed measurements, Jaeger and Frenz (2015) proposed to measure the changing local phase of beamformed echoes when changing the transmit beam steering angle. The method developed in our study is based on the study of the spatial coherence function of the backscattered echoes resulting of an ultrasound beam focusing in the medium (Mallart and Fink 1994, Lacefield et al 2002). The optimal sound speed is deduced by increasing the spatial coherence of echoes coming from a targeted focal spot. Before reaching the liver, the ultrasound beam crosses fatty and muscle layers of different thickness. These layers introduce distortions of the ultrasonic wavefront (Browne et al 2005) that hinder the robustness of SSE. In order to correct for these aberrations and account for them in the SSE calculation, firstly, a virtual point source is created in the liver and an iterative algorithm (Montaldo et al 2011) is used for phase aberration correction. By correcting the aberrations suffered by the ultrasonic beams during the first layers of wave propagation, the final SSE obtained is more precise and robust. Nevertheless, it is to a global value corresponding to the integral over the travel path. Secondly, we propose to correct the influence of the superficial layers thickness on the SSE by integrating superficial layers measurement made by the physician into the calculation. This two-steps SSE correction leads to a robust and local SSE in the liver. The technique was first tested on homogeneous and multilayers phantoms. Then a pilot clinical study of 17 patients was conducted with the MRI PDFF and biopsy as gold standards.

2.1. Ultrasound acquisitions

2.2. Experimental set-up on phantoms

2.3. Pilot clinical study design

2.4. Theory

The aim of this study is to obtain a precise estimation of the sound speed in liver. To achieve this goal the speckle noise technique for virtual source generation and aberration correction is used. This technique consists in trying to recreate inside the medium a virtual point-like reflector able to act as a virtual source generating the Green's function between this source location in the medium and the transducer elements of the ultrasonic probe. This iterative method required a reliable first SSE to be able to converge. This first SSE is obtained by using the robust van Cittert–Zernike (VCZ) Theorem applied on the ultrasound backscattered echoes.

2.4.1. Spatial coherence in random media: the van Cittert–Zernike theorem.

The technique used in this study starts with the transmission of an ultrasound pulse focused in the region of interest. Then the ultrasound field backscattered by the random distribution of scatterers is received on all the array elements. The similarity between the signals received by two distant elements of the array characterizes the spatial coherence of the received wavefield. Van Cittert and Zernike determined the degree of coherence by defining a coherence function as the averaged cross-correlation between two signals received at two points of space (Mallart 1991). The Van Cittert–Zernike theorem states that the coherence function is the spatial Fourier transform of the intensity distribution at the focus. By calculating autocorrelations between all pairs of receiver elements, the coherence function R(m) is assessed as a function of distance in number of elements m (or Element lag) (Derode and Fink 1993):

$$\begin{eqnarray}&&R(m)=\frac{N}{N-m}\frac{{\sum}_{i=1}^{N-m}c\left(i,i+m\right)}{\underset{{}}{\overset{{}}{{{\sum}_{i=1}^{N}}}}\,c(i,i)}\end{eqnarray} \tag{ 1 }$$

where $N$ is the number of elements of the array and $c(i,j)$ is defined as:

$$\begin{eqnarray}&&c(i,j)=\underset{{{T}_{1}}}{\overset{{{T}_{2}}}{\sum}}\,\left({{S}_{i}}(t)-{{\bar{S}}_{i}}\right)\left({{S}_{j}}(t)-{{\bar{S}}_{j}}\right)\end{eqnarray} \tag{ 2 }$$

Where $\left[\begin{array}{c} {{T}_{1~}} & {{T}_{2}} \end{array}\right]$ is the temporal window centered on the focal time and Si is the time-delayed signal received on transducer i.

When the ultrasonic wave is focused on a point-like target in the medium, the coherence function is equal to 1 all along the array. In the case of a homogeneous medium made of randomly distributed Rayleigh scatterers, the degree of coherence decreases as the distance between elements increases.

Therefore, a focused beam generated by a rectangular aperture will lead to a triangle coherence function of backscattered echoes coming from a random distribution of Rayleigh scatterers in the focal spot (Mallart and Fink 1994). For a fixed depth, if the speed of sound used for focusing in the homogeneous random medium is under- or overestimated, then the focal spot size at the desired depth will increase, leading to a dramatic decrease of coherence. In this study we choose a fixed depth of 60 mm as it corresponds to the focal elevation depth of probe we used.

The method consists in proposing an algorithm that will try to better focus in the medium in post-processing with different speed of sound. In order to perform any kind of transmit focusing, we acquire the backscattered echoes coming from a set of coded excitations and recreate the desired transmit focusing using a coherent recombination of these signals during the post-processing step. Such coherent recombination of ultrasonic backscattered echoes was in particular studied in the context of dynamic virtual transmit focusing (Cooley and Robinson 1994, Karaman et al 1995, Lockwood et al 1998). For every try, the coherence function corresponding to the tested speed of sound is recorded. The area under the coherence function curve is then calculated (Lediju et al 2011). To improve signal to noise ratio (SNR), the Van Cittert–Zernike theorem is applied on a 25 points grid (points separated by one wavelength, grid at 60 mm depth) and the area under the coherence function curve is averaged for every tested speed of sound. The estimated speed of sound providing the highest area under the curve value is the real cumulative speed of sound in the medium (figure 1) and enables the best focusing quality with the smallest focal spot. This technique gives access to the global SSE in the medium. It is a robust SSE that can be used as the first SSE value in the more precise iterative algorithm developed in the next section.

Zoom In Zoom Out Reset image size

2.4.2. Virtual point-like source generation and iterative focusing algorithm for phase aberration correction.

The attainment of an optimal coherence function is possible in a perfectly homogeneous medium with a constant speed of sound. However, sound speed heterogeneities in the intercostal space of difficult patients result in phase aberrations along the travel path of the ultrasonic beam and hinder the precision and robustness of the method. The aim of the second step of our method is to create a virtual point-like reflector in the homogeneous random medium (under the superficial layers) in order to assess the aberrations induced by the medium and improve the precision of our SSE. Indeed, intercostal layers act as near field phase aberrator. In our study, these aberrations are estimated by creating a virtual bright reflector from a random distribution of scattered below the aberrating layer using the concept of iterative time reversal focusing in speckle noise (Montaldo et al 2011) (figure 2) applied on the same RF data used for the previous step. Twenty-five focus beams are virtually emitted in post-treatment at different locations nearby the desired focus (60 mm depth) by recombination of the RF issued from the Hadamard encoded acquisitions (figure 2(a)). In this study, virtual backscattered signals are recorded (figure 2(b)) and 'steered' using the time delays as if all signals were coming from the same reference location (figure 2(c)). The summed corrected backscattered echoes lead to pulsed signals well compressed in time resulting from the created virtual point-like reflector (figure 2(d)). The interested reader should also note that this aberration correction based on time reversal focusing could be performed directly in blood vessels when available following (Osmanski et al 2012).

Zoom In Zoom Out Reset image size

Through an aberrating layer an iterative algorithm is necessary to converge to a well-defined wave front. In this case, we are going to take benefit on both the isoplanatic patch (Chassat 1989) and the randomness of the scattering medium. When applying the process described earlier to create a point-like reflector in a homogeneous medium, the emitted signals from this virtual bright reflector are no more matching to the heterogeneous medium. So the initial focusing beams are not perfectly focused at their desired locations, however it provides an acceptable start as a first iteration. The new iteration step start replacing the virtually emitted signal from the transducer by the time-reversed version of the signal emitted by the virtual point-like reflector obtain from the first iteration (figure 2(f)). At each iteration the emitted beam profile is improved. After 10 iterations the signal converge to a well-defined wave front. The spatial coherence in random media for the 10th iteration was studied for improving the precision of our SSE.

2.4.3. Superficial layers influence correction.

In addition to the phase aberration correction (heterogeneous local aberrations), superficial layers thickness and speed of sound in these different layers (global averaged aberrations) have to be integrated in the calculation. In this paper, one way of including fat/muscle layers influence correction is studied and discussed: superficial layers thickness measurement.

Thicknesses of fat and muscle layers were measured by the physician with conventional ultrasound. The mean sound speed used in the calculation are 1.450 mm · µs⁻¹ in fat and 1.575 mm · µs⁻¹ in muscle (Azhari 2010). In a multi-layer medium, sound speeds of the different layers are linked by the formula:

$$\begin{eqnarray}&&\frac{{{d}_{\text{tot}}}}{{{c}_{\text{cum}}}}=~\frac{{{d}_{\text{fat}}}}{{{c}_{\text{fat}}}}+~\frac{{{d}_{\text{muscle}}}}{{{c}_{\text{muscle}}}}+~\frac{{{d}_{\text{liver}}}}{{{c}_{\text{liver}}}}\end{eqnarray} \tag{ 3 }$$

Where d_tot is 60 mm, c_cum the global SSE calculate by the algorithm, d_fat the fat layer thickness, c_fat the sound speed in fat, d_muscle the muscle layer thickness, d_liver  =  d_tot  −  (d_fat  +  d_muscle), c_liver the speed of sound in liver that we aim to determine. c_liver is calculated with the formula:

$$\begin{eqnarray}&&{{c}_{\text{liver}}}=\frac{{{d}_{\text{liver}}}}{\frac{{{d}_{\text{tot}}}}{{{c}_{\text{cum}}}}~-~\frac{{{d}_{\text{fat}}}}{{{c}_{\text{fat}}}}~-~\frac{{{d}_{\text{muscle}}}}{{{c}_{\text{muscle}}}}}\end{eqnarray} \tag{ 4 }$$

In phantom experiments, superficial layers were mimicked by layers of water with controlled temperature. The CIRS phantom was used (model 054, speed of sound: 1.54 mm · µs⁻¹). Three different water temperatures (c_water) were tested: 37 degrees (c_water  =  1.52 mm · µs⁻¹), 20 degrees (c_water  =  1.48 mm · µs⁻¹), 14 degrees (c c_water  =  1.46 mm · µs⁻¹) (Boed 1999) with different thicknesses from 0 to 26 mm.

2.5. Statistical analysis

3.1. Sound speed estimation in homogeneous phantoms

A range of sound speed from 1.400 mm · µs⁻¹ to 1.600 mm · µs⁻¹ with a 0.001 step was tested for the ATS phantom (model 551) (figure 3(a)) and for the CIRS phantom (model 054) (figure 3(b)).

Zoom In Zoom Out Reset image size

For each sound speed tested, the maximum correlation was calculated as a function of the lag (in elements). The coherent sum of this maximum correlation was then calculated using equation (2) for each tested sound speed. The SSE value corresponding to the phantom was calculated as the maximum of the polynomial fit of the coherent sum of the maximum correlation. ISO-value curves are another way to visualize the influence of sound speed on the signal coherence. We plotted curves of constant minimum value of the maximum correlation for each sound speed and noticed that the sound speed characterizing the medium is assuring the highest lag for each minimum of the maximum correlation (figures 3 and 4).

Zoom In Zoom Out Reset image size

In order to find the mean SSE and deviation from the mean for both phantoms and to compare these results with the speed of sound indicated by the constructor, 20 acquisitions with probe repositioning were performed on each phantom. The mean SSE for the ATS phantom was found to be 1.449  ±  0.006 mm · µs⁻¹ for a sound speed given by the constructor of 1.450 mm · µs⁻¹. The mean SSE for the CIRS phantom was found to be 1.544  ±  0.003mm · µs⁻¹ for a sound speed given by the constructor of 1.540 mm · µs⁻¹.

These results demonstrated that in the homogeneous phantom case, the obtained SSE was in good agreement with the speed of sound indicated by the constructor.

3.2. First sound speed estimation in patients

3.3. Sound speed estimation improvement

Studying the coherence function for different sound speed give a global SSE. The aim of this section is now to calculate a local speed of sound in the liver based on the previous SSE. There are two ways of SSE improvement: by correcting the phase aberration induced by the fat and muscle superficial layers, and by taking into account the thickness of these layers.

3.3.1. Phase aberration correction.

Fat and muscle superficial layers act as near field phase screen aberrator. The technique used in this section is based on the aberration correction algorithm. The aim is to straighten the wave front coming from the virtual point like reflector to improve the quality of the focusing. In the case of a patient with a thick superficial fat layer, the phase aberration algorithm succeeds in straightening the wave front and improving the focusing (figure 5).

Zoom In Zoom Out Reset image size

The aberration correction algorithm succeeds in improving the focusing quality in every case and in some cases it improves the SSE (figure 6). From iteration 1 to iteration 10, for all patients in this study, a range from 0 to 12% increase in the area under the coherent function curve was measured. For patients with thin superficial layers and mostly composed by muscle, the aberration correction algorithm improves the coherence of the backscattered echoes and confirms the SSE (figure 6, Patient 1). For patients with thick superficial layers and mostly composed by fat, the aberration correction algorithm also improves the coherence of the backscattered echoes but modifies the SSE (figure 6, Patient 2). In both cases, this demonstrate that the first SSE calculated with the VCZ theorem is a robust value leading to the algorithm convergence. This convergence is illustrated by the increase of the area under the VCZ curve (figure 6(c)).

Zoom In Zoom Out Reset image size

3.3.2. Superficial layers influence correction.

In this section, one way of including superficial layer thickness into the sound speed calculation is detailed: by superficial layers thicknesses measurement and with a good knowledge of the sound speed in these layers.

To validate the multilayers model in the phantom experiment, we use the thickness of the superficial layer (water at different temperatures) and the sound speed in this layer to correct the SSE. Without layer correction, SSE varies depending on the sound speed in the superficial layer and the thickness of this layer (figure 7(a)). After layer correction the sound speed estimation corresponds to the sound speed given by the constructor (1.54 mm · µs⁻¹) without any dependence with the superficial layer characteristic (figure 7(b)).

Zoom In Zoom Out Reset image size

3.4. Pilot clinical study results

The aim of this section was to challenge the different steps of SSE calculation and to assess its robustness with two different gold standards, biopsy and MRI PDFF.

Ultrasound and MRI were undergone the same day and MRI and biopsy were undergone with a median of 6 months in between. The first SSE given by the VCZ theorem already gives a proportional relation between MRI PDFF and SSE (R²  =  0.595) (figure 8(a)). As our hypothesis was to find a correlation between the sound speed in the liver (SSE) and the percentage of fat in the latter, we used the coefficient of determination R² as a relative comparison tool between the different improvement steps. We noticed that aberration correction and superficial layer thickness inclusion in the calculation both improve the proportional relation between MRI PDFF and SSE (respectively figures 8(b) and (c)). With the aberration correction step, aberrations are better taken into account, leading to a higher spatiotemporal coherence of backscattered signals (see figure 6). Consequently, when adding the superficial layers correction, we reached the highest agreement we can obtain so far between MRI PDFF and SEE for this patient cohort (R²  =  0.691) (figure 8(d)).

Zoom In Zoom Out Reset image size

With the MRI PDFF technique, the cut-off value between healthy and diseased patients was 5% PDFF (Idilman et al 2013) (corresponding to the 10% used in biopsy). SSE was able to significantly differentiate healthy and diseased patients (p  <  0.0001, AUROC  =  0.942) with an associated criterion of 1.541 mm · µs⁻¹.

These results are confirmed by the comparison between biopsy and SSE (figure 9(a)). This pilot study demonstrates the SSE ability to significantly differentiate healthy and diseased patients (p  <  0.0001, AUROC  =  0.952, figure 9) with an associated criterion of 1.555mm · µs⁻¹.

Zoom In Zoom Out Reset image size

In this study, we introduced a method for precise and robust in vivo sound speed assessment. The technique was firstly validated on two phantoms. In the case of homogeneous phantom the SSE found was in good agreement with the speed of sound indicated by the constructor (Phantom 1, c_construcor  =  1.450 mm · µs⁻¹ and SSE  =  1.449  ±  0.006 mm · µs⁻¹; Phantom 2, c_construcor  =1.540 mm · µs⁻¹ and SSE  =  1.544  ±  0.003 mm · µs⁻¹).

Secondly, in the pilot clinical study, the potential of this method for noninvasive liver steatosis detection and staging was proved. Our approach enables an accurate estimate of the ultrasonic sound speed in the liver with a standard deviation below 5%. We obtained SSE ranging from 1.492 mm · µs⁻¹ to 1.604 mm · µs⁻¹ for patients with MRI PDFF from 2% to 17%, and biopsy from 0% to 80% respectively. MRI PDFF and biopsy were used as gold standards with a particular attention to MRI PDFF, as it is to date considered as the most efficient approach for steatosis estimation (Idilman et al 2013). It is even considered better than liver biopsy as this latter is hindered by the non-negligible possibility to miss a steatosic region leading to false negative values (Sumida et al 2014). Moreover, as demonstrated by Leporq and colleagues (Leporq et al 2014), by giving access to the percentage of fat protons regarding to the percentage of water protons, MRI PDFF allow a linear steatosis staging, instead a groups classification like biopsy.

Our method is based on a three steps processing in order to guarantee accuracy. The first step provides an estimate of the sound speed optimizing the spatial coherence of backscaterred speckle noise. As this estimation can be degraded by aberrations occurring during propagation in the first layers, a virtual point like reflector is created in order to find and correct these aberrations. This aberration correction technique was found to be efficient as it allows in many cases (N  =  14) to increase the AUC of the spatial coherence function of the region of interest (a ratio increase of AUC_iteration1/AUC_iteration10  =  6%  ±  1%). Finally a third correction step is introduced in order to correct the bias introduced by the intercostal layer thickness for an optimal estimation of liver sound speed. This third step, tested on multilayer phantom experiments improve SSE robustness by including superficial layers thickness in the calculation.

These successive steps lead to a very high linear correlation (R²  =  0.691) between MRI PDFF and SSE and therefore enable to grade hepatic steatosis. Moreover, SSE was excellent to significantly differentiate healthy and diseased patients (p  <  0.0001, AUROC  =  0.942).

As a second confirmation, ROC curve analysis was performed to test the ability of SSE to differentiate healthy liver and steatosis with biopsy as gold standard (p  <  0.0001, AUROC  =  0.952). This analysis leads to the conclusion that SSE is able to diagnose healthy liver and steatosis as well as biopsy. Our technique can establish cutoff value between healthy patients and patients with steatosis (SSE_criterion  =  1.555 mm · µs⁻¹).

SSE strength is its ability to grade steatosis like MRI PDFF. Moreover, with both MRI PDFF and biopsy as gold standards, SSE was demonstrated to be valuable for hepatic steatosis 10% or greater detection, which is critical in patient who undergo living donor liver transplantation because of the risks of poor graft survival, initial graft dysfunction, and other complications (Trotter et al 2007, Sharma et al 2013).

Our ultrasound based method has the advantage to provide a more easy-to-use and accessible imaging modality for the screening of liver steatosis. Moreover, it could be combined with stiffness (Young's Modulus) measurements using Shear Wave Elastography (Deffieux et al 2015) or radiation force Shear Wave Speed estimation (Palmeri et al 2005).

Some authors see the interference caused by thoracic subcutaneous fat as the limiting factor for success in obese patients (Foucher et al 2006a). However as we included the thickness of the superficial layers and the phase aberration correction, our technique was found to be robust for all patients of this preliminary proof of concept study.

A limitation of this study is the population heterogeneity. Some patients had chronic liver diseases with various degree of fibrosis or even cirrhosis. As a preliminary study the role of fibrosis was not studied as a confounding factor, but another study with a larger patient cohort is planned to study fibrosis as a confounding factor for sound speed estimation by ultrasound.

This study proposes a robust method for sound speed estimation (SSE) in homogeneous and heterogeneous medium. This technique was first validated on two phantoms. The use of sound speed as a hepatic steatosis biomarker was then proposed and demonstrated in a preliminary clinical evaluation. In vivo sound speed was assessed using an ultrasound-based method combining speckle noise coherence analysis and an aberration correction technique. It was found to be highly correlated with MRI Proton Density Fat Fraction and biopsy results. This non-invasive technique was implemented on a research ultrasound scanner and could be combined with Shear Wave Elastography for global fibrosis and steatosis grading.

The authors kindly thank J Robin for proofreading the manuscript grammar. This work was supported by the Agence Nationale de la Recherche under the program 'Future Investments' with the reference ANR-10-EQPX-15 and by LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) and the French Program 'Investments for the Future' under reference ANR-10-IDEX-0001-02 PSL^*. This work was also partly supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC Advanced grant agreement 339244 FUSIMAGINE).