Viscoelastic Parameter Estimation Using Simulated Shear Wave Motion and Convolutional Neural Networks

I. Introduction

Ultrasound imaging is one of the main imaging modalities used in healthcare. It is widely used for multiple applications due to its low cost, relative portability and safety [1]. One of the applications that has developed popularity over the last several years is ultrasound shear wave elastography (SWE). SWE typically uses focused ultrasound beams to produce acoustic radiation force (ARF) that pushes within the tissue [2], [3]. When the ARF ceases, the locally perturbed tissue returns to its equilibrium position, but also generates waves traveling perpendicular to the push direction. These shear waves propagate through the tissue, and the motion is measured using high frame rate ultrasound techniques [4], [5]. The wave motion can be analyzed to extract information about the wave velocity and attenuation. These parameters can also be used to characterize rheological properties of the tissue such as elasticity and viscosity.

The estimation of rheological parameters is dependent on mathematical models, such as the Kelvin-Voigt and Maxwell models [6], [7]. To date, methods have been developed to directly measure the wave velocity and attenuation of shear waves as they vary with frequency [7]-[15]. However, the shear wave attenuation can be difficult to directly measure [11]. With the wave velocity dispersion information alone, it is possible to estimate viscoelastic properties for different rheological models [6]-[9], [16], [17]. Most of these methods involve Fourier analysis to extract dispersion curves that are used for fitting to the rheological models. Each of the approaches that have been studied have different levels of bias and variability. Additionally, different aspects of the SWE acquisition such as the ARF distribution and beam shape and measurement noise can affect the wave motion and confound these methods [18]-[20]. Taken together, there is still a need for robust methods for viscoelastic parameter evaluation from measured shear wave motion.

Machine learning (ML) is the computer science field dedicated to creating algorithms capable of modeling complex datasets, and performs tasks such as prediction, regression or clustering. These methods commonly have high computational costs. In recent years, there has been substantial progress on both algorithms and computational power to perform data training and predictions in feasible time with available computational resources. ML can identify patterns within datasets that otherwise would be very difficult to fit by standard statistical models or by evaluating each feature on its own. The use of ML might enable SWE to reach higher applicability, by naturally modeling the intrinsic complexity of viscoelastic tissues [21]-[24].

Artificial Neural Networks (ANN), or simply Neural Networks (NN), were inspired by the central nervous system structure and its physiology. They are composed by nodes, or neurons, interconnected into multiple layers, from input to output. Each node has adjusted weights and activation functions that simulate the synaptic reinforcement plasticity. The adjustment of the weights is implemented by solving the convex problem that minimizes the loss of the networks, which evaluates the prediction assigned by the network against the ground truth output [25], in a process called supervised learning.

The convex problem is solved using numerical approaches, such as gradient descent, which are iterative, and therefore, can take time and large amounts of data. The process of weight adjustment is called optimization and it is an area of research of its own. Some of the most common optimizers used for neural networks are Stochastic Gradient Descent [26] and Adam [27].

Convolutional Neural Networks (CNN), as other NN architectures, were also inspired by neurophysiology. The architecture replicates the layer-based structure from the human visual system, where the convolutional layers work as the receptive fields present in the brain, while the max-pooling layers act as the lateral geniculate nucleus by reducing the level of information and extracting the relevant features for the task in hand [28]. The input is convolved with the feature weights to produce a feature map which is processed by a non-linear activation function (rectified linear units, traditionally). The pooling layers then reduce the spatial resolution of the feature maps to achieve spatial invariance to input distortions and translations. Max-pooling layers aggregate and propagate the maximum value in a receptive field to the next layers. After stacking multiple convolutional and pooling layers, a fully connected layer is implemented to interpret the feature representations extracted by previous layers, much like the human visual cortex. The fully connected layer is a traditional NN, with different output layer activation functions and loss functions, depending on the task [28]. The CNN basic schema is exemplified in Fig. 1. CNNs have been extensively used for image recognition tasks with high performance and generalizability [22], [28].

