Quantitative Compression Elastography With an Uncalibrated Stress Sensor

Tissue stiffness is a key biomechanical property that can be exploited for diagnostic and therapeutic purposes. Tissue stiffness is typically measured quantitatively via shear wave elastography or qualitatively through compressive strain elastography. This work focuses on merging the two by implementing an uncalibrated stress sensor to allow for the calculation of Young’s modulus during compression elastography. Our results show that quantitative compression elastography is able to measure Young’s modulus values in gelatin and tissue samples that agree well with uniaxial compression testing.


INTRODUCTION
In modern medicine, diagnostic and treatment decisions are becoming increasingly influenced by measurements of biomechanical properties of tissues. This is especially true in the fields of oncology [1], dermatology [2], and hepatology [3]. Because disease states necessarily involve the physical change of the affected tissue, biomechanical properties are often a good indicator of both disease presence and extent.
A key biomechanical property that is most often measured in clinical use is tissue stiffness. Changes in tissue structure on the microscale result in changes in tissue stiffness on the macroscale. For example, tumors often vary in stiffness from the surrounding tissue. Subsequently, tumor stiffness can be used to detect tumors [4], aid in staging them [5], determine their metastatic potential [6], and aid in the evaluation of treatment [7]. Similarly, a liver disease resulting in fibrosis or steatosis can be diagnosed via liver stiffness measurement [3,8]. Treatment of many diseases in addition to the aforementioned diseases involving soft tissues, can benefit from the analysis of tissue stiffness for diagnostic and therapeutic planning and evaluation.
The field of elastography is generally focused on the measurement of tissue stiffness via various imaging modalities. For superficial or excised tissues, this can be accomplished with noninvasive optical techniques such as optical coherence elastography (OCE) [9][10][11], which can provide highresolution elastograms. However, due to light attenuation in tissue, optical-based elastography methods cannot be used in tissues that are deeper than a few mm beneath the surface unless an invasive probe is utilized [12]. Magnetic resonance elastography (MRE) does not have this limitation and is used to make whole-organ and whole-body elasticity maps but suffers from relatively poor resolution, typically in the mm range [13]. Additionally, it is expensive and requires a large amount of space. A third option is ultrasound elastography [14], which is portable, comparatively inexpensive, and has been practiced the longest and studied the most of the three techniques. It can offer both deep tissue penetration and high resolution, though there is a tradeoff between these two factors.
Ultrahigh-frequency ultrasound transducers with frequencies in the hundreds of MHz are capable of approaching optical imaging resolution both axially and laterally at the cost of penetration depth; however, these are still at an early stage and are currently manufactured in laboratories for use in experimental setups [15]. Elastography is typically performed at lower frequencies (< 5 MHz), resulting in resolutions of hundreds of microns up to mm scale [15].
The two main branches of ultrasound elastography are static/ quasi-static compressive strain elastography [16,17] and shear wave elastography [18,19]. In strain elastography, a series of images are taken as the tissue is slowly compressed, which results in a measure of how the strain distribution changes with compression. This method typically results in semiquantitative measurements and shows the ratio between the strains in different regions. Because of the manual compression, there can be high operator error and low repeatability [20], which is compounded by the nonlinear biomechanical properties of tissues [21]. Shear wave elastography, on the other hand, involves inducing a mechanical wave inside the tissue and measuring its speed. This allows for a quantitative measurement based upon the shear wave speed, i.e., shear modulus or Young's modulus. These mechanical parameters are absolute measurements of tissue elasticity and are more valuable because they allow for comparisons between different time points during the treatment of a single patient, between different patients, and between different diseases. Strain-based compressive elastography can also obtain absolute Young's modulus in tissues but requires the use of a precisely calibrated compliant sensor for measurements of surface stress, which has been utilized widely in optical coherence micro-elastography [22,23]. However, use of a stress sensor in ultrasound elastography is not as prominent as in optical methods.
In this work, we combine strain and shear wave elastography using a stress sensor placed on the sample to obtain twodimensional quantitative stiffness maps in a new technique called quantitative compression elastography (QCE). Unlike quantitative micro-elastography, our method does not require a calibrated sensor with precisely known mechanical characteristics. Instead, shear waves are generated within the sensor, which can be converted to quantitative elastic modulus values. Using these values and the strains measured within the sensor and the tissue of interest, we measure quantitative tissue stiffness values. We also compare these values to uniaxial compression testing (UT) for validation. To our knowledge, this is the first translation of strain imaging in ultrasound elastography to the quantitative mapping of biomechanical properties through the use of an uncalibrated stress sensor.

