Ultrasound imaging of lung disease and its relationship to histopathology: An experimentally validated simulation approach

Diagnostic lung ultrasound (LUS) has become a leading real-time imaging technique that embraces both morphological and functional assessment of the lung condition, including critically ill patients.^1–4 Clinical application of point-of-care ultrasound combines differential diagnosis and treatment management in complex cases, including but not limited to acute respiratory failure,⁵ hemodynamic instability,⁶ cardiac arrest,⁷ a variety of connective tissue disorders,⁸ and hemodialysis.⁹ Bedside LUS is considered effective in prediction of successful weaning from mechanical ventilation in adults as well as premature infants with neonatal respiratory distress syndrome.^10,11 The capability of LUS to reveal changes of lung porosity,^12–14 extravascular lung water (EVLW), and lung density is proven in vivo¹⁵ and ex vivo.¹⁶ The role of LUS as a powerful prognostic tool in patients admitted to the emergency department is demonstrated in a wide spectrum of respiratory and non-respiratory diseases and syndromes.^17,18 LUS has become an essential intensive care diagnostic modality because of its relatively low cost and broad availability.^19,20 Strong agreement between LUS and computed tomography (CT) in terms of lung- and lobe-specific lesion location has been demonstrated in a diverse population of patients with acute respiratory failure on mechanical ventilation.²¹ For example, the agreement between LUS and CT regarding the localization of findings within the same lung was 92.5% for the right side and 83.6% for the left side. At the same time, overall lobe-specific agreement between these methods was 87%. Finally, performance of LUS in lesion localization was significantly higher than that of portable chest radiography.²¹ 

The interpretation of LUS is unlike conventional diagnostic ultrasound since ultrasound cannot penetrate through air. In conventional imaging, the brightness of a B-mode image directly represents tissue echogenicity, or, in other words, the travel time between the transducer surface and the underlying anatomical structures is accurately mapped to its spatial representation. LUS, on the other hand, relies primarily on interpreting visual representations of complex reflection, reverberation, and multiple scattering.^22–26 Due to the complex propagation paths, the travel time is lengthened, and diagnostic information is no longer mapped to its anatomical location in space by the beamforming process. Thus, although these acoustical interactions take place in the superficial layer of the lung, they appear deeper than the lung surface in B-mode images,^22–26 the same way as reverberation clutter appears deeper than near-field structures in soft tissue scanning.^27–29 The interpretation of these “artifacts” is a key operator skill required for accurate diagnostic assessment of the disease state of the lung.³⁰ A decrease in lung aeration is accompanied by an increase in the reverberation time constants at the lung surface. This effect results in signals that are prolonged in time in the raw RF data and, thus, appear to be deep in space in the B-mode image. These vertical artifacts (also known as “B-lines” and “comet tails”^23,31–35) in clinical LUS may extended to the bottom of the screen.^32,33 Vertical artifacts are observed in a wide variety of pulmonary and extrapulmonary diseases and syndromes, including but not limited to acute respiratory distress syndrome (ARDS),³⁶ pneumonia,³⁷ pulmonary fibrosis,³⁴ and cardiogenic pulmonary edema.³⁸ At the same time, vertical artifacts, despite their high sensitivity, have low specificity in differentiation between the lesions' causes and represent sonographic interstitial syndrome.³⁹ 

From the perspective of ultrasound propagation, lung diseases can be classified into two categories of alterations of lung tissue. First, air is replaced by fluid, which typically occurs in the acute phases of disease, e.g., in ARDS. Second, non-aerated lung components (interstitium and exudate or transudate) are replaced by fibrous tissue, which typically occurs in the chronic phases of disease. Taken on its own, the replacement of air by fluid has a more significant acoustical effect, since air is practically totally reflective due to its high impedance mismatch with surrounding tissue. However, both categories of disease modify the alveolar microstructure and its properties, such as mean size, distribution, and organization.^40,41