In this proof-of-concept study, we propose to analyze, in a controlled setting, the capabilities of ML to model wave propagation and retrieve elasticity and viscosity values from wave motion images simulated using the staggered-grid finite-difference (SGFD) approach to implement a Kelvin-Voigt (KV) material model [29]-[31]. Due to the two-dimensional (2D) image characteristics of the task, the CNN architecture was deemed the best for the objectives of this proof-of-concept. We also compared the ML performance to the two-dimensional Fourier transform (2D FT) analysis approach for evaluating dispersion curves and fitting to the KV theoretical dispersion curve [6], [7], [16].

The modeling of viscoelastic wave propagation is useful for understanding wave phenomena in complex media. Finite-difference modeling approaches have great importance for geophysical exploration, reservoir engineering, and military applications because they provide complete wavefield responses at reasonable computational cost [30], [32]. The SGFD method was used to simulate wave propagation in viscoelastic homogeneous media. SGFD adopts a finite difference operator with fixed-order accuracy to calculate spatial derivatives for a homogeneous or heterogeneous medium. The method’s accuracy, efficiency and stability are improved by the use of staggered grid for spatial discretization [30], [31].

II. Methods

For this study, the simulations were created using SGFD simulations implemented in MATLAB (Mathworks, Natick, MA, USA) and computed on GTX 1080 and RTX 2080Ti GPU cards (Nvidia, Santa Clara, CA, USA). Each simulation took approximately 1-8 hours of computing time depending on the GPU card used. The SGFD simulations implement a KV material in the Navier-Stokes equations to calculation the motion caused by acoustic radiation force from ultrasound push beams. The shear elasticity and viscosity ranges for the KV model were chosen so they would cover a wide range of values for normal and pathological soft tissues [33]. The shear elasticity, μ₁, ranged from 1-25 kPa with increments of 1 kPa, and shear viscosity, μ₂, ranged from 0-10 Pa·s with increments of 0.5 Pa·s. The 25 values of shear elasticity and 21 values of shear viscosity yield 525 unique combinations of viscoelastic parameters. The medium is assumed to be homogeneous, isotropic, linear, and incompressible [29].

We used two different ARF excitation focus profiles with f-numbers (F/N = focal length/width of aperture) 1 and 2 (Fig. 2) with a focal distance of z_f = 2.5 cm, adding up to 1050 simulations. These values of F/N are used frequently in practice and provided cases of a very focused beam (F/N = 1) and a less focused beam (F/N = 2). The ARF was simulated using FOCUS using the model of an L7-4 linear array transducer (Philips Healthcare, Andover, MA) [33], [34]. The element pitch was set to 0.0308 cm, the element kerf was 0.0025 cm, the element height was 0.7 cm, the elevation focus was 2.5 cm, and the frequency was 4 MHz. The medium was modeled using a density of 1000 kg/m³, speed of sound of 1540 m/s, and ultrasound attenuation of 0.5 dB/cm/MHz. A focal depth of 2.5 cm was used and for F/N = 1 and 2, such that 82 and 40 elements were used, respectively. The pressure was simulated, and the intensity was calculated as it is proportional to the ARF in absorbing soft tissues. For the models, the force was applied for 200 μs. For each combination of μ₁, μ₂, and F/N, a dataset was produced.

An important factor of finite-difference methods is the model’s boundary conditions. The boundary conditions have to absorb and attenuate the waves without creating artifacts such as reflections. It is also important that boundary conditions do not overly increase computational costs. Perfectly matched layers (PML) are considered the optimal absorbing boundary condition for finite-difference methods, and therefore, were implemented in this study [35].

Another important aspect of SGFD simulations is the Courant-Friederichs-Lewy criterion (CFL), as it describes the conditional stability that allows convergent behavior of the solution. The space-time discretization is dictated by the smallest shear elasticity (μ₁) and the highest shear viscosity (μ₂) to be simulated, and so, the CFL criterion determines the time increment and thus the number of time steps to achieve a certain simulation length for a given shear wave velocity.

