**Magnetic Resonance Elastography**

Six normal adult volunteers n = 6 (age 24–29 yrs, 3 male, 3 female) and three pwCF (1 male and 2 females; 21–35 years old) were scanned after obtaining written informed consent. The study was approved by the Institutional Review Board (IRB) at both the Ohio State University and Nationwide Children’s Hospital. All imaging was performed on a 1.5 T MRI scanner (Avanto, Siemens Healthcare, Erlangen, Germany). All subjects were positioned supine and headfirst into the scanner. For MRE scans, mechanical vibrations at a frequency of 50 Hz were introduced separately into the right and left lungs using an in-house built acoustic driver system. The passive driver was positioned anteriorly on the apex of the lungs and the passive driver was connected to the active driver (speaker system) using a plastic tube placed outside the scan room.

**MRI and MRE Image Acquisition**: A modified spin echo-echo planar imaging (SE-EPI) sequence was used to acquire five axial slices of the lung. To minimize T2* effects and to achieve an adequate MR signal in the lungs, a SE based sequence was used. To achieve a strong signal in the lung, the shortest achievable TE of 11.6 ms was used. The sequence parameters used were: FOV: 40x40cm2; TR: 400ms; acquisition matrix: 128x64 interpolated to 256x256; slice thickness: 10mm; voxel size 1.56x1.56x10mm; echo train length: 9 (i.e. 7 shots to fully cover a single k-space); and 4 MRE phase offsets. For MRE, a 1:1 motion encoding gradient (MEG) was split into two unipolar lobes around the 180° refocusing pulse to achieve the minimum possible TE. Each Motion Encoding Gradient (MEG) lobe had a period of 2 ms (i.e. 250 Hz combined). In-plane and through plane mechanical motion were encoded with a 26 second breath hold for each subject at residual volume (RV).

**Lung Density Imaging**: The lung’s density (LD) changes during respiratory cycle as described previously(Holverda et al. 2011; Theilmann et al. 2009; Mariappan et al. 2011; Mariappan et al. 2014; Fakhouri, Dong, and Kolipaka 2019) and quantifying lung density is crucial for accurately estimating lung stiffness. Therefore, LD scans were performed using a fast GRE sequence(Holverda et al. 2011; Theilmann et al. 2009) with TR of 10 ms. To calculate T2* at RV, four different TEs of 1.07 (minimum TE achievable by the sequence), 1.5, 2, and 2.5ms were selected. Due to the low density of the lung at TLC relative to RV, the two shortest TEs were selected (i.e. 1.07 and 1.5ms) because of the rapid decay of the signal at TLC(Holverda et al. 2011; Theilmann et al. 2009). This scan involved a single breathhold for each, RV and TLC. A whole-body coil was used with the following parameters: FOV: 50x50cm2, slice thickness: 10mm, acquisition matrix: 64x64, and number of averages: 4.

**Calculation of Shear Stiffness**

The left and right lung images were first segmented by drawing a region of interest (ROI) in each of the 5 axial slices that defined the lung outline. Next, to eliminate longitudinal and reflected waves, a 4th order Butterworth bandpass directional filter was applied in 8 directions with cutoff values of 42 − 40 waves/FOV at RV(Manduca et al. 2003). Filter cutoff values were selected based on the number of pixels manually measured in a wavelength. Third, lung shear stiffness was calculated for each slice individually by using 2D direct inversion(Manduca et al. 2001) using MRElab (Mayo Clinic, Rochester, Minnesota, USA). Finally, the resultant stiffness maps were median filtered (3x3 kernel)(Huang, Yang, and Tang 1979; Frieden 1976), and the 95th percentile(Langford 2006) of the data was used by removing outliers due to noise imbedded in the lung’s stiffness maps to report the mean stiffness value.

**Calculation of Lung Density**

Lung density (LD) was estimated relative to a Gadolinium-doped water phantom that was placed on the volunteers’ chest during LD scans(Holverda et al. 2011; Theilmann et al. 2009). The MR signal received from the phantom, which is mostly water, was considered as a good representation of lung’s signal without any air. Therefore, an absolute measure of lung’s density can be estimated relative to the Gadolinium-doped water phantom. The following equation was used to obtain the initial signal (I0) of the lung for each pixel in the ROI,

$$Sm={I}_{0}{e}^{TE*T{2}^{*}}$$

1