Pulmonary fibrosis is a well-defined consequence of ARDS, which is manifested by irreversible morphological changes in lung parenchyma and decline in its compliance and respiratory function and is finally associated with poor long-term prognosis.^42,43 Despite dramatic reduction of long-term pulmonary complications of ARDS with the advent of lung-protective ventilation, patients diagnosed with COVID-19 and developed ARDS constitute a high-risk group for pulmonary fibrosis.⁴⁴ Fibroproliferative phase of ARDS is characterized by fibroblast proliferation, which is a common response to lung injury without dependence on its dissemination or initial injurious agents.^40,45 The patterns of reaction observed histologically and radiologically in patients with ARDS are similar, regardless of the underlying cause.⁴¹ Since lung fibrosis has gradual development and may coexist in time with exudative phase, and specificity of mentioned LUS artifacts is low, sonographic assessment of fluid-to-fibrotic tissue ratio in these patients is challenging.^46,47 Pleural line measurements (its thickening) and subpleural nodules on LUS are defined as the most efficient sonographic features in early diagnosis and assessment of progression in interstitial lung diseases.² However, CT remains the gold standard in clinical diagnostics of pulmonary fibrosis because, unlike ultrasound, the CT scan can image the entire lung.⁴⁸ A variety of modern technical studies in clinical ultrasound is focused on the quantitative sonographic assessment of lung in pathology and its automation.^31,49–57 Spectroscopic imaging and analysis performed on lung phantoms and patients with lung fibrosis revealed difference in vertical artifacts' intensity depending on the fundamental frequency.^58–61 Ultrasound multiple scattering characterization of the normal rat lung and modeled acute edema using calculation of the diffusion constant and transport mean free path facilitated reproducible differentiation of the presence and degree of pulmonary edema.⁶² Histopathology and high-resolution CT have been widely used as validation methods for diagnostic features observed with ultrasound imaging and non-imaging characterization techniques.^2,13,62–64

Due to its complexity, ultrasound propagation in the lung and its relationship to diagnostic assessment remain poorly understood, and resolving this relationship is a long-standing goal for the field. Numerical simulations of ultrasound propagation in homogeneous media with artificially coded air inclusions revealed the impact of domain geometry and wavelength on the visualization of the vertical artifacts.⁶⁵ Using segmented slides of human upper abdomen and synthetic air inclusions, basic LUS views were modeled, demonstrating the signs characteristic for normal and diseased lung with exudative lesion and a quantitative link between the percentage of fluid in a lesion and the contrast of vertical artifacts established.⁶⁶ Two-dimensional (2D) finite-difference numerical simulations based on an animal histopathology have been performed, and significant correlation between the lung parenchyma aeration and calculated transport mean free path was observed.⁶² 

The aim of this work is to establish a quantitative link between human lung anatomy, its microstructural organization, and the resulting ultrasound images using a physics-based wave propagation approach. This includes (1) demonstrating the ability to accurately simulate ultrasound imaging pulses and their propagation in normal and diseased human lung, (2) generating in silico (numerically simulated) conventional B-mode images that are comparable to clinical in vivo ones, and (3) quantitatively and qualitatively characterizing the impact of independent changes of aeration and fluid-to-fibrotic tissue ratio in the lung parenchyma based on ultrasound images. This is achieved by using an acoustic simulation tool called Fullwave-2⁶⁷ in conjunction with the generation of highly realistic acoustical properties derived from the National Library of Medicine's Visible Human data set⁶⁸ and human lung histology. The anatomical chest wall maps and microscopic lung structures are combined to generate acoustical maps of the human body. The porosity, or ratio of fluid and connective tissue, is quantified in both pathological and normal lung specimens. Simulations of diagnostic ultrasound imaging pulses propagating in the tissue are then conducted.

Digital modification of the alveolar microstructure is used to model disease progression in a way that offers fine control over the acousto-morphological relationship at the microscopic level. This approach offers the ability to control the microstructure and determine relationships between LUS and disease in a way that would not be feasible clinically. Gradual decrease in porosity and fibrosis development are modeled in silico, artificially changing the portion of air, fluid, and connective tissue elements (pixels) in the acoustical maps. Based on the calculated radio frequency (RF) data, conventional B-mode images are generated, and pathological cases are scored in accordance with current clinical recommendations.^37,69 Pulse propagation depth is assessed in normal lung via simulation of hydrophone measurements at the surface and inside the parenchyma. Each grade of fluid retention and modeled progressive lung fibrosis are characterized by imaging [regional lung ultrasound score, brightness of region of interest (ROI) in subpleural region] and non-imaging reverberation measures (reverberation time constant, decay coefficient, and characteristic length).

Accurate modeling of ultrasound wave propagation in biological materials incorporates physical phenomena, such as diffraction, reflection, multiple scattering, frequency-dependent attenuation, and nonlinearity. Acoustic fields in a nonlinear attenuating heterogeneous medium with large contrast can be determined using staggered-grid finite-difference method that is characterized by higher stability and accuracy and lower computational cost compared to the methods implemented on a conventional grid.^70–73 Here, we use Fullwave-2,⁶⁷ which is an evolution of the Fullwave solver.⁷¹ The Fullwave solver was extensively used to simulate propagation of imaging ultrasound pulse through abdominal wall,⁷⁴ focused pulse in transcranial brain therapy,⁷⁵ and modeling of traumatic brain injury mechanisms.⁷⁶ The tool was experimentally validated in vivo in terms of reverberation and phase aberration metrics⁷⁷ along with effects of attenuation, scattering, and absorption.⁷⁸ The physical equations are summarized briefly here.