Experimental Setup
The experimental setup is shown in Figure 1. It consists of a Vantage 256 (Verasonics, Kirkland, WA, United States) ultrasound system with an L11-5 V transducer. The imaging and push frequencies were 7.8 MHz and the pulse duration was 128 μs. The transducer was placed into a custom holder to ensure stability and prevent movement. Each sample was placed on a linear z-axis stage which served to produce the compression. Ultrasound gel was degassed using a vacuum chamber and a thin layer was placed between the ultrasound transducer and the sample. The transducer was placed above the surface of the sample and guided into place by raising the z-axis stage.

Phantom Preparation
Gelatin phantoms were created by mixing gelatin (gel strength 300, type A. Sigma-Aldrich Corp, MO, United States) with distilled water and pouring it into a cylindrical mold with a diameter of approximately 8 cm. Silica powder was added to all phantoms and sensors to promote acoustic scattering and to make the boundaries easy to visualize on ultrasound imaging. Prior to measurement, all phantoms were removed from their molds via a custom lift mechanism to avoid tearing the phantom. This lift mechanism consists of a circular piece of thin acrylic cut with a laser cutter to match the diameter of the cylindrical mold. Two thin strips of metal were bent at a 90-degree angle and glued to the bottom of the acrylic to allow for the whole acrylic disk and the solidified phantom to be lifted as a single unit. When possible, the sensor was created by pouring the sensor gelatin mixture on top of the sample to avoid acoustic impedance mismatch that could create artifacts in the presence of an air gap or the insertion of bubbles. For each concentration, a separate small cylindrical phantom with a diameter of approximately 3.5 cm was also cast for validation via uniaxial mechanical compression testing.

Homogeneous Phantoms
Homogeneous phantoms with concentrations of approximately 10% gelatin (w/w) and 16% gelatin (w/w) were created. Silica powder at a concentration of approximately 0.8% (w/w) was added for scattering. They were poured into the aforementioned cylindrical mold allowed to solidify. A sensor layer of 9% gelatin (w/w) containing 0.5% silica powder (w/w) with a thickness of approximately 4 mm was poured on top and allowed to solidify.

Heterogeneous Phantoms
A half and half phantom was created to assess the capability of the technique to assess spatial variations in stiffness. Half of the phantom consisted of 16% gelatin (w/w) with 0.8% silica powder (w/w), and the other half was 10% gelatin (w/w) with 0.6% silica powder (w/w). The silica powder was varied to allow for visualization of the boundary between the halves. The gelatin was poured into the aforementioned cylindrical mold, allowed to solidify. A sensor layer of 9% gelatin (w/w) with 0.4% silica powder (w/w) with a thickness of approximately 4 mm was poured on top and allowed to solidify.

Tissue Preparation
Beef steak was acquired fresh from a local grocery store and cut into a small strip of approximately 5 cm by 10 cm. Half of the sample was cooked on a hot plate at 150°C for~1-2 min per side to change the local biomechanical properties of the region. Both sides of the cooked region were heated to ensure biomechanical changes were consistent through the entire depth, and the tissue was allowed to cool to room temperature before measurements.
An image of the tissue is shown in Figure 2, which clearly shows the delineation between the cooked and uncooked regions. When placed under the transducer, a copper wire was aligned with the visible cooked boundary for visualization on the ultrasound B-mode image before a 16% gelatin (w/w) sensor with a thickness of approximately 5 mm was placed on top. Prior to the acquisition, the wire was carefully removed from beneath the sensor.

Uniaxial Compression Testing
Phantoms of each concentration underwent uniaxial mechanical compression testing (Model 5,943, Instron Corp, Norwood, MA, United States) to measure Young's modulus by the "gold standard" for comparison. The height and radius of the cylindrical phantoms were measured with digital calipers, and the sample was coated with water to prevent friction between the compression testing frame plates. After coming into contact with the phantom, the test began and continued until either 20% strain or 50 N of force was reached. Compression was performed at 0.25 mm/s. The phantom was removed, and rewet after each measurement, and a total of five measurements were performed on each sample. The raw mechanical testing data was loaded into MATLAB (Mathworks, Inc., Natick, MA, United States), and stress and strain were calculated using the measured dimensions of the phantom. A polynomial best fit line to the stress-strain curve was obtained, and the derivative was calculated. This derivative equation corresponded to the tangent Young's modulus as a function of strain and was used for validation of the ultrasound compression elastography method.

