Quantifying Regional Lung Deformation Using Four-Dimensional Computed Tomography: A Comparison of Conventional and Oscillatory Ventilation

Mechanical ventilation strategies that reduce the heterogeneity of regional lung stress and strain may reduce the risk of ventilator-induced lung injury (VILI). In this study, we used registration of four-dimensional computed tomographic (4DCT) images to assess regional lung aeration and deformation in 10 pigs under baseline conditions and following acute lung injury induced with oleic acid. CT images were obtained via dynamic axial imaging (Siemens SOMATOM Force) during conventional pressure-controlled mechanical ventilation (CMV), as well as high-frequency and multi-frequency oscillatory ventilation modalities (HFOV and MFOV, respectively). Our results demonstrate that oscillatory modalities reduce intratidal strain throughout the lung in comparison to conventional ventilation, as well as the spatial gradients of dynamic strain along the dorsal-ventral axis. Harmonic distortion of parenchymal deformation was observed during HFOV with a single discrete sinusoid delivered at the airway opening, suggesting inherent mechanical nonlinearity of the lung tissues. MFOV may therefore provide improved lung-protective ventilation by reducing strain magnitudes and spatial gradients of strain compared to either CMV or HFOV.


INTRODUCTION
Ventilator-induced lung injury (VILI) may inadvertently occur in critically ill patients receiving mechanical ventilation, due to the harmful stresses and strains associated with gas flows driven under positive pressure (Slutsky and Ranieri, 2013). Patients with injured, inflamed, and/or edematous lungs, such as those with the acute respiratory distress syndrome (ARDS), are especially at risk for VILI due to increased ventilation heterogeneity (Carrasco Loza et al., 2015). Maintaining normal levels of arterial oxygen and carbon dioxide tensions in patients with ARDS may be extremely difficult using conventional mechanical ventilation (CMV), since increasing the tidal volume or driving pressure to compensate for gas exchange deficiencies may be detrimental to the mechanically overburdened lung (Acute Respiratory Distress Syndrome Network, 2000; Amato et al., 2015;Gattinoni et al., 2016a). The goal of lung-protective ventilation is to provide life-sustaining gas exchange without exacerbating existing injury within vulnerable parenchymal tissues.
In some instances, high-frequency oscillatory ventilation (HFOV) has been proposed as a rescue treatment for refractory hypoxemia in ARDS, given its theoretically ideal qualities for lung-protection: small tidal volumes that mitigate the risk of dynamic strain injury (i.e., volutrauma), and high mean airway pressures that prevent cyclic recruitment/derecruitment (i.e., atelectrauma) (Sklar et al., 2017). Despite an extensive history of scientific and clinical research over several decades, optimal strategies for HFOV initiation and management remain a subject of controversy (Malhotra and Drazen, 2013;Kneyber and Markhorst, 2016;Nguyen et al., 2016). As it is currently delivered, HFOV may not be an appropriate ventilatory modality in many patients with ARDS for several reasons. First, the use of high mean airway pressures may result in hemodynamic compromise due to impaired venous return and low cardiac output (Meade et al., 2017). In addition, high strain rates within the lung tissues during HFOV may also contribute to VILI, by increasing the mechanical power dissipated across the parenchyma (Gattinoni et al., 2016b;Protti et al., 2016;Herrmann et al., 2019b). Finally, high frequency oscillatory flows are distributed in a heterogeneous and frequencydependent manner, predisposing overventilated regions to excess mechanical strain, and underventilated regions to derecruitment and atelectasis (Amini and Kaczka, 2013;Herrmann et al., 2016).
Multi-frequency oscillatory ventilation (MFOV) has thus also been proposed as an alternative modality for lungprotective ventilation, by delivering multiple frequencies of oscillatory flow and pressure simultaneously (Kaczka et al., 2015;Herrmann et al., 2019b). MFOV has been shown to improve gas exchange and mechanical function compared to traditional "single-frequency" HFOV (Kaczka et al., 2015;Herrmann et al., 2018). Its postulated mechanism of benefit relies on the frequency-dependence of ventilation distribution and gas transport during conventional and oscillatory modes of ventilation , which can be exploited by using flow waveforms with enhanced harmonic content. Given such frequency-dependence of gas flow throughout the airway tree, mechanically disparate regions of an injured lung may thus selectively filter out "less desirable" frequency components of a broadband oscillatory waveform delivered at the airway opening. The corresponding homogenization of intrapulmonary gas transport with MFOV may then contribute to enhanced gas exchange, and simultaneously distribute the mechanical burden of ventilation more uniformly throughout the lung. Nonetheless, specific experimental evidence to support such a mechanism is lacking, and it remains to be seen whether MFOV improves the regional distribution of parenchymal strain throughout the lung compared to CMV or traditional HFOV.
The purpose of this study was to characterize the distribution of regional intratidal lung deformation during conventional and oscillatory modes of ventilation in pigs under baseline conditions and following acute lung injury with oleic acid. We hypothesized that such variations would be less heterogeneous during MFOV, compared to CMV or traditional singlefrequency HFOV. Specifically, we quantified intratidal variations in dynamic regional lung aeration, volumetric strain, and volumetric strain rate during CMV, HFOV, and MFOV using frequency-selective four-dimensional computed tomographic (4DCT) image reconstruction  and registration (Zhao et al., 2016).