The wave equation is represented as a pressure-velocity system of equations, including the acoustic pressure characteristic for specified location and time instant; the material density and compressibility, which are the functions of space; a velocity wavefield vector; and the non-linearity parameter (B/A). Dispersion and attenuation are modeled by coordinate transformation of the spatial derivatives, specifically by the partial derivative scaling and summation of the convolutions of N relaxation functions. The linear scaling of the derivative outside the convolution is represented by parameter $κ x 1$⁠, which changes the sound speed in a particular dimension and is a function of space. Transformations are performed along each of two dimensions, and multiple relaxation parameters are determined as maps as a function of space to impose spatially and frequency-dependent attenuation laws and perfectly matched layers. The numerical solutions and their validation are thoroughly described in Ref. 67. The 2D simulations described in this work employed numerical solver Fullwave-2, while air inclusions were represented by elements of constant zero pressure in the initial condition maps.

A linear transducer L12-5 50 mm (ATL, Bothell, WA) is simulated, emitting a 6.25 MHz, two-cycle, 2.5 MPa focused pulse, which can be used with the research ultrasound system Vantage 256 (Verasonics) or a variety of Philips machines. This array is consistent with clinical LUS standards and recommendations, although a variety of arrays can be used.^32,37,79,80 The focal depth was set to the pleural line level (20 mm). Additionally, an investigation of the focal depth was performed across a range of 15–40 mm to determine the reduction in image quality when it does not match the pleural depth. The ultrasound imaging sequence was composed of 128 independent transmit-receive events; in each, the subaperture of 64 elements is shifted laterally by 0.195 mm (equivalent of a walking aperture scheme). RF signals were received with the same 64 transmitting elements of each subaperture. Lateral and axial resolution values were estimated as 0.375 and 0.246 mm, respectively. Calculations were based on the wavelength of 0.246 mm and f-number restriction to 1.6 (for 20 mm, depth subaperture was 12.5 mm) and took into account the two-cycle pulse length. The simulation field for each transmit-receive event was 25 mm wide and 35 mm deep. The total duration of the simulations (131.8 μs) is determined in compliance with a recommended imaging depth range of 4–8 cm starting from the pleural line,³² so it is possible to visualize at least the second reflection of the initial pulse from the lung surface (in normal LUS, first reflection corresponds to the pleural line, and the second one is the first horizontal artifact, which is an equally spaced replica of the pleural line³⁰). The spatial step size is calculated as 20.5 μm, relative to the 6.25 MHz transmit fundamental frequency (λ = 0.246 mm) and arbitrarily chosen discretization of 12 points per wavelength. The temporal step is set to 6.7 ns for a reference speed sound of 1540 m/s, which corresponds to a Courant–Friedrichs–Lewy condition of 0.5.⁸¹ 

Simulations of focused imaging ultrasound pulse propagation in a normal and diseased human lung were performed by modeling oblique lung scanning (with the transducer located along the intercostal space).⁶⁹ The ability to generate in silico parasagittal views, including the ribs, was demonstrated earlier.⁶⁶ In these simulations, air was modeled with elements of constant zero pressure (no sound wave propagation), which imposes a total reflection of the acoustical wave from the tissue/air interface in accordance with experimental data of Oelze et al.^82,83 Indeed, it releases strict limitations of spatial element size and enables simulation of the pulse propagation in tissue only, not in low sound speed and density air.

Scanned body tissue may be roughly divided into two complex layers (Fig. 1): chest wall [Fig. 1(A)] and the lung parenchyma covered with the pleura [Fig. 1(B)]. The former comprises skin, adipose tissue, blood vessels and lymphatics, ribs (not shown here), muscles, and fascia, while the latter in subpleural regions in healthy individuals is a highly aerated complex network of thin alveolar walls (capillaries covered with respiratory epithelium) and blood vessels and lymphatics surrounded by an elastic connective tissue framework.^84,85 Superimposition of these two layers was employed in this study [Fig. 1(C)]. While pleural effusion takes place in most cases of acute lung injury,⁸⁶ fluid of water acoustical properties was added between parietal and visceral pleura in pathological cases, as is shown in Fig. 5(B).