where the dimensionless number C is called the Courant number, V_p and V_s are the compressional and shear wave velocities, Δt is the time step, and Δx the spatial step in the x dimension. In this study’s simulations, the Δx and Δz (spatial step in the z dimension) were set to 0.1 mm, in a field of 44 mm width and 32 mm height. To optimize the simulation run time, two different values of C were used in the simulations depending on the media properties: 0.2 if μ₁ > 5 kPa and 0.05 if μ₁ ≤ 5 kPa and μ₂ ≥ 4 Pa·s. The simulation Δt was then calculated based on Eq. (1), for a total period of 25 ms. To meet the conditions described above, the temporal resolution was very high (5-7 MHz). After completion, the simulated data was resampled in time to 16 kHz, to ensure that all models would have the same time discretization. This sampling frequency could be reached with interpolated data that would be upsampled for precise time-of-flight estimates for evaluation of shear elasticity estimations [17].

The three-dimensional simulations (z, x, t) were then cropped and sampled in the z-direction in order to create the 2D x-axis vs. time wave motion dataset for training. The excitation generates two identical wave fronts propagating perpendicular to the push, therefore, the simulation region-of-interest (ROI) used to create the dataset was composed solely by the right hand-side wave front (Fig. 3a). The ROI has a width of 2 cm with its origin at 0.2 cm laterally from the push center, to eliminate the high amplitude motion of the push from the dataset images. The time range adopted was 20 ms. The data that was used for the CNN input was the motion from particular depths in the z-direction. To augment the dataset the right hand-side lateral propagation was used from five different vertical positions relative to the push focus (Fig. 3a), at −0.90, −0.45, 0.00 (center), 0.45 and 0.90 cm (Fig. 3). Another technique implemented to augment the number of samples, was adding different levels of Gaussian noise to the lateral propagation images to generate motion with signal-to-noise ratios (SNRs) consistent to literature, at 20, 10 and 5 dB (Figs. 3b and 3c) [29], [36], [37]. The segmentation and augmentation of the dataset was also implemented in MATLAB. With the five depth positions and three noise levels and 1050 unique combinations of the viscoelastic parameters, the total number of images used for the study was 15,750.

The CNN architecture was implemented in Python (Python Software Foundation, Wilmington, DE, USA), using Keras [38] with Tensorflow (Google Brain, Mountain View, CA, USA) backend, and trained on a Z820 workstation (Hewlett-Packard, Palo Alto, CA, USA) equipped with a GTX 1080 GPU (Nvidia, Santa Clara, CA, USA). The data in each image was normalized by subtracting the minimum amplitude from the whole ROI and dividing all motion amplitude values by the maximum amplitude after subtraction, creating images where all the values had an amplitude between 0 and 1. The labels were also normalized with a normalizing factor of the maximum label value (μ₁ = 25 kPa for elasticity and μ₂ = 10 Pa·s for viscosity) plus a 20% margin, to leave a margin for model overestimation in the top of the range. The same margin was not used at the bottom of the range (μ₁ = 1 kPa for elasticity and μ₂ = 0 Pa·s for viscosity), because negative values are not physically possible.

The training, validation and test sets were defined using a group shuffle split of 60/20/20, which means that images derived from the same simulation (μ₁, μ₂, F/N combination) could not be present in more than one set simultaneously. This ensures that simulation subjects presented to the CNN in the validation and test datasets were never seen during training in any of their augmented versions. Due to the regression characteristics of the task, sigmoid and mean squared error (MSE) were implemented as output layer activation function and loss function, respectively.

The elasticity and viscosity determine different characteristics of the wave propagation. The elasticity is related to the time-of-flight (e.g., slope of the wave motion distribution in the x-t plane, Fig. 3b), whereas the viscosity affects the widening of the wave at a particular location.