Data Acquisition
Custom acquisition software was written in MATLAB to control the Vantage 256 system in conjunction with the z-stage. For each compression position, 120 plane-wave ultrasound images were taken and saved as IQ data. Following this, a shear wave was generated approximately 0.5 mm below the sensor layer, and an additional 120 plane-wave ultrasound images were taken to monitor the wave as it passed through the sensor and sample. An image was taken every 100 µs. The shear wave was generated by using 32 elements of the transducer, and the apodization and delays were calculated via software provided by the ultrasound manufacturer. This ensured that the wavefronts from each transducer element arrived at the shear wave excitation focal point at the appropriate time to generate a shear wave. The z-stage was raised by 100 μm, and the process was repeated until the desired number of compressions was reached.

Data Processing
A flow chart of the processing steps is shown in Figure 3. The shear wave data sets were processed as follows. First, the particle velocity was obtained in each of the 120 frames using the vector method [24] with a kernel size of 20 pixels, which corresponded to 1.35 µm. The particle velocity, i.e., the temporal derivative of the displacement in the axial direction, for each pixel was summed in time to obtain displacement in each frame. Then, the displacement images were smoothed using a Savitzky-Golay filter with polynomial order of one and a frame length of 21. The data were interpolated in time by a factor of 5. The Fourier transform of the displacement data in time was taken at the level of the sensor in the path of wave propagation to determine the frequencies corresponding to the shear wave. The displacement data stack was then bandpass filtered to isolate the frequencies corresponding to the shear wave. A cross correlation-based method was utilized to obtain the temporal lags of the wave propagation in the sensor. A shear wave speed map was generated from the wave propagation time lags and spatial coordinates using a sliding window with a kernel size of approximately 1.69 mm.
The compression data was processed as follows. The IQ data for each compression data set was loaded, and the 120 frames were averaged in time to remove noise for each compression step. The result was a stack containing complex B-scan images of the sensor and sample. The complex images were converted to intensity grayscale images, which were then segmented using MATLAB's volumeSegmenter tool. The vector method was used on the complex image stack to generate particle velocity maps, which were then cumulatively summed in time to generate displacement maps. Following this, each displacement image was normalized such that the displacement at the interface between the sensor and the transducer was zero. The strain was calculated from the displacement frame by using the leastsquares regression fit of the axial displacement with a window size of approximately 0.225 mm [25].
The Young's modulus of the sensor was calculated from the speed of the shear wave. The area of the excitation pulse was removed to avoid near-field effects, which overestimates the velocity of the waves generated in or near the shear wave excitation focal point. The wave speed map of the sensor was converted to Young's modulus, E, via the following relationship [10]: where ρ was the mass density of the sensor material and c s was the shear wave speed. The mass density was assumed to be 1,000 kg/ m 3 . The speed in the sensor was averaged to obtain Young's modulus.
The stress in the sensor was then calculated via the following relationship: σ Eε where σ was the stress in the sensor and ε was the strain in the sensor. Assuming that the stress was uniform axially between the two rigid boundaries of the setup, i.e., the transducer and lower plate, Young's modulus was computed for each pixel in the strain map of the sample by The resulting images were then smoothed via cubic spline interpolation. The average of the segmented area was then taken after outliers were removed to obtain the final average Young's modulus value for the sample or region of interest.

Gelatin Phantoms
The results from the homogeneous phantoms are shown in Figure 4 and in Table 1. For the 10% phantom, the average value was 22.36 ± 1.59 kPa, and the average for the 16% phantom was 41.39 ± 3.15 kPa. The slope of the stress-strain curve from uniaxial testing resulted in Young's modulus of 22.73 kPa for the 10% phantom and 40.14 kPa for the 16% phantom. The average percent errors for the 10% phantom and 16% phantom were 5.91% and 6.68%, respectively. The results from the half and half gelatin phantom are shown in Figure 5. For the softer side, the average value of the measurements using the ultrasound compression method was 23.33 ± 1.66 kPa. The slope of the stress-strain curve from uniaxial testing resulted in Young's modulus of 23.01 kPa. Overall, there was an average percent error of 5.60% by QCE as compared to mechanical testing. For the stiffer side, the average elasticity was 40.83 ± 1.75 kPa as measured by compression elastography. The slope of the stress-strain curve from mechanical testing was 41.44 kPa. Overall, there was an average percent error of 3.98% by QCE as compared to mechanical testing.