Acoustical maps of chest wall [Fig. 1(A)] were derived from optical data of the National Library of Medicine's Visible Human data set,⁶⁸ which have a resolution of 330 μm. These digitized images were segmented into four types of tissue (connective, adipose, muscle, and blood) using tissue-specific thresholds.⁸⁷ The maps of lung parenchyma and pleura in pathology [Fig. 1(B)] were derived from deidentified educational lung histopathology images of higher resolution (0.35 μm) provided by the Department of Pathology and Forensic Medicine, Dnipro State Medical University (Dnipro, Ukraine). The maps of normal inflated swine lung parenchyma and pleura were derived from histological bright-field microscopy images (section thickness 5 μm, hematoxylin and eosin staining) scanned with slide scanner SlideView^TM VS200 (Olympus Corporation, Tokyo, Japan) at 10× magnification, resolution 0.55 μm (Hooker Imaging Core, University of North Carolina at Chapel Hill). Swine lung tissue was harvested as part of the Tissue Sharing Program (North Carolina State University) and processed by Pathology Services Core, University of North Carolina at Chapel Hill. The images from all three sources were interpolated using nearest neighbor interpolation, taking into account the difference in their original spatial resolution and the size of the simulation elements.

Lung histopathology images were binary segmented (air/non-air components). The segmentation into different tissue types (connective tissue, muscle, fat, blood) was performed for chest wall images, and the corresponding acoustic properties (speed of sound, density, attenuation, and nonlinearity coefficients)^88,89 were assigned to each of the spatial elements of the formed 2D acoustical maps. To realistically model the acoustical scattering in soft tissue, scatterers of 20.5 μm diameter and average impedance mismatch of 6% relative to the background tissue were generated and randomly added to the sound speed maps. The scatterer density was set at 12 scatterers per resolution cell of size 0.375∗0.246 mm [lateral∗axial]. The simulation calibration (reverberation and speckle brightness values) was extensively described earlier.⁷⁷ 

The acoustical maps derived from lung histopathology were analyzed, and lung aeration was calculated planimetrically as a percent of air elements out of all lung parenchymal elements. Size of alveoli and alveolar wall thickness were estimated via linear intercept using custom matlab (Mathworks, Natick, MA) code reproducing the standard method.⁹⁰ To study the link between spatial distribution of acoustical channels and vertical artifacts, linear (axial and lateral acoustical channel size) and planimetric (area of non-aerated medium) indices were calculated [Fig. 2(A)]. For the manipulations with tissue geometry and its acoustic properties, the tissue/air interfaces and each air inclusion were localized automatically. Gradual increase in lung aeration was modeled through expansion of air inclusions in all directions by a single element increment (20.5 μm). In the basic simulations of LUS in ARDS, the parenchymal non-air elements were assigned generic fluid acoustical properties, corresponding to EVLW (Table I). In all the simulations, pleura has been segmented and assigned connective tissue properties (Table I).

To model and analyze gradual evolution of LUS images characteristic for progressive fibroproliferative disease (lung fibrosis), parenchymal elements were transformed from fluid to connective tissue properties, while aeration was kept at constant level. These substituted connective tissue elements are chosen to be located at equal distance between adjacent tissue/air interfaces, which leads to location of the remaining fluid elements (if any) around the air inclusions and gradual expansion of fibrosis from central part of each tissue wall. Such distribution of fibrotic tissue reflects extracellular matrix deposition as a result of fibroblast proliferation and activation.^41,45

The RF data were generated, and a conventional delay-and-sum beamforming algorithm was used. The received signals were sampled at a rate of 150 MHz and then downsampled to 25 MHz, which is appropriate for clinical ultrasound imaging systems. No filtering or interpolation was applied to B-mode images, except for time gain compensation (TGC) with an exponential depth/time-dependent gain based on an attenuation coefficient of 0.5 dB/cm/MHz, assuming the tissue sound speed of 1540 m/s. The image processing was done in matlab. The simulations were performed on a Linux-based computer cluster running Intel Xeon^® E5–2630 v4 processors at 2.20 GHz, using the mentioned custom Fullwave-2 simulation tool. Each of 256 individual simulations necessary to form a single B-mode image took 6 h and was run in parallel. The source code is located in the GitHub repository.⁹¹ 

To assess the relations between acoustical channels' size and vertical artifacts' detectability on the simulated B-mode images (taking into account their distribution), Spearman's correlation coefficient was calculated. For assessment of detectability of the artifacts or structures, the generalized contrast-to-noise ratio (gCNR) was calculated as $gCNR = 1 − OVL$⁠, where OVL is the overlap area between probability density functions of false detection Pf and miss Pm.⁹² For quantification of changes in B-mode images associated with shifting of the parenchymal acoustic properties, averaged brightness of subpleural ROI was calculated, and hydrophone measurements were simulated at the pleural surface. The simulated LUS B-mode images were scored in accordance with Refs. 37, 79, and 80.