Therefore, for each dataset, a pair of CNN models were implemented, one for elasticity regression, and one for viscosity regression. To facilitate the compilation of results, a coding system was applied to name the models and results in terms of what the training/validation and test f-numbers were used (Table I). The letter designates what label was used, elasticity (E) or viscosity (V); the first number specifies the F/N used for training and validation, whereas the second number designates the F/N used for testing.

The CNN models trained, validated and tested on data from both f-numbers were denoted by the X1+2 designation. The CNN models were optimized for multiple combinations of 1, 2 and 3 convolutional layers; 2, 3 and 4 dense layers; 16 and 32 nodes per layer (convolutional and dense); 0.001, 0.0005 and 0.0001 learning rates; totaling 54 models. Tensorboard (Google Brain, Mountain View, CA, USA) was used to monitor and evaluate the models during the hyperparameter grid search. By evaluating the validation loss and mean absolute error (MAE) during 200 epochs, it is possible to determine when the learning reached an asymptote or demonstrated overfitting. MAE is defined by the equation:

where N is the number of test subjects, ${\widehat{x}}_i$ is the estimated value of the subject and x_i is the true value of the parameter. Each 200-epoch training session took approximately 1 hour of GPU computing time on the GTX 1080 GPU card. For each grid search, the top 20% CNN models with the lowest validation MAE were selected and re-trained for a subsequent round of evaluation until one hyperparameter set remained (Fig. 4). This approach was used to mitigate hyperparameter overfitting, because each iteration of cross evaluation is performed on different random group splits. The final model was then trained, also for 200 epochs, using the optimal hyperparameters. Checkpoint recall was used to store the best model generated and prevent overfitting. The models were then evaluated against the test set.

The datasets were also evaluated using 2D FT analysis and the estimated dispersion curves were fit to the theoretical Kelvin-Voigt model to obtain the viscoelastic parameters. This analytical approach was implemented in MATLAB. The 2D fast Fourier transform was performed with a zero padding factor of 4096 for optimum frequency and wavenumber resolution of the frequency domain representation (k-space). The dispersion curves were than estimated from the k-space using the relationship c_p(ω_p) = ω_p/k_max, where c_p is the phase velocity at ω_p frequency from 100-600 Hz, and k_max is the wavenumber of maximum k-space energy at ω_p. Additionally, c_p values below 0 m/s and above 10 m/s were discarded from further analysis. The dispersion curves were then fit to the theoretical Kelvin-Voigt model (Eq. 3), with the density ρ = 1000 kg/m³, to obtain μ₁ and μ₂ values. The estimated values were then compared to the CNN results.

III. Results

The testing results were evaluated using the MAE between the true value of the viscoelastic property and the CNN prediction output. A linear regression was performed to compare the results from the 2D FT analysis as well as the CNN prediction to the identity line, and the slope, intercept and R² of the regression line were also calculated and reported. The optimal hyperparameters were also noted and reported in Tables II and III for the elasticity and viscosity CNN models, respectively. All elasticity models trained and tested on the same f-numbers showed very good fitting with MAE below 0.079 kPa with slopes and intercept of approximately one and zero, respectively. The viscosity models trained and tested on the same f-numbers also showed very good agreement with MAE below 0.091 Pa·s with slopes and intercept of approximately one and zero, respectively. Alternatively, CNN models tested on f-numbers different from the training/validation F/N had diminished performance with MAE of at least 1.542 Pa·s (V12) and as high as 2.834 kPa for E12.