Animal Preparation
All experimental procedures were approved by the University of Iowa Institute for Animal Care and Use Committee (Protocol Number 5061428). Ten healthy pigs were used in this study, weighing between 9 and 13 kg. Following intramuscular injection of 2.2 mg kg −1 telazol, 1.1 mg kg −1 ketamine, and 1.1 mg kg −1 xylazine, general anesthesia was induced during spontaneous breathing with inhaled isoflurane delivered via nosecone. Each pig was then intubated with a cuffed endotracheal tube (4.5-5.5 mm inner diameter), and mechanically ventilated (Uni-Vent Eagle Model 754, ZOLL Medical Corporation, Chelmsford, MA). Capnography, peripheral oxygen saturation (S p O 2 ), and electrocardiogram waveforms were obtained using a Philips patient monitor equipped with the M3001A measurement module (Philips Healthcare, Andover, MA). Surgical incision into the neck was performed to allow cannulation of the internal jugular vein and carotid artery, as well as relocation of the endotracheal tube into an incision in the trachea just below the larynx. Endotracheal tubes were manually shortened (to final length 12-15 cm depending on tube diameter) prior to tracheostomy, to reduce apparatus deadspace. Anesthesia was maintained using an intravenous infusion of propofol (7-9 mg kg −1 hr −1 ), and muscular relaxation was achieved using intermittent intravenous boluses of either rocuronium (1-2 mg kg −1 ) or pancuronium (0.01-0.15 mg kg −1 ).
Lung injury was induced by slow infusion of 0.08 cm 3 kg −1 oleic acid into the internal jugular vein over 15 min. Maturation of injury required ∼90 min, and was confirmed by a ratio of arterial oxygen tension to inspired oxygen fraction (P a O 2 : F i O 2 ) <300 mmHg during ventilation with 5 cmH 2 O of positive end-expiratory pressure (PEEP). Additional oleic acid (0.04 cm 3 kg −1 ) was administered if P a O 2 : F i O 2 remained >300 mmHg 90 min after the initial dose. If necessary, normal systemic arterial blood pressure (systolic/diastolic ≥ 90/60 mmHg) was maintained with intravenous crystalloid solutions and intermittent doses of phenylephrine (1-2 µg kg −1 ).

Experimental Protocol
All measurements and ventilation modalities were performed under baseline conditions and following lung injury. Each subject was mechanically ventilated with CMV, HFOV, and MFOV in random order over 30-min intervals, using a FabianHFO hybrid oscillator/ventilator (Acutronic Medical Systems AG, Switzerland). Mean airway pressure (P aw ) was set to 12 cmH 2 O for all modalities. During baseline conditions, F i O 2 was set to 40%, but was increased to maintain S p O 2 ≥ 90% following lung injury. Pressure-controlled CMV was delivered at rates between 20 and 32 min −1 (0.33 to 0.53 Hz), with inspiratory:expiratory ratio of 1:2. Single-frequency HFOV was delivered at 5 Hz, while MFOV was delivered using a combination of 5, 10, 15, and 20 Hz superimposed sinusoids with uniform flow amplitudes (Kaczka et al., 2015). Example ventilator waveforms are provided in Supplementary Figures S-1, S-2, S-3. Sampling frequency for recorded ventilator waveforms was 200 Hz. Ventilator driving pressure or pressure amplitudes were adjusted to obtain arterial CO 2 tension (P a CO 2 ) in the target range of 30-60 mmHg. Each 30-min experimental interval was followed by an arterial blood gas analysis and 4DCT scan sequence, without interrupting mechanical ventilation . Between each experimental ventilation interval, a 15-min wash-out period of CMV and a 30-s recruitment maneuver to 35 cmH 2 O of airway pressure were used to restore a control mechanical and physiological state. After completion of the experimental protocol, subjects were euthanized with an intravenous solution of pentobarbital sodium and phenytoin sodium (1 mL + 0.2 mL kg −1 ).

Assessment of Gas Exchange and Mechanics
Oxygenation was assessed using the oxygenation index (OI), defined as (Ortiz et al., 1987): The efficiency of CO 2 elimination during each modality was assessed using a ventilatory cost function (V C ), defined as (Kaczka et al., 2015): where V rms is the root mean-square of the volume waveform V (t) delivered to the airway opening: T denotes the time duration over which the integration in Equation (3) is performed (i.e., one breath for CMV or one second during HFOV/MFOV), and V is the mean of the volume waveform over the interval 0 to T. 1 Peak-to-peak range of volume (V pp ) was computed as the difference between the largest and smallest values of time-varying volume: Dynamic respiratory system elastance (E rs ) during pressurecontrolled CMV was computed as the quotient of driving pressure ( P aw ) and tidal volume: where: Measurements of respiratory impedance (Z rs ) were obtained under baseline conditions and after maturation of lung injury. A pseudorandom waveform consisting of nine sinusoids ranging in frequency from 0.078 to 8.9 Hz was generated by a digital-toanalog converter (NI USB-6356, National Instruments, Austin, TX), low-pass filtered at 12 Hz (858L8B-2, Frequency Devices, Haverhill, MA), and used as the input driving signal to a custom servo-controlled pneumatic pressure oscillator (Kaczka and Lutchen, 2004). Airway pressure was measured with a transducer placed at the proximal end of the endotracheal tube (Celesco LCVR-0050, Canoga Park, CA), while flow was measured using a screen pneumotach (4700A, Hans Rudolph, Shawnee, KS) coupled to a differential pressure transducer (Celesco LCVR-0002, Canoga Park, CA). The airway pressure and flow signals were sampled at 40 Hz over ∼90 s during the application of forced oscillations, allowing for three repetitions of the 25.6-s periodic control signal. Each Z rs spectrum and its corresponding coherence function (γ 2 ) was computed using the Welch periodogram technique (Welch, 1967), with rectangular windowing and 80% overlap. The complex values of Z rs were evaluated only at the nine frequencies used in the control signal (Suki and Lutchen, 1992). Corresponding resistance R rs and reactance X rs spectra were defined by the real and imaginary parts of Z rs , respectively. The resonant frequency f res was estimated from the zero-crossing of X rs using spline interpolation. Each Z rs spectrum was then characterized by a 4-element constant-phase model consisting of parameters for Newtonian airway and chest wall resistance (R), airway inertance (I aw ), tissue hysteresivity (η), and tissue elastance (H): whereẐ rs denotes the model-predicted Z rs and α = 2 π tan −1 1 η . Model parameters were estimated by constrained nonlinear gradient descent technique (MATLAB, Natick, MA).