The feasibility of using the merged histological and anatomical acoustical maps and Fullwave-2 numerical simulation tool to generate LUS images with a quality that is comparable to clinical ones of normal and diseased human lung was assessed. The acoustical map of normal lung derived from a swine histology image was analyzed planimetrically, and the baseline parenchyma aeration was calculated as 70.7%. For reference, the average porosity of the normal human lung parenchyma is within the range of 70%–90%.^93–95 Alveolar size estimated by linear intercept was 132 ± 24 μm, and alveolar wall thickness was estimated as 15 ± 1.5 μm, which is consistent with literature data.⁹⁶ Non-alveolar spaces (such as bronchovascular bundles) were excluded in both estimations. While the normal alveolar wall thickness in the used swine lung image was 15 ± 1.5 μm and the other segmented structures in both histology sources were significantly larger and took into account reasonable trade-off between spatial resolution and computational cost, a spatial step size of 20.5 μm was considered as appropriate for performed simulations.

To investigate the acoustical behavior of normal lung on the image characteristics and on the penetration length scales into the parenchyma, simulations were performed using focused imaging at a depth of 20 mm. Received echoes at the transducer surface were used to generate B-mode images. A conventional approach consisting of a single transmit per A-line of generated B-mode image was used. To compare the simulation based on the histology of a normal inflated swine lung and horizontal specular reflectors, both types of air inclusions were combined in one physical map [Fig. 2(B)]. The histology map that was inserted in the anatomical map of tissue properties was shown in Fig. 1(B). To the left and right of the histology inset, the tissue/air interfaces are flat or specular [Fig. 2(B)]. Within the inset box, the acoustical properties are described by the microstructural geometry.

The beamformed B-mode image [Fig. 2(C)] exhibits the hyperechoic pleural line, and equidistant horizontal artifacts are created due to the dominant superficial sound reflection. This is consistent with clinical LUS of healthy tissue,⁵ with previous ex vivo¹⁶ and numerical simulation studies.^65,66 The resulting generated B-mode image, if assessed by a sonographer, would show an insignificant difference between the appearances of characteristic LUS signs, including the pleural line and horizontal artifacts in both its central (histology) and lateral (straight tissue/air interfaces) aspects [Fig. 2(C)]. In addition, on the signal plot of a superficial (red) hydrophone, decrease in peaks' amplitude in time and their elongation can be appreciated [Fig. 2(A)]. This pattern is visible in B-mode image [Fig. 2(C)] as decrease in brightness and contrast of the horizontal artifacts (appearing at 40 and 60 mm depth) too. Both findings are consistent with existing ideas about acoustical channels and traps and reflect the deeper propagation of the ultrasound pulse inside the lung parenchyma as well as its depth-dependent retention and attenuation.³⁰ 

The same acoustical map and focused imaging transmits were used to determine the penetration of ultrasound in this healthy lung specimen. A 16 mm wide line of virtual hydrophones was placed at a depth of 20 mm and at 21 mm as shown in Fig. 2(B). These two depth locations measure the sound entering the histology specimen at the pleural surface [Fig. 2(B), red] and then reaching 1 mm depth of the lung specimen as a function of time [Fig. 2(B), blue].

The averaged signal amplitude at 21 mm of depth was 67.9 dB lower than at 20 mm of depth. This indicates that in these simulations, the sound penetration in healthy lung is, as expected, low due primarily to the significant air content and high impedance mismatch. Furthermore, although the tissue layer between virtual hydrophones is thin, it is sufficient to model the acoustical behavior in healthy lung tissue. Subsequent results will demonstrate that for diseased tissue, larger thicknesses are necessary.

For quantitative analysis of subpleural lung region and specifically connections between pleura and connective tissue septa (entrances to the acoustical channels), lateral and axial size of entrance were measured 0.2 mm deeper than pleura surface [Fig. 3(A)]. These two indices show the length of an incident beam segment before its first reflection inside the lung parenchyma and were calculated using unidirectional segmentation of non-air elements [Fig. 3(B)]. Acoustical channel entrance area was calculated as a sum of all non-air elements adjacent to the beam elements at depth range from 0.2 mm deeper than pleura surface to the entrance axial size. This area shows projection of the acoustical channel entrance on the elevational plane. The plotted acoustical channel entrance area [Fig. 3(C)] demonstrates four dominant acoustic channels spreading deeper than normal pleural thickness (0.2 mm). Note the higher axial size and entrance area characteristic for the channel that is coaxial to the incident beam. To assess vertical artifact detectability, each A-line of B-mode image at a depth range from 0.5 to 18.5 mm deeper than the pleural line was analyzed [Fig. 3(D)]. Averaged normalized signal level of A-lines did not demonstrate consistent relation to lung specimen and acoustic channels [Fig. 3(E)]. However, a plot of gCNR [Fig. 3(F)] calculated for each A-line at the mentioned depth range shows the highest peak of its detectability (gCNR = 0.71) at lateral position (+3 mm) matching the channel of the largest entrance area (0.52 mm²). The general increase in gCNR in a specimen region may be explained by a different distribution of brightness (signal level) compared to the regions of specular reflectors, although average signal levels are comparable [Fig. 3(E)]. The similarity and difference between healthy lung and specular reflectors in terms of reflectivity and ultrasound pulse retention also can be visually appreciated on the pressure field snapshots at the moments before [Fig. 4(A)] and after [Fig. 4(C)] pulse reflection. Here, it can be seen how the pulse is reflected by both synthetic specular tissue/air interfaces and irregular lung surface. At the same time, shallow propagation into the acoustic channels of lung parenchyma (lateral position −3 mm and 3 mm) is evident [Fig. 4(C)]. Higher level of signal propagating into the coaxial connective tissue septum (lateral position 3 mm) can be visualized. The latter observation corresponds to the results of quantitative assessment of this case.