where Sm is the mean signal of the lungs at a given TE (i.e. 1.07, 1.5, 2, 2.5 ms). The resultant initial signal of the lung (I0) was then used to calculated LD in each pixel by

$$LD={I}_{0}*CR*{I}_{ph}$$

2

where Iph is the mean signal of the Gadolinium-doped water phantom, and CR is a correction factor that had to be determined due to the long decay constants of the phantom and to correct for the steady state signal that the phantom might reach at a TR of 10ms. By using the same fast GRE sequence, the CR was obtained by scanning the phantom twice with two different TRs of 10ms and 6 seconds. Then the CR was obtained by dividing mean phantom signal at TR = 6 seconds by mean phantom signal at TR = 10ms which resulted in a CR of 1.873(Holverda et al. 2011; Theilmann et al. 2009; Mariappan et al. 2014). LD was used to correct the shear modulus(Mariappan et al. 2011) and the average LD in each subject was used in the FEM models described below.

**Finite Element Modeling**

3D FEM of lung deformation during normal tidal volume breathing was developed by inputting MRE-derived shear stiffness values at RV obtained in normal (n = 6) and cystic fibrosis (CF) (n = 3) subjects. Notably, both lungs were scanned in pwCF, and only the right lung was scanned in healthy adults. For each normal and CF lung, an outline of the lung in each axial MR scan was obtained using ImageJ and a custom written MatLab code that uses non-uniform rational basis splines (NURBS)(Spink) was used to create a 3D geometry of the lung section scanned during the MRI protocol (Fig. 1A). Due to edge-effects in the MRE-derived stiffness maps, the size of the MRE shear stiffness maps were cropped and smaller than the actual lung geometry outlines. As a result, separate inner outlines were drawn using the MRE measurements (Fig. 1B) and used to create an inner 3D geometry using NURBS. The outer and inner 3D geometries were imported into the COMSOL multi-physics finite element package and converted into a 3D solid domain (Fig. 1C). The outer geometry was used for specifying boundary conditions (Fig. 1D) while the inner geometry allowed for direct specification of spatial distributions in shear stiffness (G, referred as shear modulus in FEM). Specifically, the MRE measurements of LD corrected shear stiffness obtained in each plane (Fig. 1B) were interpolated using a 3D linear interpolation function and this function was used to directly specify shear modulus in the inner region. For these studies, the shear modulus in the outer region was assumed to be equal to the average shear modulus within the entire model. The resulting model was then meshed with quadratic tetrahedral elements (Fig. 1C) and the Solid Mechanics module in COMSOL was used to simulate lung tissue deformation during normal tidal volume respiration by specifying appropriate material models and boundary conditions.

The time-dependent equations governing tissue deformation solved in COMSOL were

$$\begin{array}{cc}\begin{array}{cc}{\rho }\frac{{\partial }^{2}\varvec{u}}{\partial {t}^{2}}={\nabla \bullet \left(FS\right)}^{T}& F=I+\nabla \varvec{u}\end{array}& \begin{array}{cc}S=\frac{\partial W}{\partial ϵ}& \epsilon =\frac{1}{2}\left({F}^{T}F-I\right)\end{array}\end{array}$$

3

where ρ is tissue density, **u** is the displacement field, *F* is the deformation gradient, *I* is the identity tensor, *S* is the second Piola-Kirchhoff stress tensor, *W* is the strain energy density and ε is the strain tensor. Previous studies indicated that finite element models of lung deformation that utilize a two-parameter Mooney-Rivlin hyperelastic material model yielded the highest accuracy in capturing experimentally measured motion of lung tumors (Tehrani et al. 2015). Therefore, in this study we implemented this material model by specifying

$$\begin{array}{ccc}W={c}_{1}\left({I}_{1}-3\right)+{c}_{2}\left({I}_{2}-3\right)+{\kappa \left(J-1\right)}^{2}& {c}_{1}={c}_{2}=\frac{G}{4}& \kappa =\frac{2G(1+v)}{3(1-2v)}\end{array}$$

4

Here, *I**1* and *I**2* are the first and second invariants of the isochoric elastic right Cauchy-Green tensor, κ is the bulk modulus and *J* is the elastic volume ratio. The shear modulus, *G*, was either specified as a constant value or specified based on MRE-derived measurements as described above. Since the lung has been modeled with a range of Poisson’s ratios, 0.2–0.5 (Zhang et al. 2004; Tehrani et al. 2015; Villard et al. 2005; Werner et al. 2009), we conducted a sensitivity analysis using the input values from Table 1 to determine the influence of Poisson’s ratio on strain and strain gradient magnitudes. Results (Supplemental Fig. 1A) indicate that strain and spatial strain gradients magnitudes are insensitive to changes in ν and therefore for all models in this study we specified ν = 0.2 since that is the most common value used in the literature.(Villard et al. 2005)