The errors for each CNN model were also aggregated and plotted versus the true value (Fig. 5). The data were reported with boxplots where the red line represents the median, and the edges of the box represent the interquartile range (IQR, 25^th and 75^th quantiles). The whiskers represent 1.5*IQR above or below the 75^th or 25^th quantile, respectively. Values denoted by markers outside of whiskers are denoted as outliers. For cases E11, E22, and E1+2, the level of bias was minimal over the range of true shear elasticity. For E12, μ₁ was overestimated when the true μ₁ < 15 kPa, and underestimated for larger values of μ₁. For E21, a consistent underestimation of μ₁ was observed. In the cases of V11, V22, and V1+2, the bias was relatively low, where the MAE was less than 0.069 Pa·s. For V12, μ₂ was consistently underestimated and for V21 μ₂ was consistently overestimated.

When compared to the 2D FT analysis, the CNN results were found to be significantly more reliable. The Fourier analysis showed reduced performance at higher viscosity and lower elasticity combinations, driving elasticity underestimation and viscosity overestimation. In these types of materials, the ARF excitation was not able to generate energy at higher frequencies (300-600 Hz), and therefore it limited the quality of Kelvin-Voigt fitting. Although the linear regression technique was able to retrieve acceptable values for slope and intercept for the elasticity estimations, the errors still drove MAE values above 2.2 kPa, and R² below 0.78. The viscosity estimations were more inferior to the performance achieved with the CNNs, with slopes above 1.5, MAE above 2.9 Pa·s and R² below 0.7.

The errors were also plotted versus vertical position in Fig. 6. Similar trends were observed for these additional parameters as compared to those discussed above. The prediction of μ₁ was largely insensitive to the vertical position of the wave motion used, though the results for model E1+2 exhibited more outliers. For E12, the bias is positive (average of 0.14 kPa), and 1.5*IQR of over 0.63 kPa showing a very large variability, but insensitive versus vertical position. For E21, there was a larger negative bias for data from the pre-focal locations (−0.27 kPa at −0.90 cm and −0.1 kPa at the focus), while larger numbers of outliers for post-focal datasets. For V11, V12, and V1+2, the bias was small, but the number of outliers slightly increased for the V1+2 case. The errors in μ₂ for V12 were largely insensitive to the position, though for V21 the error increased both in bias and variability as the data farther from the focus was used, with 0.1 ± 0.18 kPa at the focus, 0.15 ± 0.22 Pa·s at −0.90 cm and 0.28 ± 21 Pa·s at 0.90 cm. When compared to the 2D FT analysis, the CNN is more reliable with respect to the depth position, when trained with the proper f-number. The same sensitivity to off-focus position was also observed using the 2D FT analysis for F/N = 1 (Fig. 6 – E1-2D FT and V1-2D FT).

The errors for the models as a function of noise level are shown in Fig. 7. The model results for the elasticity prediction were not largely affected by the level of noise. The cases of E1+2 and E21 did demonstrate a larger number of outliers. The models trained/validated and tested with different f-numbers showed bias invariance to the noise levels, nevertheless, E11 showed an increase of 21% and 82% in 1.5*IQR from 20 dB to 10 dB and 10 dB to 5 dB, respectively, while V11 showed an increase of 60% and 38% in 1.5*IQR from 20 dB to 10 dB and 10 dB to 5 dB, respectively. The results for model V1+2 showed more outliers than V11 and V22 models. The 2D FT analysis showed a larger number of outliers, when compared with the CNN. While the number of outliers were relatively constant, the 2D FT results were similar across different levels of noise.

For models trained, validated and tested with both f-numbers, the errors were also compiled by f-number (Fig. 8). No relevant difference was found between the errors for both elasticity and viscosity. Because all the data is aggregated, there are a high number of outliers present, but the error has a similar overall range for all models and F/N values (Tables II and III). The outliers cover a similar range for both F/N values for V1+2. Models trained and tested on the same f-numbers generally showed very low absolute errors with no relevant bias and standard deviation of less than 0.5 kPa or 0.5 Pa·s throughout the true value range.

Models E21 and V12 showed a pattern of higher bias and IQR with the increase of true value, whereas E12 and V21 showed more complex patterns of error across the true value range. In some cases, even with the overestimation margin applied to the label normalization, the higher value estimations were saturated and drove the bias negatively, e.g. E11 and V11. It is also notable that model V1+2 presented a higher level of outliers and positive bias at 0 Pa·s.