Diseased lung tissue was used to relate the acoustical behavior, B-mode image characteristics, and penetration characteristics into the parenchyma as a function of disease state. In the case of pathology, lung parenchyma demonstrates patterns that are characteristic for specific lesion.⁴⁰ The human lung histopathology [Fig. 5(A)] in this sample exhibits severe combined interstitial and alveolar edema, vascular congestion, and decreased aeration. Parenchyma aeration was estimated using segmentation analysis performed with matlab and was estimated to be 21.0%, which is considered a pronounced decrease in lung porosity compared to healthy lung.⁹⁴ Alveolar size was calculated as 131 ± 23 μm (normal 216 ± 28 μm). Interalveolar (border-to-border) distance was 741 ± 81 μm (normal 15 ± 1.5 μm). This combination of linear and planimetric indices shows not only significant alveolar de-recruitment, but also complete collapse or fluid filling of some alveoli and alveolar ducts. The pleura represented 0.9% of all (air and non-air) lung elements in this image. Binary segmented acoustical maps were generated based on this image. To model LUS in exudative phase of ARDS, the non-air elements were assigned generic fluid acoustical properties, corresponding to EVLW, the pleura was assigned connective tissue properties, and air was assigned constant zero acoustic pressure throughout the simulation (Table I). The lung histopathology map was inserted laterally in the center and at the same depth of 20 mm with the addition of the specular reflector, consisting of flat interfaces, to the left and right. Simulation of the same focused imaging scheme was performed (see Sec. II B).

The simulated B-mode image using this map demonstrates multiple vertical artifacts and irregularity of the pleural line, corresponding to a LUS score of 2 [Fig. 5(C)]. The capability to model beam focusing at various depths and relevant changes of obtained B-mode images is demonstrated in Fig. 5(D), where doubling the focal depth leads to a limited ability to detect separate vertical artifacts visually. Area under the curve calculated in gCNR plots in the specimen's lateral position (red plot) was 0.73 in the case of 2 cm focal depth [Fig. 6(F)] versus 0.63 (p < 0.05) in the case of substantial mismatch between pleural location and focal depth [Fig. 6(H)], which indicates lower detectability of vertical artifacts in the second case. Evolution of B-mode images related to stepwise change of focal depth and their contrast metrics are shown in the  Appendix. These data are consistent with literature data received in in vitro experiments.^60,97 In both cases, the image can be interpreted as sonographic interstitial syndrome.^39,80 Quantitative analysis of the subpleural lung region and simulated B-mode image was performed as described in Sec. III A. Its results are illustrated in Fig. 6(C). In comparison with the normal lung case, these data indicate significantly higher channel size that exceeds lateral resolution of the simulated sonographic system, wider channel entrance area throughout the specimen, and higher detectability of vertical artifacts (gCNR from 0.5 to 0.95) related to its spatial position [from –9 to 9 mm; Fig. 6(F), red plot]. The correlation analysis of the data received in normal and severely abnormal cases with the same focused imaging scheme and chest wall showed moderate correlation between gCNR and lateral size of channel entrance (r = 0.69, p < 0.05), moderate correlation between gCNR and axial size of channel entrance (r = 0.7, p < 0.05), and strong correlation between gCNR and estimated area of channel entrance (r = 0.71, p < 0.05).

Synthetic parenchyma re-aeration via expansion of air inclusions in the acoustical maps described in Sec. II B was used to model the alveolar recruitment and its effects on LUS images (Fig. 7). Parenchymal porosity ranging from 21.0% [Fig. 7(A)] to 80.3% [Fig. 7(F)] was considered. The lower end of this range represents the original state of the lung, which is severely diseased, whereas the upper aeration range is consistent with levels in healthy lung.