To simulate lung deformation during the respiratory cycle, several boundary conditions and loads were applied to the model (Fig. 1D). Because the patient specific geometries only represent the section of the lung scanned during the MRE protocol, a rolling/sliding boundary condition was applied to the upper and lower surfaces to restrict apical/basal movement to be in-plane only. To account for the restriction of lung deformation due to organs/tissues in the mediastinum cavity, we applied a spring foundation boundary condition to the mediastinal surfaces of the lung. Specifically, a restoring force that is linearly related to the local deformation field was applied and the spring constant in the normal and shear/tangential directions was specified as

$$\begin{array}{cc}{k}_{n}=\frac{{E}_{s}\left(1-{v}_{t}\right)}{{d}_{s}\left(1+{v}_{t}\right)\left(1-2{v}_{t}\right)}& {k}_{s}=\frac{{E}_{s}}{{d}_{s}\left(2\right(1+{v}_{t}\left)\right)}\end{array}$$

5

Where tissue stiffness was set to Es=800kPa, Poisson’s ratio to \({v}_{t}\)=0.4 and tissue thickness to \({d}_{s}\)=10 cm. These values allowed for medial movement of the lungs without affecting the strain magnitudes within the inner region of interest. In addition, a sensitivity analysis was conducted by varying Es to investigate the effect of this parameter on simulation results, and it was determined that varying Es in the range shown in Table 1 did not affect the computed strains or strain gradients (Supplemental Fig. 1B). To simulate negative pressure ventilation, a sinusoidal vacuum pressure was applied to the outer surfaces of the lung section, and in this study, a vacuum pressure that ranged from 0 to a maximum vacuum pressure of -p was used with a sinusoidal frequency of 0.2 Hz. The maximum vacuum pressure p was determined for each patient specific model, which was based on the resultant volume ratio. Specifically, for an average normal adult, inspiration leads to an increased lung volume of ~ 0.5L (i.e. the tidal volume) about the functional residual capacity of ~ 3L. As a result, normal tidal ventilation resulted in a volume ratio change of 3.5L/3L = 1.167. As a result, p was determined in each normal adult patient specific model to give a volume ratio change of 1.167. Although the average maximum vacuum pressure in normal adults was 751.8 ± 47.3 Pa, pwCF had a 50% decreased tidal volume as compared to healthy subjects(Colasanti et al. 2004). Therefore, p was set in these models to achieve a volume ratio change of 1.083 and this resulted in an average maximum vacuum pressure in pwCF of 467.7 ± 29.0 Pa.

Simulation of lung deformation (Fig. 1F) was then used to calculate both local strain and spatial strain gradients at end inspiration in each lung. First, since we are primarily interested in tensile strain, the strain magnitude was quantified using the 1st principle strain \({\epsilon }_{p}^{1}\). Second, the derivatives of \({\epsilon }_{p}^{1}\)were used to calculate the strain gradient magnitude as

$$\nabla {\epsilon }_{p}^{1}=\sqrt{{\left(\frac{d{\epsilon }_{p}^{1}}{dx}\right)}^{2}+{\left(\frac{d{\epsilon }_{p}^{1}}{dy}\right)}^{2}+{\left(\frac{d{\epsilon }_{p}^{1}}{dz}\right)}^{2}}$$

6

**Statistical Analysis**

Distributions of shear stiffness, 1st principal strain magnitude and spatial strain gradient were not found to consistently follow normal or log-normal distributions by Kolmogorov–Smirnov test. Therefore, for each subject, the median values of shear stiffness, strain magnitude and spatial strain gradient as well as the inter-quartile range (IQR) were calculated and paired t-tests or 1-way analysis of variance (ANOVA) was used to document statistical differences between median and/or IQR values. For this analysis, each lung was considered a “subject” and although both lungs were scanned in n = 3 pwCF, only the right lung was scanned in n = 6 healthy adults. Therefore, there were n = 6 “subjects”/lungs for both the healthy and CF groups. Statistical significance was determined at a p = 0.05 level.