IV. Discussion

The results shown are encouraging and exhibit the potential of the use of convolutional neural networks in ultrasound SWE for viscoelastic characterization. The results demonstrated that the CNN models were robust to noise, vertical position and partially to F/N. When tested on datasets with the same excitation F/N for training and testing, the elasticity and viscosity estimations were very good and showed robustness to noise and vertical position.

The CNN models were marginally reliable when presented with test data from a different F/N than used for the training, with a trend of overestimation for E12 and V21, and underestimation for E21 and V12. Vertical position showed some level of variability in models trained with data using F/N = 2 and tested on data with F/N = 1 (X21), possibly due to the planar characteristics of higher F/Ns, when compared to a lower F/N, which produced more curved wave fronts.

Higher f-numbers have less focused beams that might cause wider lateral shear wave profiles that could simulate the widening effect of viscosity in the wave propagation. When trained on data using F/N = 1 and tested on data when F/N = 2 (X12), the results showed a consistent underestimation profile increasing linearly with the true viscosity (V12, Fig. 5), while the opposite effect occurred with V21, showing consistent viscosity overestimation. The V21 overestimation is exaggerated with the distance from the focus, as F/N = 1 has a more curved intensity profile to produce the wave front while F/N = 2 has a more planar wave front (Fig. 2). This behavior is well exemplified in Figs. 6 and 9. This characteristic is important because it suggests the possibility of reliable estimation correction for cases of F/N for which the model was not trained.

When different F/N values were used, the elasticity estimations showed the opposite behavior from the viscosity, overestimation for E12 and underestimation for E21. It is possible that the more curved intensity profile in the z-direction used for the ARF push beam, and the presence of more pre-focal and post-focal intensity for F/N = 1 also leads to errors in elasticity estimation, as it might change the perceived slope of the wave motion (Fig. 10).

Similarly to viscosity, the elasticity showed correlation between absolute bias and the distance to push focus when trained with F/N = 2 and tested on F/N = 1 (X21). Models trained with F/N = 2 are tuned to the planar wave front characteristics, therefore, when presented with data from a lower F/N for testing, show higher levels of error off-center (Fig. 6, V21 and E21).

Noise did not affect the estimations significantly in any of the scenarios tested, which suggests the robustness of the technique to SNR changes.

When trained with both f-numbers the CNNs were able to overcome the shortcomings discussed above and performed similarly as the pairs trained, validated, and tested with a single F/N. The results indicate that the CNNs can discern between wave front profiles with high accuracy, as long as they are trained properly.

Overall, this proof-of-concept study showed that the CNN models were able to successfully retrieve the viscoelastic properties from the wave motion images. The approach proposed required minimal preprocessing, without the employment of filters or transforms. The conventional 2D FT analysis approach showed inferior performance on the same datasets, with worse MAE and R² levels found compared to those found with CNN regression, even in cases where the CNN was trained and tested with different f-numbers. The 2D FT analysis was especially affected by materials where the effects of the viscosity was more prominent compared to the elasticity (e.g., high μ₂ and low μ₁), whereas the CNN did not show the same behavior and was reliable over the full viscoelastic parameter range. As noted in the Results, in the cases where the viscosity effects were substantial, energy at higher frequencies was diminished and curve fitting was more error prone. These findings with the CNN models are particularly encouraging, because changes in viscoelastic properties are related to pathological changes to the tissue, e.g. fibrotic tissue is correlated with higher μ₁ [10], [16]. Therefore, patients and practitioners can benefit from ML classification of tissue state based on SWE, which is an affordable and non-invasive measurement.

Although encouraging, the methods discussed still need further investigation. The method needs to be extended to allow for heterogeneous cases and other noise characteristics. The practical use of ML in ultrasound SWE has specific challenges such as the noise due to ultrasound detection and physiological motion, the standardization of images, and simulation model limitations that will have to be addressed in future work.