Simulations of ultrasound propagation in these maps and the resulting B-mode images showed dramatic changes of vertical artifacts accompanied by stepwise formation of horizontal artifacts as the aeration increases. This is also reflected in a decrease in the LUS score from 2 [Fig. 7(A)] to 0 [Fig. 7(F)].

Dramatic differences in reflection, absorption, and back emission of the diagnostic ultrasound pulse depending on the lung parenchyma aeration can be appreciated in the pressure field images at the same time instances before and after lung tissue insonification (Fig. 8).

Progression of proliferative and fibrotic phases of ARDS was modeled in the acoustical maps, as described in Sec. II B. Parenchymal aeration was maintained at the initial level of 21.0%. The EVLW in the parenchyma was gradually replaced with connective tissue in the range of 0%–78.1% as measured by the total fraction. In other words, at the highest percentage, the EVLW was totally replaced by connective tissue, and the parenchymal acoustical map consisted of either connective tissue (78.1%) or air (21.0%).

The 12 simulations were performed, and corresponding images are shown for side-by-side comparison (Fig. 9). They demonstrate faster attenuation of the vertical artifacts as the connective tissue fraction increases. This can be explained by a combination of two factors: (1) the significant distance that reverberating wave travels inside of the parenchymal acoustic trap and (2) the higher attenuation of connective tissue compared to the EVLW (Table I). No difference in the lateral distribution of vertical artifacts was found compared to the images based on acoustic maps without artificial addition of connective tissue (Fig. 5). The pleural line irregularity and thickening, as known LUS signs of pulmonary fibrosis,^2,98,99 are observed in the cases with the wider connective tissue spreading [Figs. 9(G)–9(L)].

One of the strengths of the simulation approach is that it allows the local quantification of reverberation lengths. This quantitative measurement was performed by inserting a virtual hydrophone at the reverberating location and quantifying fundamental physical properties of the decay rates. As such, this physical quantity is a quantitative measure that is dependent on the local wave physics and its interaction with the material (unlike quantification in B-mode images, which suffers from substantial propagation confounds). The subpleural ROI brightness was quantified, and hydrophone measurements at the lung surface were calculated for each simulation. The ROI of size 10 × 10 mm was arbitrarily chosen in the lateral center of beamformed images at depth 30 mm [Fig. 10(A)]. For quantitative characterization of back emission related to reverberation inside the lung parenchyma, the simulated hydrophone data were fitted, and exponential decay of signal amplitude was observed. Finally, the corresponding decay coefficient, reverberation time constant, and estimated characteristic length were calculated for each simulation (Table II). It is shown that the growth of connective tissue representation in parenchyma is accompanied by the decrease in averaged brightness of the vertical artifacts by 3.7 dB (Fig. 10), consistent decrease in reverberation time constant (Fig. 11), and characteristic length. Both plots were inherently mirrored because parenchyma in these simulations is constituted of only three components (connective tissue, fluid, and air), and aeration is kept constant; thus, change of connective tissue amount led to proportional opposite change of fluid amount.

Ultrasound imaging in the lung was modeled based on the fundamental physics of wave propagation in anatomically realistic maps of body wall and tissue properties. The simulation tool that solves the full wave equation in heterogeneous media (Fullwave-2) established this synthetic data set that connects modifications of lung microstructure based on the diseased state to ultrasound imaging. The simulations of 6.25 MHz linear array transducer L12-5 50 mm were performed consistently with published clinical LUS standards and recommendations. These resulting B-modes are comparable to clinical observation in terms of key diagnostic indicators, such as anatomical features (bright pleural line, visible chest wall muscle, fat, and connective tissue components) and artifactual LUS signs (horizontal artifacts in highly aerated lung and vertical artifacts in injured lung).