Image Acquisition and Processing
CT scans were acquired using a Siemens SOMATOM Force (Siemens Healthineers, Germany) in an axial scanning mode, with 5.76 cm of axial coverage and 0.6 mm slice thickness. Subjects were continuously scanned at 80 kVp tube voltage and 150 mA tube current, with 250 ms scanner rotation period. Each scan lasted a total duration of 30 s, resulting in a radiation exposure of 345 mGy and generating a continuous set of projection images. The 4DCT image sequences were FIGURE 1 | Schematic illustrating the jointly estimated 3D transformations φ n (X) = x n mapping between all images in the periodic 4DCT image sequence I n (xn) and the target average image I avg (X). Solid black arrows indicate the transformation from the target reference frame X to the general reference frame x n . Dashed black arrows indicate the inverse transformation. Gray arrows represent the periodic sequence of N ventilatory phases indexed from 0 to N − 1. Any I n (xn) may be deformed to achieve structural alignment with any other I m (xm) by applying a sequence of two transformations x n = φ n φ −1 m (xm) , warping x m through X.
reconstructed by retrospective binning of x-ray projection data according to mechanical ventilation phase using a frequency-selective reconstruction algorithm , yielding between 13 and 21 volumetric images per 4DCT sequence with isotropic 0.6 mm spatial resolution. The corresponding temporal sampling frequencies for 21-phase image sequences were 7 Hz during CMV with fundamental 0.33 Hz, and 105 Hz during HFOV or MFOV with fundamental 5 Hz. Each sequence was periodic in the temporal (i.e., ventilation phase) dimension, such that the choice of "initial" image in the sequence was arbitrary. Each volumetric image in the temporal sequence was labeled I n (x n ), where n indexes the number of images in the sequence 0 through N − 1, and x n is a vector representing 3-dimensional spatial position. Voxels corresponding to spatial positions within the lungs were identified by a fully automated segmentation algorithm using a deep convolutional neural network (Gerard et al., 2018(Gerard et al., , 2020, generating a distinct lung mask M n (x n ) for each image phase. The neural network was trained using manually segmented lungs in CT images obtained from multiple datasets of experimental lung injury models, including a subset of images from the current study. The periodic motion of respiratory structures was estimated using a deformable 4D image registration technique. This procedure produced N 3-dimensional transformation functions denoted as φ n (X) = x n , where X is a vector representing 3dimensional spatial position within the reference frame of a phase-averaged image I avg (X), and x n represents spatial position in image I n (x n ). Image I avg (X) was used as the target for groupwise image registration of the images I n (x n ), n = 0 . . . N − 1. Thus, φ n (X) = x n is a spatial mapping from the coordinate system of I avg (X) to that of I n (x n ), such that a deformed image I n→avg (X) was generated with respiratory structures aligned to the target image as: Figure 1 shows a schematic of the transformations mapping between all images in the periodic sequence and I avg (X).
The transformations φ n (X) were coupled together using 4dimensional cubic B-splines and jointly estimated to ensure smoothness across both spatial and temporal dimensions (Metz et al., 2011). The parameters of the transformation function were iteratively adjusted using the Elastix library (Klein et al., 2010) to minimize the sum of squared tissue volume differences (SSTVD) between each I n→avg (X) and I avg (X), such that each deformed image preserved the volume of tissue contained within each voxel (Gorbunova et al., 2008;Yin et al., 2009;Zhao et al., 2016). Thus, changes in CT voxel density due to variations in fractional gas content were adjusted using this similarity cost function. Following registration, φ n (X) and its inverse mapping φ −1 n (x n ) were used to align images between arbitrary phases. Deformed images were denoted as I n→m (x m ), indicating the nth image deformed into the mth spatial reference frame: For the trivial case of m = n, Equation (9) simplifies to I n→n (x n ) = I n (x n ). Conceptually, aligning images from any two phases can be achieved by warping one spatial reference frame to another, passing through the target reference frame X (see Figure 1). Using this approach, all images were warped to align structures with a single, arbitrarily selected reference phase at m = 0. Phasic variations in regional aeration and strain could then be associated with the tissue contained within the region of interest at the reference phase. Figure 2 schematizes the use of a single reference phase to align respiratory structures across several ventilatory phases, thus tracking tissue-correlated changes in aeration and regional volume. Changes in regional aeration were then assessed by the range of voxel intensity values in the deformed image: while changes in regional volumetric strain were estimated using the Jacobian matrix (J n→m ) of transformation spatial derivatives: where ∇ x m is the spatial gradient operator. Note that J n→m is expressed in the mth spatial reference frame and describes the gradient of the "pullback" transformation φ n φ −1 m (x m ) . This pullback transformation is interpreted as "pulling back" point x n in ventilation phase n to the corresponding point x m in phase m. Regional volume changes relative to the reference phase (i.e., m = 0) were computed by the determinant of J n→ m : where V n (x 0 ) corresponds to the phase-varying volume of a region which, in the reference phase 0, occupies a single voxel centered at position x 0 . Accordingly, V 0 (x 0 ) is everywhere equal to the volume of a single voxel δV, determined by the image spatial resolution. Thus, Equation (12) may be simplified to: Note that |J n→0 | < 1 when the corresponding region deflates relative to the reference phase, |J n→0 | > 1 when it expands, and |J n→0 | = 1 when there is no volume change (i.e., isovolumetric deformation). Figure 3 shows an example image sequence for each ventilation modality in one injured subject, demonstrating phase-varying aeration I n→0 (x 0 ) and volume change V n (x 0 ) for each voxel in a sample region of interest in x 0 . Animated image sequences for one subject under baseline and injured conditions are also provided (Supplementary Animations S-1, S-2). Phasevarying strain (ε) for each lung region with respect to its respective minimum inflation state was defined as: Changes in regional volumetric strain were assessed by the range of ε n (x 0 ): Since the minimum value of ε n (x 0 ) is zero by definition, Equation (15) may be reduced to: Similarly, regional volumetric strain rate (ε) was computed by the temporal rate of change in ε n (x) using the forward difference scheme: where δt is the time difference between adjacent phases. The modulo operator is used to indicate the periodicity of ventilatory phase. The value of δt is computed from the ventilatory fundamental frequency f 0 and the number of images in the 4DCT sequence (N): Finally, changes in regional volumetric strain rate were assessed by the range ofε n (x 0 ): Note thatε n (x 0 ) generally will be positive during inflation of a specified lung region, and negative during the corresponding deflation. Thus, ε (x 0 ) will reflect the sum of the fastest inflation and deflation rates of the region. Voxels within 6.0 mm of the image boundaries along the rostral-caudal axis were excluded from the mask for the purposes of image analysis, due to potential registration artifacts caused by motion of lung tissue into or out of the axial field of view. Regional aeration ( I), strain ( ε), and strain rate ( ε) were compared across subjects, lung conditions, and ventilation modalities by their mean value throughout the lung mask, FIGURE 2 | Schematic of a periodic 4DCT image sequence during conventional mechanical ventilation (A), with inspiratory and expiratory phases (B). Image registration is used to map each image in the sequence to a single reference phase (C), after which intratidal variations in aeration and regional volume may be associated to specific regions of lung tissue (D). coefficient of variation, and spatial gradients in the right-left, dorsal-ventral, and rostral-caudal directions. Spatial gradients for each variable were determined from the slopes computed by linear regression with respect to position along each of the principal anatomic axes. Spatial gradients were examined unnormalized and normalized by the spatial mean, in each case comparing both signed values of the gradient (i.e., positive and negative) and gradient magnitude (i.e., absolute value).
Given the harmonic nature of MFOV, the individual frequency components of I n , ε n , andε n were assessed using the Discrete Fourier Transform (DFT), to relate these periodic time-varying properties at each spatial position into harmonic amplitude and phase components. The regional heterogeneity of volumetric strain ( ε) was assessed using octree and supervoxel decomposition, each of which recursively subdivide regionsof-interest (ROIs) within the lung mask. Both decomposition methods were initialized with a single ROI containing the entire lung mask. At each recursive step, the designated ROI was subdivided into multiple new ROIs if the standard deviation of the contained voxel intensity was greater than a fixed FIGURE 3 | Four-dimensional image registration was used to deform images such that respiratory structures were aligned to a single reference phase. After registration, regional aeration and volume change throughout the lung were assessed by the variations in deformed image voxel value and determinant of the Jacobian matrix of the deformation. This technique was applied during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Ventilatory cycle duration is indicated by spacing between tick marks. In this example, each ventilatory cycle corresponds to a sequence of images representing 21 distinct ventilatory phases. Note that multiple cycles of each ventilation waveform are shown for clarity, however all cycles are identical replicates. See also Supplementary Animations S-1, S-2.
threshold, set to a standard deviation of mean-normalized ε ≥ 0.3. 2 Octree decomposition, a three-dimensional extension of quadtree decomposition, was used to partition each ROI into eight octants according to bisecting coronal, sagittal, and transverse planes (Perchiazzi et al., 2014). Supervoxel decomposition was used to partition each ROI into two new ROIs by a weighted k-means clustering of voxels according to intensity similarity and spatial proximity (Conze et al., 2017). Recursive subdivision was continued until all ROIs either contained belowthreshold values of standard deviation, or were <0.2% of the lung mask. Regional heterogeneity of ε was then quantified by the mean ROI volume (V ROI ) as a fraction of the lung mask.

Statistics
Non-parametric Kruskal-Wallis rank sum tests were performed for each outcome, testing for main effects of lung condition (i.e., baseline, injured) and ventilation modality (i.e., CMV, HFOV, MFOV) at the 0.05 significance level. The effect of ventilation modality was tested separately under baseline and injured conditions. For outcomes with a significant main effect of ventilation modality, differences among modalities were identified by a Dunn post hoc comparison with adjustment for rank ties and further adjustment by the Benjamini-Hochberg procedure for reduced false discovery rate. A web-based tool was used for statistical calculations (Vasavada, 2016). Figure 4 shows a summary of the respiratory impedance spectra Z rs across all subjects measured under baseline and injured conditions, while Table 1 provides a summary of the constant phase model applied to Z rs . Under baseline conditions, average model-based tissue elastance H was 121 cmH 2 O L −1 at 1 rad s −1 (0.16 Hz) with a coefficient of variation of 0.31. The increase in H after oleic acid injury was highly variable, resulting in a mean of 435 cmH 2 O L −1 at 1 rad s −1 with a coefficient of variation of 0.94. By comparison, measurement of dynamic elastance (E rs ) during CMV (Equation 5) was 148.9 ± 26.9 cmH 2 O L −1 and 257.5 ± 64.2 cmH 2 O L −1 under baseline and injured conditions, respectively. Resonant frequency (f res ) increased after lung injury, in some cases above the measurement range allowed by the sampling parameters. Table 2 contains a summary of all significant effects identified by non-parametric Kruskal-Wallis tests, as well as post hoc comparisons as appropriate. These results are described further in the context of specific figures below. FIGURE 4 | Respiratory impedance measured under baseline (black) and injured (white) conditions, represented by the mean value (circles) and standard deviation (error bars) at each frequency, across subjects. The complex-valued impedance spectrum is represented by its real part or in-phase component (i.e., resistance), and its imaginary part or out-of-phase component (i.e., reactance). Coherence provides an estimate of signal quality and linearity of the measured relationship between pressure and flow oscillations.  (5). Paw, mean airway pressure (cmH 2 O); R, Newtonian airway and chest wall resistance (cmH 2 O L -1 s); Iaw, airway inertance (cmH 2 O L -1 s 2 ); η, tissue hysteresivity; H, tissue elastance (cmH 2 O L -1 at 1 rad s -1 ); fres, resonant frequency (Hz); Ers, dynamic respiratory system elastance (cmH 2 O L −1 ). Figure 5 shows gas exchange outcomes related to mechanical ventilation, including oxygenation index (OI = P a O 2 −1 · F i O 2 · P aw ) and ventilatory cost function (V C = P a CO 2 · V 2 rms · Wt −1 ). A significant main effect of lung condition (baseline vs. injured) was found for P a O 2 , P a O 2 · F i O 2 −1 , OI, and P a CO 2 (p < 0.001). A significant main effect of ventilation modality was found in P aw , V pp · Wt −1 , V rms · Wt −1 , and V C (p < 0.001) for both lung conditions, as well as in OI for the injured condition (p < 0.05). In particular, median OI increased by 300% for injured subjects compared to baseline (p < 0.001), but was 48% lower for injured subjects during MFOV compared to CMV (p < 0.05). MFOV and HFOV were associated with lower median V C for both baseline and injured conditions compared to CMV (HFOV 84% lower than CMV, p < 0.01; MFOV 93% lower than CMV, p < 0.001). Pairwise comparison of V C for MFOV against HFOV also achieved a statistically significant difference, with 60% lower median V C for MFOV compared to HFOV (p < 0.05). Median V rms · Wt −1 was 41% lower for MFOV vs. HFOV (p = 0.0505 baseline; p < 0.05 injured) and 75% lower for MFOV vs. CMV (p < 0.001). Figure 6 shows registration-based estimates of regional ventilation, as intratidal changes in voxel aeration ( I), volumetric strain ( ε), and volumetric strain rate ( ε). A significant main effect of condition (baseline vs. injured) was found for I spatial mean (p < 0.001), ε coefficient of variation (p < 0.001) and ε coefficient of variation (p < 0.01). A significant main effect of ventilation modality was found for the spatial means of all three variables-I, ε, and ε-in both baseline and injured lung conditions (p < 0.001). In particular, MFOV and HFOV both produced smaller I, smaller ε, and larger ε compared to CMV under baseline and injured conditions (p < 0.05). For the injured condition, MFOV resulted  in significantly lower I and ε compared to HFOV (p < 0.05), but no difference for ε (p = 0.20). Figure 7 shows unnormalized spatial gradients of I, ε, and ε in the right-left, dorsal-ventral, and caudal-rostral axes. The largest spatial gradient magnitudes observed occurred along the dorsal-ventral axis. In particular, the injured lung condition was associated with significantly more positive dorsal-ventral gradients of ε (p < 0.001) and ε (p < 0.05). Ventilation modality had a significant effect on all dorsal-ventral gradients under baseline conditions ( I, p < 0.001; ε, p < 0.001; ε, p < 0.01), with MFOV producing the least negative dorsal-ventral gradients of I and ε, but the most positive dorsal-ventral gradients of ε. Although modality was not a significant predictor of variability in the actual value (positive or negative) of any dorsal-ventral gradient for injured subjects, it was a significant main effect for the magnitudes of all dorsalventral gradients ( I, p < 0.001; ε, p < 0.01; ε, p < 0.05). In general HFOV and MFOV reduced the dorsal-ventral gradient magnitudes of I and ε under both baseline and injured conditions, while increasing the corresponding FIGURE 5 | Gas exchange outcomes, represented by box plots showing minimum, maximum, and quartiles across subjects under baseline (black) and injured (white) conditions, during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Asterisks indicate significance levels (*p < 0.05, ***p < 0.001) for main effects (C = lung condition, MB = modality for baseline condition, MI = modality for injured condition). See Table 2 for additional clarification. P a O 2 = arterial oxygen tension; F i O 2 = fractional inspired oxygen; P a CO 2 = arterial carbon dioxide tension; P aw = airway pressure amplitude; V pp = peak-to-peak volume range; Wt = body weight; V rms = root-mean-square volume.

RESULTS
ε gradient. Figure 8 shows corresponding mean-normalized spatial gradients. In general, normalizing each spatial gradients by the respective spatial mean accounted for much of the variability observed in the unnormalized spatial gradients (Figure 7) with respect to varying ventilation modality. This is also evidenced by the loss of significant differences between ventilation modalities for the corresponding spatial gradients as indicated in Table 2. Figure 9 shows an example of regions-of-interest defined by recursive supervoxel and octree decomposition performed on mean-normalized ε in one subject. Clustered regions are illustrated using contours and filled regions in two-dimensional FIGURE 6 | Summary of image registration-based regional ventilation measures for intratidal changes in aeration ( I), volumetric strain ( ε), and volumetric strain rate ( ε). Mean and coefficient of variation throughout the lung mask are shown for each variable, represented by box plots showing minimum, maximum, and quartiles across subjects under baseline (black) and injured (white) conditions, during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Asterisks indicate significance levels (**p < 0.01, ***p < 0.001) for main effects (C = lung condition, MB = modality for baseline condition, MI = modality for injured condition). See Table 2 for additional clarification.
views of a slice, as well as translucent surface renderings in threedimensional projections. Note that supervoxel decomposition produces larger contiguous ROIs, despite using the same decision criteria for recursive cluster subdivision. Figure 9 also provides a summary of mean cluster size (V ROI ) across subjects, condition, and ventilation modality. A significant effect of lung condition was found for the octree technique (p < 0.001), with 32% smaller median V ROI in injured lungs. Figure 10 shows example distributions of strain amplitude and phase at the first four harmonics for CMV, HFOV, or MFOV, computed using the discrete Fourier transform of time-varying Jacobian determinants in aligned images. Figure 11 provides a summary of the harmonic strain amplitudes during HFOV and MFOV across all subjects, correlated with corresponding amplitudes of the ventilator volume waveform V (t). The oscillatory waveforms were a 5 Hz sinusoid (HFOV) and a combination of 5, 10, 15, and 20 Hz sinusoids with uniform flow amplitudes (MFOV). The relative amplitude distributions in the flow waveforms delivered by the ventilator were largely preserved throughout the lung, according to the registration-based metrics (r 2 = 0.78). However, harmonic distortion measurements during HFOV indicated nearly 20% of the spectral power in regional strain was concentrated in the higher-order harmonics of the 5 Hz waveform (i.e., 10, 15, 20 Hz), whereas the ventilator volume waveform exhibited only 5% harmonic distortion (Zhang et al., 1995;Amini et al., 2017). FIGURE 7 | Unnormalized spatial gradient values (A) and magnitudes (B) for image registration-based regional ventilation measures of intratidal changes in aeration ( I), volumetric strain ( ε), and volumetric strain rate ( ε). Linear spatial gradients in each of the three principal anatomic directions throughout the lung mask are shown for each variable, represented by box plots showing minimum, maximum, and quartiles across subjects under baseline (black) and injured (white) conditions, during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Asterisks indicate significance levels (*p < 0.05, **p < 0.01, ***p < 0.001) for main effects (C = lung condition, MB = modality for baseline condition, MI = modality for injured condition). See Table 2 for additional clarification.

DISCUSSION
In this study, we demonstrated that oscillatory ventilation improves the average regional lung strain, as well as the spatial gradients of lung strain, compared to a conventional pressure-controlled modality in pigs with heterogeneous lung injury. This study revealed a mechanism by which enhanced harmonic content in MFOV waveforms may improve gas exchange compared to either CMV or single-frequency HFOV waveforms.
The increase in respiratory system elastance following lung injury (i.e., 183% increase in E rs and 355% increase in H) is consistent with substantial lung derecruitment and/or surfactant dysfunction. Consequently, increased driving pressure was required to maintain eucapnia at comparable respiratory rates. Despite similar V pp · Wt −1 and ventilatory cost function (V C ), the driving pressure ( P aw ) after lung injury increased to 18.7 cmH 2 O, compared to 13.0 cmH 2 O at baseline (Figure 5). Increased P aw during CMV after injury exceeded the 15 cmH 2 O threshold identified by Amato et al. (2015), above FIGURE 8 | Mean-normalized spatial gradient values (A) and magnitudes (B) for image registration-based regional ventilation measures of intratidal changes in aeration ( I), volumetric strain ( ε), and volumetric strain rate ( ε). Linear spatial gradients in each of the three principal anatomic directions throughout the lung mask are shown for each variable, represented by box plots showing minimum, maximum, and quartiles across subjects under baseline (black) and injured (white) conditions, during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Asterisks indicate significance levels (*p < 0.05, **p < 0.01, ***p < 0.001) for main effects (C = lung condition, MB = modality for baseline condition, MI = modality for injured condition). See Table 2 for additional clarification.
which the mortality odds ratio was greater than unity in patients with ARDS. The use of P aw > 15 cmH 2 O suggests substantial potential for VILI despite average measured V pp · Wt −1 of 7.5 mL kg −1 (Figure 5). However, such heuristics are only applicable during CMV, for which P aw fluctuations are well-correlated with lung strain. During oscillatory modes of ventilation, P aw as computed from Equation (6) does not represent elastic pressures distending peripheral alveoli, but also includes substantial pressure losses due to resistive and inertial loads imposed by the endotracheal tube and conducting airways, especially with increasing frequency. The fact that P aw was greater during MFOV than HFOV indicates that larger pressure amplitudes were required to overcome the increased resistive and inertial losses at frequencies up to 20 Hz compared to 5 Hz, but does not indicate that these larger airway pressures were transmitted to distal alveoli (Pillow et al., 2002). We found that the spatial mean of strain throughout the entire lung during CMV was unchanged after injury, despite increased P aw and FIGURE 9 | Regional patchiness of intratidal variations in volumetric strain ( ε) during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). Regional patchiness was assessed by the size, number, and distribution of regions-of-interest or clusters (randomized colors) obtained by recursive supervoxel or octree decomposition. Examples from a representative subject (A), and cluster volumes as a fraction of the lung mask (B) represented by box plots showing minimum, maximum, and quartiles across subjects under baseline (black) and injured (white) conditions. Asterisks indicate significance levels (***p < 0.001) for main effects (C = lung condition).
I spatial mean (Figure 6). Some evidence to reconcile this discrepancy is apparent in the greatly increased coefficient of variation for ε (Figure 6). Upon closer examination of the regional strain distributions in injured lungs, ε was nearly zero in consolidated regions of the lung with I >-100 HU. By contrast, ε was increased by 25% in normally aerated regions (i.e., −900 HU < I < −500 HU) relative to the average throughout the entire lung mask. Thus, the increased P aw required to deliver similar V pp · Wt −1 (Figure 5) resulted in concentrated strains within the recruited lung. This finding is consistent with the concept of the "baby lung, " which describes reduced functional recruited volume in ARDS (Gattinoni and Pesenti, 2005;Gattinoni et al., 2016a).
Spatial gradients of intratidal strain also provide some insight into the redistribution of mechanical flows after lung injury. As expected, the magnitudes of right-left gradients were small relative to those of the dorsal-ventral and caudal-rostral gradients (Figure 7). Most noticeably, lung injury was associated with dramatic changes in dorsal-ventral ε gradients. Ventilation in healthy subjects during CMV tended to distribute toward the dorsal and caudal lung regions. However, varying degrees of lung injury resulted in consolidation or derecruitment of the dependent lung, such that dorsal-ventral ε gradients were reduced (or even reversed) depending on the extent of edema, particularly for CMV and HFOV. Interestingly, intratidal changes in aeration tended to preserve the dorsalventral I gradient before and after injury, despite the positive shift observed in dorsal-ventral ε gradients. It should be noted that actual spatial distributions of I and ε are likely not well-described by linear functions of spatial position following lung injury, and may even be non-monotonic (Johnson et al., 2017). For example, negative dorsal-ventral gradients may be possible in the recruited lung, but counteracted by zero ventilation in the derecruited dependent lung such that the resulting overall linear gradient may be negative, zero, or even positive. Nonetheless compared to either CMV of HFOV, MFOV reduced the magnitude of spatial gradients in ventilation regardless of the lung condition. Other strategies for improving ventilation distribution and eliminating large dorsalventral gradients include prone position ventilation (Scholten et al., 2017;Xin et al., 2018), which was not used in this study. Prone posture serves to eliminate large-scale gravitational gradients in lung expansion (Hoffman, 1985), but would not be expected to affect the local, region-to-region mechanical differences which MFOV targets. Finally, both HFOV and MFOV were associated with more positive caudal-rostral I gradients compared to CMV, regardless of injured condition (p < 0.001). CMV-associated findings indicate that in contrast to oscillatory modalities, CMV preferentially ventilates the basal regions of the lungs, consistent with a previous positron emission tomography study (Venegas et al., 1993).
Much of the variability in spatial gradient magnitudes across different ventilation modalities was due to changes in the spatial means. By contrast, the normalized gradients (i.e., the spatial gradients divided by the spatial mean) exhibited fewer significant differences across modalities (Figure 8, Table 2). This finding indicates that, at least for the particular frequencies and waveforms used in this study, the relative spatial gradients were not substantially altered by modality. This may be due in part to the predominance of the 5 Hz component in both HFOV and MFOV waveforms, which was below the resonant frequency in almost every case (Table 1), and therefore FIGURE 10 | Frequency-domain representations of tissue-correlated regional strain amplitude (A) and phase variations (B) in one subject during conventional mechanical ventilation (CMV), high-frequency oscillatory ventilation (HFOV), and multi-frequency oscillatory ventilation (MFOV). For each ventilation modality, strain amplitudes and phases were obtained from the discrete Fourier transform of time-varying Jacobian determinants within each voxel, using transformations mapped to a single reference image. In each cross-section of the reference image, voxel color in (A) indicates the strain amplitude at that frequency, corresponding to the amplitude of volumetric strain relative to the minimum volume of that particular voxel. Voxel color in (B) indicates the amount of phase shift (leading or lagging) relative to the median phase of all lung voxels at the specified frequency. below a threshold frequency for developing regional ventilation heterogeneity due to resonant amplification (Herrmann et al., 2019a). Our results then indicate that broadband MFOV waveforms, while producing the same relative ventilation heterogeneity as HFOV, still reduce the average strain and magnitude of strain gradients that are associated with risk for VILI.
Based on our supervoxel and octree analysis, regional strain heterogeneity did not differ greatly across ventilation modality, although the injured lungs did require smaller sizes to describe the spatial clustering of ε compared to the baseline conditions (Figure 9). Lung injury may thus result in increased smallscale strain heterogeneity (Kaczka et al., 2011;Perchiazzi et al., 2014). It should be noted that our estimates of regional strain FIGURE 11 | Harmonic amplitudes (A) from imaging-based regional strain measurements and corresponding amplitudes from the ventilator volume waveform measured at the airway opening under baseline (circles) and injured (triangles) conditions. Harmonic distortion index (B) mean and standard deviation for each set of measurements, grouped by baseline (black) and injured (white) lung conditions. derived from image registration are inherently smooth due to the use of B-splines to define the spatiotemporal deformations. The discrepancy between the supervoxel-and octree-based results may also be due in part to the discrete anatomic structure of lungs, which consist of an assumed continuum of spongy parenchymal tissues penetrated by branching networks of airways and vessels. Airways and vessels are not expected to undergo large volumetric strains during mechanical ventilation, and therefore give rise to irregularities in the spatial field of ε. Such irregularities may be approximately contoured by supervoxel boundaries, which are not rigidly constrained. By contrast, octree decomposition is strictly planar, and therefore may require greater degrees of recursive subdivision to isolate such irregular spatial distributions of ε (Perchiazzi et al., 2014). The supervoxel decomposition incorporates spatial proximity as well as image intensity, such that reducing the spatial proximity weighting may facilitate the identification of larger clusters albeit with irregular boundaries. Therefore, it is possible that our findings would be altered if smaller airways and blood vessels were further segmented from the lung parenchyma prior to cluster analysis, or if a different spatial proximity weighting was used for the supervoxel decomposition.
Separation of the harmonic content of regional ventilation illustrates a potential influence of nonlinear mechanics during HFOV at 5 Hz (Figure 11). For example, image registration resulted in estimation ofε amplitudes at 10, 15, and 20 Hz of roughly one third those at 5 Hz, despite the ventilator waveform consisting of only 5 Hz oscillation (Supplementary Figure S-4). By contrast,ε amplitudes during MFOV approximated the input distribution of uniform amplitude flow at all four harmonics (Supplementary Figure S-4).
Intrapulmonary harmonic distortion of distributed flows during HFOV (with only a single frequency) may result from the use of larger-amplitude flows exacerbating inherent mechanical nonlinearities, but may also result from the complex interaction of mechanicallyinterdependent lung tissues (Suki and Bates, 1991;. During MFOV, the use of small-amplitude flows may contribute to relatively linear mechanical behavior of the lung tissues. It is also possible that the observed nonlinearity may be due in part to imaging artifact and errors in image registration resulting in misidentified high-frequency, small-amplitude motion. Assuming that our registration results are accurate however, theε amplitudes measured at the higher harmonics during HFOV are surprisingly large, almost comparable to those intentionally delivered during MFOV. It is plausible that the mechanisms by which MFOV may enhance gas transport already occur, at least to some extent, within the lung periphery during HFOV. However, MFOV may induce these higher harmonics directly, without incurring the overhead of generating high flow rates at the fundamental frequency. The consequences of such overhead may include, for example, increased spatial gradients of strain and strain rate due to gravitational dependence (Figure 7).

Limitations
MFOV may improve ventilation distribution and the efficiency gas exchange efficacy (Kaczka et al., 2015;Herrmann et al., 2018Herrmann et al., , 2019b. However, given the technical limitations of the mechanical ventilator as well as our dynamic CT image reconstruction technique which assumed strict periodicity of motion , the MFOV waveform used in this study comprised only harmonic components of the fundamental frequency. Cardiogenic motion was not synchronized with respiratory motion and therefore presents a source of artifact in the reconstructed 4DCT images of respiratory motion, depending on the relative periodicity between the ventilator frequency and the heart rate. Such artifact may manifest as rapid oscillations in CT density , which may influence registration-based estimates of strain rate, especially in regions adjacent to the mediastinum. In addition, the interpretation of regional strain in the proximity of recruitment/derecruitment is challenging for several reasons. Standard image registration techniques (e.g., using a similarity cost function defined as the sum of squared intensity differences) assume a uniform intensity transform to predict voxel-wise correspondence, but cannot account for local changes in intensity associated with changing gas fraction. Consequently, these approaches may result in the prediction of non-physical deformations, such as contraction of nonaerated lung regions. The sum of squared tissue volume differences (SSTVD) similarity cost function used in this study addresses this problem using a non-uniform intensity transform that accounts for local variations in CT intensity that are due to changing proportions of gas-to-tissue content in each lung region during breathing or ventilation (Zhao et al., 2016). Another challenge to registration of injured lungs is the lack of contrast or structural landmarks in consolidated regions, which renders validation of estimated voxel-wise correspondence nearly impossible, especially when large changes in recruitment occur between two images. In this study, high temporal resolution imaging (7-105 Hz) was used to minimize the time between adjacent ventilation phases, which resulted in small deformations and small changes in recruitment between adjacent images being registered. Linear elastic regularization produces plausible predictions of voxelwise correspondence in regions with little contrast, assuming the lung parenchyma deforms as a linear elastic material. An important limitation of most registration techniques is the assumption of smoothly varying deformation, which may not accurately reflect discontinuous motion at distinct structural boundaries between aerated and nonaerated lung tissue. Due to the limited spatial resolution in this study (0.6 mm), each voxel may contain up to dozens of alveoli. Thus, it is difficult, if not impossible, to distinguish the underlying cause of observed changes in aeration (Cereda et al., 2019). For example, the same change in aeration of a single voxel may represent uniform deflation of its alveoli, complete derecruitment of only a portion of alveoli, or any combination thereof. Thus, prediction of recruitment/derecruitment on the basis of voxel aeration is speculative.
Finally, while our study focused on regional distributions of aeration, strain, and strain rate as the primary indicators of ventilation and VILI, the exact relationship between these mechanical variables and VILI is difficult to ascertain. VILI is a complicated process associated with a diverse range of mechanical stimuli (Hussein et al., 2013;Carrasco Loza et al., 2015;Güldner et al., 2016;Smith et al., 2017;Tonetti et al., 2017). So-called "mechanical power" (i.e., the rate of energy dissipation across parenchymal tissues) may provide a better prediction of VILI compared to either strain or strain rate alone (Gattinoni et al., 2016b;Serpa Neto et al., 2018). However, measurement of regional mechanical power is not practical, given the requirement of regional gas pressure within the lung for its estimation. Alternatively, VILI may also be assessed using fluoro-deoxyglucose (FDG) uptake (Wellman et al., 2014), serum or alveolar inflammatory cytokines (Liu et al., 2013), protein content in bronchoalveolar lavage fluid (Smith et al., 2017), and histopathology. All of these quantitative measurements of lung injury require extensive durations of mechanical ventilation prior to measurement, some of which may only be obtained post mortem.

Conclusion
In an oleic acid model of porcine ARDS, high-frequency and multi-frequency oscillatory ventilation resulted in improved gas exchange efficiency and reduced regional lung strain compared to conventional mechanical ventilation. MFOV also reduced spatial gradients in lung strain compared to CMV or HFOV, resulting in a more uniform overall distribution of lung strain. By contrast, mean-normalized spatial gradients exhibited less dependence on ventilation modality, indicating that the use of higher harmonic frequencies in MFOV does not substantially alter relative ventilation heterogeneity. Thus, the reduced spatial strain gradients observed with MFOV, compared to HFOV or CMV, may be explained by reductions in the average parenchymal strain. In summary, MFOV supports gas exchange with reduced lung strain, and may provide additional benefit over HFOV by reducing ventilation heterogeneity associated with spatial gradients in absolute strain, and by reducing the strain associated with conventional ventilation or single-frequency oscillation.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The animal study was reviewed and approved by University of Iowa Institute for Animal Care and Use Committee (Protocol Number 5061428).

AUTHOR CONTRIBUTIONS
JH and DK conceived the study, collected the data, and wrote the manuscript. JH, SG, WS, MH, JR, GC, EH, and DK analyzed the data and revised the manuscript.