Tissue Sample
The measured elasticity values for the tissue sample can be seen in Figure 6 and in Table 2. The elasticity of the uncooked steak increased from 16.25 kPa to 31.84 kPa over a compressive range of 4.94%-8.87%, as estimated by QCE.  For the strain ranges shown, the cooked half of the steak had a much more linear slope with increasing compression than the uncooked half of the steak. The solid blue line is the derivative used to approximate the stress-strain curve of the uniaxial testing data. The raw tissue did not have a linear trend of tangent Young's modulus beyond~7%. The mechanical testing uniaxial data is plotted in Figure 6 as the solid red and blue lines for the uncooked and cooked steak, respectively.

DISCUSSION
Our results show good agreement between Young's modulus measured with QCE and Young's modulus measured by uniaxial testing. These results validate QCE as a simple, multimodal, accessible, and relatively low-tech method to obtain quantitative information on tissue stiffness.
Gelatin was chosen because, for small strains, Young's modulus is approximately constant as a function of strain [26,27]. Some studies have shown an approximately linear stressstrain curve up to approximately 20% strain [28]. We found that this remained a good approximation in the strain region we considered, which was below 10% strain. When measuring the homogeneous phantoms, the errors remained below 7% when compared to uniaxial compression testing. When estimating the elasticity of each side of the heterogeneous phantom using QCE, the errors remained below 6% as compared to uniaxial mechanical compression testing, indicating that QCE is capable of a high degree of accuracy for tissue stiffness assessment.
Hooke's law is valid only in the case of materials with linearelastic behavior; i.e., in situations where the stress-strain curve is linear. Here, we assume that we are within the linear-elastic region, that the strain increases slowly, and that there is no permanent modification to the material.
The tangent Young's modulus values measured in the tissue sample (steak) via QCE were in good agreement with mechanical testing as well. There was a significant preload on the raw side prior to elastography measurements being taken. This is because the tissue shrank during cooking, and thus, the initial strain on the thicker uncooked side was greater as  compared to the thinner, cooked side. Though the raw and cooked tissues were under different loads during the experiment, accurate results were still obtained without the use of a calibrated stress sensor. As in the gelatin phantom, the error remained low at~6%, proving further that QCE is able to accurately measure tissue stiffness. The tissue did not exhibit linear elasticity, which is as expected [21,29,30]. While most tissues show strain hardening, this behavior varies by tissue type [31,32]. This is characterized by a large initial compression with low levels of stress followed by strain hardening, where the stress increases in a nonlinear fashion with increasing strain [31,33]. The raw tissue exhibited this to a greater degree than the cooked tissue, but QCE can accurately measure tissue stiffness regardless of whether it is in the linear or nonlinear regime since the stress measurements are made for a given strain.
Typically, to perform quantitative strain elastography, a compliant sensor with precisely known properties is needed [22,23]. This involves uniaxial testing beforehand to measure these properties and precise measurements of the pre-stress before imaging. Uniaxial testing machines are expensive and are unlikely to be found in most hospitals or clinical settings, unlike ultrasound imaging systems. Additionally, if the sensor breaks or needs to be changed, calibration must again be performed on the sensor, which is particularly noteworthy in cases where sterility is necessary. One major novelty of the presented method is that it eliminates the need for the calibration step. The sensor need only be linear in the strain region that will be utilized in the sensor and sufficiently scattering to visualize shear wave propagation. If the sensor enters the nonlinear region, the elasticity is no longer accurately estimated by Eqs 4-1, which means that Young's modulus value will not be correct. An alternative form of the relationship in Eq. 3 shows the problem more clearly: If E sensor is nonlinear with compression, then this nonlinearity will be transferred into E sample and give an incorrect measure of elasticity despite having the proper strain ratio.
An additional constraint on the sensor is that it must not be too soft or stiff. If the sensor is too soft, then during compression, most of the deformation will go into the sensor instead of being distributed between the sensor and sample, and there will be no measurable strain in the sample. If the sensor is too stiff, it will not deform during compression, and there will be no measurable strain and, thus, stress. In both cases, the method will likely fail due to the lack of a measurable strain. The strain must be accurately measured in both the sensor and the sample for this method to be successful. In the clinic, acoustic coupling gel is often used, but acoustic coupling gel pads are also commercially available, which may serve as stress sensors. One such gel pad is Aquaflex, which has been studied previously [34] and is widely available. However, their mechanical properties are not often known, and thus, this technique would be useful in such use cases assuming the stiffness of the gel pad is not too low or high. One difficulty in using commercial gel pads is that they are designed to be acoustically transparent, which would limit visualization of the shear wave propagation and, therefore, the estimation of shear wave speed. A commercially manufactured gel pad with a low concentration of acoustic scatterers would be preferable for this application.
Clinical ultrasound imaging machines typically do not allow the transducer to act as a shear wave source, which means that this method cannot be used in those settings. Though it is most convenient when the transducer can act as both a shear wave generator and an imaging modality, this method could be performed with an additional wave source, such as a mechanical actuator utilized in MRE [35]. The main concern 2 | Values of stress and strain for the cooked and uncooked sides of the steak experiment.