The realistic acoustical maps were established by combining photographic cryosections of the human body wall with human histology and pathology slides. This approach has the advantage of preserving the soft tissue anatomy in the body wall while representing the microscopic alveolar structure that is critical to modeling the propagation and multiple scattering interactions in the lung. It enabled the description of tissue and air geometry as close as possible to realistic scenarios, both in normal and pathological conditions. High spatial resolution of performed models made it possible to compare simulated B-mode images based on normal inflated lung histology and specular reflector (straight tissue-air interfaces). The porosity of the lung parenchyma was calculated planimetrically for both normal and pathological cases. Ultrasound propagation in normal lung parenchyma is modeled and characterized by scanning a virtual hydrophone at different depths and locations inside the lung. Such measurements would not be realizable in vivo. The demonstrated-in-simulations changes of general reflectivity and pulse propagation in lung tissue, depending on its aeration, are consistent with significant fluctuation of pressure reflection coefficient and acoustic impedance values determined experimentally.^82,83 The impact of the focal depth on diagnostically important features (detectability of vertical artifacts and pleural line geometry) was also demonstrated in silico. Alveolar recruitment was modeled along with corresponding LUS images. The observed horizontal and vertical artifacts and their distribution were consistent with changes of modeled lung porosity and made it possible to score the images according to the current point-of-care ultrasound recommendations.⁸⁰ Particularly, the highest aeration characterized by regional LUS score = 2 was 35.4%, and the next shift from score = 1 to score = 0 took place at a minimum lung aeration of 63.9%, which was consistent with existing clinical and experimental data.^{12,13,35,39,64,66} However, the regional LUS score as a semiquantitative metric of EVLW has limited sensitivity and even lower specificity.^{2,5,19,46,61,100} Moderate correlation between lineal size of acoustic channel entrance and vertical artifact gCNR was observed (r = 0.69, p <0.05).

It was possible to model gradual development of fibroproliferative lung disease as a late phase of ARDS, since the simulation tool provided complete control over the medium, represented by acoustical maps. Signal amplitude at the surface of lung (and corresponding brightness of vertical artifacts) was lower by 3.7 dB in the case of maximal connective tissue representation compared to the case of maximal fluid portion in parenchyma, which corresponds to clinical observations and published experimental data.^3,61,98 Back emission arising from reverberation and multiple scattering inside the diseased lung parenchyma was quantitatively characterized. It was shown that reverberation time (related to the length of visible vertical artifacts) decreases consistently from 56 to 30 μs with increase in the connective tissue portion from 5.3% to 78.1% (and proportional decrease in fluid amount) in the lung tissue.

The main limitations of the described approach are (1) difficulty of receiving large pathology specimens/slides of the human lung containing intact pleura; (2) performance of all analysis planimetrically (2D), although the lungs have a complex three-dimensional architecture; (3) performance of simulations on stationary acoustical maps, not reflecting respiratory dynamics (changes of position and alveolar geometry) that takes place in real-time scanning of lung in vivo; and (4) analysis of only two histological images, although there is a broad structural variability for both normal swine and diseased human lung and their corresponding ultrasound images.

In conclusion, this work establishes a link between diagnostic ultrasound images, the ultrasound propagation physics, vertical artifacts' detectability, quantitative characteristics of superficial acoustical channels, and underlying lung pathology, including parenchyma aeration and acoustic properties of tissue. This assessment would be difficult or impossible to obtain in clinical or in vivo experimental settings, where the ground truth represented by the lung microarchitecture is not measurable.

The ability to diagnose and assess gradual development of fibroproliferative lung disease on the background of acute exudative lung lesion is relevant not only in ARDS, but also in cases of preexisting heart and lung conditions (interstitial lung diseases; chronic left ventricular congestive heart failure, especially with preserved myocardial contractility; congenital and acquired obstructive left heart diseases; chronic obstructive pulmonary disease). The obtained results demonstrate one of the possible efficient ways to research imaging and non-imaging ultrasonographic characteristics in a variety of lung diseases in silico, especially related to changes of lung aeration and its tissue composition. A prospective opportunity for the described method is generation of multiple serial diagnostic RF data sets using human lung histology as a representation of specific known disease. These data could be used for development and testing of new sequences in diagnostic ultrasound and specialized image-formation techniques to improve their diagnostic value and performance. Employment of machine learning in these tasks would empower faster transition from qualitative and semiquantitative, operator-dependent diagnostic metrics (such as LUS score) to quantitative reproducible indices for assessment of EVLW, morphological extent of the lesion, and its tissue composition.

We thank Gabriela De la Cruz (Animal Histopathology Core, School of Medicine, University of North Carolina at Chapel Hill) for assistance with processing of normal swine lung histology specimens. This research was supported by the National Institutes of Health (Grant No. R21-EB033150).

Here, we demonstrate in detail visual (Fig. 12) and quantitative (Fig. 13) changes of B-mode images related to stepwise increase in focal depth from 15 to 40 mm in a single case of severe combined lung disease (aeration 21%). Focusing at pleural line depth (20 mm) led to the highest visibility, contrast, and contrast-to-noise ratio of vertical artifacts.

Contrast was calculated as $C = 20 log 10 ( S i / S o )$⁠, where S_i and S_o are the brightness levels measured inside and outside the ROI (histology location), respectively. Contrast-to-noise ratio (CNR) was calculated in accordance with Patterson and Foster¹⁰¹ as $CNR = ( S i − S o ) / ( σ i 2 + σ o 2 )$⁠, where σ_i and σ_o are the standard deviations of the image brightness inside and outside the ROI, respectively.