In this study, we used simulated motion from the SGFD based on a Kelvin-Voigt rheological model. Although it has been shown in literature, in multiple works, that the Kelvin-Voigt model does describe shear wave velocity dispersion over the typical elastography analysis ranges of frequency, such as 100-600 Hz [39]-[42], the use of the KV model on SGFD simulations might impair the use of simulations for ex vivo and in vivo transfer learning as not all tissues may be well-described by a KV rheological model. We also chose to use the 2D x-t wave motion images to directly estimate the values of μ₁ and μ₂. Other ML approaches or NN architectures could be utilized to estimate the dispersion curves for fitting to different rheological models.

We did not incorporate the ultrasound detection in these simulated data, which may add additional bias, but this proof-of-concept approach shows the promise of using shear wave motion and the overall shape of the wave fronts to estimate the viscoelastic properties. We did explore adding Gaussian noise to the wave motion data, and the models were largely robust to different noise levels. Higher levels of noise could be explored in future work. Physiological motion will also introduce noise as well that is not accounted for in this work. Each of the images used for these models used a fixed range of lateral space and measurement time. For in vivo implementation, this would have to be standardized for a given application or transducer. We also limited the range of depths over which we made the evaluations. Typically, viscoelastic measurements are made within or near the depth-of-focus of the ARF push beams, but an evaluation of parameter estimation versus depth will be investigated in our ongoing work.

Ultrasound SWE implementations have a wide variety of characteristics, such as F/N, transducer profile, motion detection approaches, all of which with specific use cases, strengths and weaknesses, and therefore, need to be evaluated separately. Ideally, for in vivo implementation, an acquisition approach with the aforementioned parameters would need to be thoroughly modeled and likely used in the training and validation and later testing for a robust technique. The entire pipeline including ARF push beam modeling, wave motion modeling with different material models, ultrasound detection, and motion estimation and filtering would need to be incorporated to all for transfer learning with experimental data from phantoms or ex vivo and in vivo soft tissues. This work represents a first step to demonstrate the promise of this type of approach as well as identify some of the next steps to be taken in this technology development.

V. Conclusions

In this study, we proposed to evaluate the capabilities of CNNs to model wave propagation and retrieve elasticity and viscosity values from wave motion images sourced from staggered-grid finite-difference (SGFD) simulations. A total of 1050 SGFD models were simulated using a wide range of elasticity and viscosity values for two different values of excitation F/N. The dataset was then augmented using five different vertical positions to generate the spatiotemporal images and using three different levels of noise. The dataset was used to train, validate and test CNN architectures configured for data regression using mean squared error as loss function. The viscoelastic parameters were also estimated using the conventional 2D Fourier analysis approach with Kelvin-Voigt dispersion curve fitting and compared against the CNN results.

The overall results showed the CNNs’ potential to be an alternative to complex mathematical analyses such as Fourier analysis and dispersion curve estimation used currently for shear wave viscoelastic parameter estimation. The CNN approach does not impose the biases and variabilities associated with these methods and requires minimal amount of preprocessing, such as filters and transforms. The methods proposed might enable simpler and more reliable non-invasive evaluation of tissue injuries that alter rheological parameters such as liver and kidney fibrosis. The use of transfer learning is also a possibility to enable ex vivo and in vivo viscoelastic estimations [43], and this will be pursued in ongoing work.

More evaluations regarding other f-numbers are necessary, although higher values will yield even more planar wave fronts and should provide easier translation to other CNN models. The initial robustness to noise was also encouraging as it is a strong characteristic of ultrasound applications, but other noise settings need to be evaluated. This work serves as a proof-of-concept to demonstrate the use of CNN models for estimation of viscoelastic properties directly from shear wave motion. We will continue to develop this technology to test on experimental data from phantoms and soft tissues.