Steak
Step Cooked strain (%) Cooked stress (kPa) Raw strain (%) Raw stress (kPa) is that wave propagation with a sufficient signal-to-noise ratio [36] must be observed in the sensor. Shear wave propagation in the sample is not necessary when using this method. This means that the method can work for any deformable tissue even in the presence of high dispersion or attenuation, provided that strain in the tissue can be accurately measured. If the underlying structure is not deformed enough, then strain cannot be determined, and therefore the value of E sample is unreliable. In this case, it is unlikely that either strain or shear wave elastography would be possible, e.g., in bone. However, having shear wave propagation in the sample is advantageous, as it allows for additional analysis and comparison, such as for multimodal elastography [37,38], which has shown improvements in disease detection as compared to ultrasound imaging alone.
Two key processing steps have the potential to change the accuracy of QCE. First is the estimation of the shear wave speed in the sensor. If the speed in the sensor is poorly estimated due to poor wave propagation, poor scattering, or other factors limiting wave propagation or its detection, then the accuracy of this method will decrease. The algorithm with which speed is estimated also determines accuracy to some degree. If using windowed cross-correlation, for example, the choice of window size will influence the accuracy [39]. The second key step is the estimation of strain within both the sensor and the sample. The algorithm used here will determine the key factors that must be considered. Like the estimation of speed, the window across which the strain is estimated will also influence the accuracy [40,41].
Studies in optical coherence elastography have shown that quantitative strain elastography has better contrast and mechanical spatial resolution than wave-based elastography methods [42,43]. These methods rely upon high SNR because they depend upon being able to track small movements either in via phase-sensitive methods or via speckle tracking [40]. Additionally, the use of a sensor decouples the measurement from the optical (in the case of OCE) or acoustic properties of the tissue, meaning that stiffness can be estimated in samples with poor imaging visibility because the underlying mechanical properties can be measured by their effects on the sensor. This may also mean the ability to detect the effects of structures that lay beyond the imaging depth [22].
There are several limitations to this technique. First, this technique is only capable of capturing a 2D representation of tissue stiffness. These experiments were performed in a very controlled manner by using the custom holder. Ideally, this technique would be performed clinically by manually controlling the transducer on patients where discrete 100 µm steps would not be possible or practical and where significant outof-plane motion is expected. It may be possible to overcome this problem with transducers that provide 3D information or by incorporating a fixed transducer or sample holder. Additionally, the effect of compression step size on decorrelation between image frames was not considered. This method relies on the correlation between successive stepped frames, and too much movement within the frame will degrade this correlation. Future work will focus on applying QCE to various tissue types, investigating the relationship between window size and accuracy, assessing the repeatability of QCE across different tissue samples, examining the effect of the rate of loading on measurements, and evaluating the technique on highly anisotropic and highly viscoelastic tissues.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JR created the setup, designed and manufactured the custom holder, performed the experiments and analysis, and wrote the paper, MS and AN provided expertise throughout and helped guide processing methodology. MS and AN also helped with strain estimation algorithm implementation. SA and KL provided technical guidance, supervised data analysis, and initiated the research. JR, MS, SA, and KL came up with the concept. All authors contributed to manuscript revision, read, and approved the submitted version.