Mechanical properties of the airway tree: heterogeneous and anisotropic pseudoelastic and viscoelastic tissue responses.

Airway obstruction and pulmonary mechanics remain understudied despite lung disease being the third cause of death in the United States. Lack of relevant data has led computational pulmonary models to infer mechanical properties from available material data for the trachea. Additionally, the time-dependent, viscoelastic behaviors of airways have been largely overlooked, despite their potential physiological relevance and utility as metrics of tissue remodeling and disease progression. Here, we address the clear need for airway-specific material characterization to inform biophysical studies of the bronchial tree. Specimens from three airway levels (trachea, large bronchi, and small bronchi) and two orientations (axial and circumferential) were prepared from five fresh pig lungs. Uniaxial tensile tests revealed substantial heterogeneity and anisotropy. Overall, the linear pseudoelastic modulus was significantly higher axially than circumferentially (30.5 ± 3.1 vs. 8.4 ± 1.1 kPa) and significantly higher among circumferential samples for small bronchi than for the trachea and large bronchi (12.5 ± 1.9 vs. 6.0 ± 0.6 and 6.6 ± 0.9 kPa). Circumferential samples exhibited greater percent stress relaxation over 300 s than their axial counterparts (38.0 ± 1.4 vs. 23.1 ± 1.5%). Axial and circumferential trachea samples displayed greater percent stress relaxation (26.4 ± 1.6 and 42.5 ± 1.7%) than corresponding large and small bronchi. This ex vivo pseudoelastic and viscoelastic characterization reveals novel anisotropic and heterogeneous behaviors and equips us to construct airway-specific constitutive relations. Our results establish necessary fundamentals for airway mechanics, laying the groundwork for future studies to extend to clinical questions surrounding lung injury, and further directly enables computational tools for lung disease obstruction predictions. NEW & NOTEWORTHY Understanding the mechanics of the lung is necessary for investigating disease progression. Trachea mechanics comprises the vast majority of ex vivo airway tissue characterization despite distal airways being the site of disease manifestation and occlusion. Furthermore, viscoelastic studies are scarce, whereas time-dependent behaviors could be potential physiological metrics of tissue remodeling. In this study, the critical need for airway-specific material properties is addressed, reporting bronchial tree anisotropic and heterogeneous material properties.

NEW & NOTEWORTHY Understanding the mechanics of the lung is necessary for investigating disease progression. Trachea mechanics comprises the vast majority of ex vivo airway tissue characterization despite distal airways being the site of disease manifestation and occlusion. Furthermore, viscoelastic studies are scarce, whereas timedependent behaviors could be potential physiological metrics of tissue remodeling. In this study, the critical need for airway-specific material properties is addressed, reporting bronchial tree anisotropic and heterogeneous material properties. airway characterization; anisotropy; heterogeneity; lung; mechanical properties; viscoelasticity MOTIVATION The billions spent annually and the rank of lung disease as the third leading cause of death in the United States, with higher fatal incidence elsewhere, motivate the necessity for pulmonary mechanics research (1, 7, 10a, 67). Trial and error diagnoses are often a long process facing pulmonary patients (19,22). The rapid evolution of pulmonary research has influenced recent shifts in medical practice, such as a sharp increase in lives saved by reducing the volume of ventilator-pumped oxygen (3). However, the mechanisms relevant to this groundbreaking discovery or to the fatalities linked to lung expansion stresses are still not well understood.
Many chronic pulmonary diseases require an understanding of the fundamental science governing airway mechanics to inform diagnoses, design interventions, and formulate preventative measures. Respiration is an inherently mechanical process, and respiratory complications, such as airway collapse in asthma and bronchitis, involve altered mechanical behavior. However, the field lacks critical measurements of healthy tissue structure and function, let alone detailed information about remodeled and evolved states caused by disease progression and endurance (2,5). To date, the majority of tissue characterization studies reporting the mechanical properties of the airway walls have focused on the trachea, the largest, most accessible passageway (12,45,60,63). Force generation of the active smooth muscle has been examined at different airway levels, but proximal to distal airway stiffness parameters and mechanical tissue responses were not reported (23). Other studies investigated the mechanics of the lung as an entire functioning organ, and only a few explicitly considered the airways beyond the trachea, which support pulmonary function and are the sites of obstruction upon collapse due to remodeling and disease (25,51,59). Although research regarding tracheomalacia and tracheal stent replacements benefit from these explorations, the vast majority of lung diseases are intraparenchymal, impacting the smaller bronchioles observed to occlude during expiation, for which no mechanical property data are currently available (32,36,69). Furthermore, the time-dependent mechanical responses of the airways have been minimally investigated, whereas they have been thoroughly documented and shown to be physiologically relevant in other organs (18,56). As respiration can involve both voluntary and involuntary (due to obstruction or ventilation) variations in breathing volume, frequency, and duration, examination of time-dependent viscoelastic behaviors can potentially introduce a metric for tissue remodeling and disease categorization.
The predictive capabilities of pulmonary mechanics computational models are limited by the unavailability of airwayspecific constitutive relations. Previous studies highlighted the need for relevant experimental data to inform physiological models, especially for the bronchial tree (17,45,48,57). Past studies extensively investigated the buckling and folding behavior of the airways due to inflammation and bronchoconstriction in diseases such as asthma and bronchitis, since the effect is directly tied to airway resistance (4,16,53). Complex models have evolved from simplified analytical models to computational optimization to personalized magnetic resonance images informing patient-specific geometries (15,29,35,38,40,65); yet all of these studies highlighted limitations due to the lack of regional airway material properties and resorted to utilizing generic constitutive laws. The material properties of the trachea cannot necessarily be extrapolated to model the mechanics of distal bronchi, and the current lack of relevant mechanical property data hinders modeling of clinically relevant pathologies and development of predictive bronchial tree occlusion diagnostic tools.
The goal of the present study was to address this gap in knowledge of functionally relevant biomechanical property data through experimental characterization of the heterogeneous, anisotropic material properties of porcine extraparenchymal and intraparenchymal airways (notably airways characterized as cartilaginous bronchi, but not smaller bronchioles terminated by alveoli) (13). This was done through uniaxial tensile testing of specimens from three airway levels (trachea, large bronchi, and small bronchi) and two orientations (axial and circumferential) (see Fig. 1). The data revealed significant anisotropy and regional differences in measured pseudoelastic and viscoelastic tissue behaviors of proximal and distal bronchi airways. These bronchial mechanics findings lay the foundation for future studies to build upon. Furthermore, utilization of airwayspecific material properties reflecting regional and directional variations in properties will enable computational tools to better analyze pulmonary occlusion in the case of disease and excessive distention in pulmonary complications such as ventilated induced lung injury (VILI).

METHODS
Reproducible experimental protocols were developed to allow proper comparative analysis of proximal and distal axial and circumferential tissue samples. Various sample preparation strategies, mechanical testing protocols, and data analysis procedures were explored in preliminary testing before a final procedure was selected for consistent and robust data collection. The justification for the final methodology is presented here.
Specimen preparation and protocol development. Porcine airways are composed of a stiff outer cartilaginous layer with a soft inner tissue layer connecting the freely deforming cartilaginous rings and scales (6). The intact, full thickness of the airway wall was tested to determine the effective mechanical properties without peeling the soft tissue from the cartilage layers (12,45), as preliminary attempts to remove the soft tissue compromised connective fibers. Effectively, the freely deforming, intact soft tissue continuously spanning the cartilaginous surface and connecting the two cartilage rings was tested, as shown in Fig. 2. This had the further advantage of allowing samples to be clamped at the stiffer cartilage, avoiding artifacts due to excessive clamping or damage of the soft tissue. Anatomic airway landmarks (bifurcation of the trachea into the right and left main bronchus) were selected due to their shared formation with human lungs. Therefore, the tracheal bronchus (a porcine-specific additional upper-third main bronchus) was not examined. To balance the statistical model, the same number of samples was excised from the right and left lobe of every pig despite the opportunity to dissect more samples from the right lobe due to the larger volume and location of the heart on the left side of the chest cavity.
Samples were cut by hand with a scalpel and were~3-4 mm wide and 8 -9 mm long, as initial attempts to extract samples with a punch resulted in slippage and inconsistent sample shapes. Excising by hand also allowed navigating and avoiding the complex branching regions of airways where the tissue could not be reliably described as either axial or circumferential, and cartilaginous rings were reinforced and non-uniform. Sample thickness was dependent on the airway wall dimension and decreased in distal airways (Table 2). Proximal to distal airways were examined by excising samples from the trachea to the main left and right bronchi. Specimens were prepared from the trachea and from large and small bronchi (distinguished based on average diameters) and two orthogonal direction-defined specimen orientation (axial and circumferential) within each region. Six samples were taken per specimen category (region and direction) from each of 5 donors. ID, inner diameter. For mechanical testing, the porcine airway tissue was clamped, the crosssectional geometry was measured, and the grip-to-grip distance was used as the initial tissue length (elongated sample grasped at cartilaginous section shown). Samples were loaded at 1%/s deformation rate, preconditioned 5 times to a deformation of 35%, and the response of the 6th loading ramp was used to derive the pseudoelastic response. Samples were then held at peak elongation for 300 s to investigate the viscoelastic stress relaxation response and then unloaded. Image of the opened cylindrical airway illustrates tissue topology and the monopodial nature of porcine airways; branching sections of the tissue were not used as samples.
Uniaxial tensile tests were performed using an Instron 5848 Microtester and 10-N load cell. Samples were mounted by clamping the specimen first to the upper grip and zeroing the load cell, followed by clamping the lower grip while avoiding load application. The gripto-grip distance (5-6 mm) was set as the sample's initial length, and the applied displacement profile was adjusted as a consistent percentage of this initial length. To check slipping and strain distribution, preliminary digital image correlation results of 20 samples were analyzed with the Ncorr package (8) running within MATLAB software (v. 1.2, Mathworks, Natick, MA) found a fairly homogeneous strain distribution across the specimens, with some boundary effects near the grips.
Samples were preconditioned through loading and unloading cycles using a triangular waveform to a peak 35% deformation (protocol illustrated in Fig. 2). Preliminary tests indicated substantial preconditioning effects over the first few cycles but decreased changes between cycles 5 and 10 (Fig. 3A). To maintain practical testing times, experimental samples were preconditioned for five cycles and the sixth loading curve was used to capture the preconditioned pseudoelastic behavior representative of a repeatable, cyclic loading response, a common approach to characterizing the loading response of biological tissues as effectively elastic despite the known presence of inelastic (e.g., viscoelastic) phenomena (18). Whereas previous reports based on the steady-state stress following 2-3 min of stress relaxation (41,45) provide information relevant to sustained deformation, the quasi-static loading response examined in the present study is more reflective of airway behavior during the inspiration phase of periodic respiration. A deformation rate of 1%/s was chosen as minimal change was observed in preliminary tests between a slower (0.2%/s) and faster (5%/s) rates [in agreement with previous observations (60)], as shown in Fig. 3B. To capture viscoelastic stress relaxation, the deformation was held at 35% for 300 s following the sixth loading step, as in previous studies (41,45,56).
Defining a consistent unloaded reference state and establishing its relation to the physiological naturally occurring state of the tissue in vivo is one of the main challenges of ex vivo tissue characterization. No universal technique has been adopted in the biomechanics literature to do so (11,55). Some studies have established a reference preceding data collection via a tare load while others have done so through post-testing data analysis. Regardless, the classical J-shaped curve (linear-concave) is obtained by physically or virtually removing the initial portion of the S-shaped loading curve (convex-linearconcave; Fig. 4, exhibited by nontared samples). A tare load preceding data collection (45,60) is the most common approach, with stress and strain determined relative to this tared reference state. However, a small difference in the low-strain modulus can result in widely varying reference states under a consistent tare load. As previous work suggested that the low-strain modulus of airway tissue could vary substantially with location and orientation (58), we chose to avoid applying a tare load and to instead determine the reference state after testing. For each sample, the linear portion of the load deformation curve was determined by finding the first derivative using two-point forward finite difference and identifying the deformation range where the first derivative was within 10% of the minimum slope value. The transition point between the convex and linear portions of the curve was identified by fitting a line to the linear portion of the curve and identifying the lowest displacement for which the stress deviated from the best-fit line by no more than 1% error (relative to peak stress value). The sample length at this transition point (typically Ͻ5% deformation relative to the unloaded starting point) was treated as the reference for strain calculation, and the load at the transition point was subtracted as a "virtual" tare load (Fig. 5). This posttesting virtual tare load technique was consistent, individualized to each sample's load deformation behavior regardless of region and orientation and was reproducible. One notable limitation of this technique is the shifting of the reference length for strain calculation, resulting in a slight variation in calculated strain rate and peak strain among samples. All Fig. 3. Two representative tissue curves demonstrated the force-displacement range, typical preconditioning response, strain rate behavior, and notable hysteresis as indicated by the area between the upper and lower loading curves. A: 10 cycles of preconditioning show greater changes among initial cycles than later cycles. Cycles 1, 5, and 10 are solid lines. B: the 6th loading cycle for the same sample at 0.2, 1, and 5%/s demonstrates minimal rate sensitivity in this range. To maintain testing time feasibility, experimental samples were tested at 1%/s deformation rate through 5 preconditioning cycles followed by a 6th loading curve and a 300-s hold. Fig. 4. Example load-displacement curve (black) illustrating S-curve typical of biological tissues as a convex portion followed by linear and concave regions (indicative of fiber unkinking, straightening, and load bearing). To replicate in vivo settings, application of a tare load crops the convex portion of the curve, revealing the J-curve typically examined for biological tissues. However, a consistent tare load applied to tissues from varying regions, orientations, and thicknesses could introduce artifact due to different stiffnesses in the small deformation, low stiffness regimen. As an alternative, the transition point (arrow) was identified as the start of the linear portion (gray, extended for visualization) of the loading curve. The force at this point was subtracted as a virtual tare load, and strains were calculated relative to the sample length at the transition point. samples had peak calculated strains of at least 0.25, allowing comparisons of behaviors over that range.
Experimental design. Fresh lungs from five pigs older than 6 mo were obtained from a slaughterhouse (Allied Pringles Food Sales, Oakland, CA) on the day of slaughter. Tissue samples were isolated from three regions (trachea, large bronchi, and small bronchi) and at two orthogonal directions (axial and circumferential), resulting in six sample categories (Fig. 1). Furthermore, for each of these categories, a total of six samples were isolated from each donor, with three samples coming from each of the left and right lung lobes for the large and small bronchi. While some samples were lost due to damage during preparation, or testing, or slipping during loading, 150 samples were analyzed (see Table 2 for sample distribution). Airways were excised from the parenchymal lobes and kept hydrated in the intact specimen's natural blood before uniaxial testing; the introduction of foreign chemicals (even preservation solution) was avoided, and dissection from surrounding tissue was performed just before testing so as to limit the extent of structural changes. Testing was done within 36 h post mortem over two consecutive days (66).
The cylindrical inner diameters of the airways were measured with Vernier calipers along the largest cross-sectional length (average of three measurements) before they were cut open to excise axial and circumferential samples. The width and thickness of each sample in the unloaded state were also measured with Vernier calipers and used to calculate the Piola engineering stress from the measured load. For simplicity, uniaxial strain was defined as the fractional change in length relative to the sample-specific reference length determined through the virtual tare procedure (e.g., stretch ratio minus 1).
Data analysis. Several measurements resulted from this experimental protocol. The pseudoelastic linear modulus was derived as the slope of the linear region of the preconditioned stress-strain curve during loading (analogous to inflation) and was interpreted as a measure of material stiffness during the loading phase of cyclic loading-unloading (analogous to respiration).
The tangent moduli of the stress-strain curve, which provide insight into the relative nonlinearity of the stress-strain curve across sample groups, were determined as the slope of the stress-strain curve for discrete strain values of 0.05, 0.15, and 0.25. A line was fitted to the stress-strain data for a 2% strain range centered at the target strain (e.g., 0.04 -0.06 to determine the tangent modulus at 0.05 strain). For comparison, tangent moduli at the same strain values and the pseudoelastic strain modulus were derived from previously published airway data ( Fig. 6 and Table 1) (12,45).
To compare the viscoelastic stress-relaxation responses, a two-term Prony series exponential decay (g o ϩ g 1 exp ͑Ϫt ⁄ 1 ͒ ϩ g 2 exp ͑Ϫt ⁄ 2 ͒ ) was fitted to the 300-s stress relaxation hold, where t is time, g symbolizes relaxation coefficients, and 1 and 2 are the characteristic time relaxation constants) (18). The percent relaxation was determined as the ratio of the measured stress reduction over the 300 s hold to the initial peak stress. Fig. 5. Representative stress-strain responses for trachea, large bronchi, and small bronchi samples from axial and circumferential orientations. The pseudoelastic linear modulus, defined as the slope of the linear portion of the stress-strain curve (inset graph), demonstrates both anisotropic and heterogeneous responses of airway tissue. Circumferential samples offered less resistance to increasing strain than axial samples, and axial samples transitioned out of the "toe region" to a higher modulus regimen within the tested strain range. Fig. 6. The pseudoelastic linear modulus for trachea, large bronchi, and small bronchi samples from axial and circumferential orientations. Anisotropy was evident as greater pseudoelastic linear moduli for axial samples than for corresponding circumferential samples. Heterogeneity was found within circumferential samples, as the modulus of circumferential samples from small bronchi were significantly larger than that of trachea or large bronchi counterparts. No significant differences in pseudoelastic linear modulus were found among axially oriented samples. Although sample sources and preparations were different in previous studies (12,45) (including determining modulus after stress relaxation), modulus values and reported anisotropy were similar to those in this study. Pseudoelastic linear modulus values are listed in Table 2; *, **, and *** denote P values Ͻ 0.05, 0.01, and 0.001, respectively. All measured or derived results were subject to an optimal Box-Cox transformation and analyzed via a two-factor (airway level, orientation) general linear model with significance at P Ͻ 0.05, with Bonferroni's test for pairwise comparisons within significant effects (Minitab v. 17, Minitab, State College, PA). Tangent moduli were analyzed via a three-factor general linear model (airway level, orientation, strain value). Initial statistical models included lobe (left vs. right) and day of testing as fixed factors, but these factors did not have significant effects and were removed. Although donor animal was initially included as a random variable to create a repeated-measures analysis, this factor was rarely significant and was removed from all models except for the tangent modulus analysis for parsimony. Results are presented as means Ϯ SE; *, **, and *** denote P values Ͻ 0.05, 0.01, and 0.001, respectively. Fig. 3) demonstrate the varying force displacement range, typical preconditioning response, insensitive strain rate behavior, and notable hysteresis indicated by the area between the upper and lower loading curves. Figure 3A demonstrates the impact of preconditioning on the tensile loading and unloading behavior; the pseudoelastic linear modulus decreased from the first to the sixth loading cycle by 2-2.5 times. Circumferential tissue samples exhibited a greater modulus reduction due to preconditioning than axial samples (P ϭ 0.025). For deformations up to 35%, the maximal tissue force response ranged from Ͻ1 N to Ͻ15 N and depended on the region and orientation of the specimen.

Two example tissue curves (shown in
Pseudoelastic response: linear modulus and tangent moduli. Figure 5 displays representative stress-strain responses for samples with closest to the average pseudoelastic linear modulus from each group. Axial samples displayed higher stresses than circumferential samples for the same strain. Axial tissues furthermore transitioned to the expected nonlinear stiffening behavior at smaller strains than their circumferential counterparts. Representative stress-strain curves for axial trachea and large bronchi samples were similar; stress-strain curves for circumferential trachea and large bronchi were also alike. However, the small bronchi axial and circumferential stressstrain behavior differed from other groups. Inset reveals drastic differences in the linear slopes used to define the pseudoelastic linear modulus. The pseudoelastic linear modulus (Fig. 6) varied with both airway size and orientation. Overall, the modulus was significantly greater for axial samples (30.5 Ϯ 3.1 kPa) than for circumferential samples (8.4 Ϯ 1.1 kPa, P Ͻ 0.001). Axial samples exhibited no significant variation among regions in the pseudoelastic linear modulus. In contrast, the pseudoelastic linear modulus of circumferential samples from small bronchi (12.5 Ϯ 1.9 kPa) was significantly greater than those from trachea (6.0 Ϯ 0.6 kPa, P ϭ 0.007) or large bronchi (6.6 Ϯ 0.9 kPa, P ϭ 0.021). For comparison, Fig. 6 also presents values of pseudoelastic linear moduli for trachea axial, trachea circumferential, and bronchi circumferential samples derived from previously reported data (12,45). Sample preparations and species differed: previous studies had excised soft tissue from the inner layer of the airway wall and measured the stress values after tissue relaxation. Despite these differences, our findings are consistent with the magnitude and limited information on anisotropy found in previous reports.
Average tangent modulus values and SEs at strain values of 0.05, 0.15, and 0.25 for the six tissue groups are reported in Table 1. Consistent with the expected nonlinear response of the tissue, tangent values for the 0.25 strain were significantly increased comparied with the 0.05 and 0.15 strains for axial trachea (P Ͻ 0.001), axial large bronchi (P Ͻ 0.001), axial small bronchi (P ϭ 0.003), circumferential large bronchi (P Ͻ 0.001), and circumferential small bronchi (for 0.05 strain P Ͻ 0.001, for 0.15 strain P ϭ 0.003); no significant differences in tangents for increasing strain value were found for circumferential trachea tissue. Consistent with the pseudoelastic linear modulus, mean axial tangents for trachea and large bronchi samples at each strain level were multiple times greater than their circumferential counterparts (P Ͻ 0.001), indicative of pronounced anisotropy. Small bronchi exhibited greater axial tangents than circumferential tangents at 0.05 and 0.15 strains (P Ͻ 0.001 and P ϭ 0.004), with no significant difference at the 0.25 strain as circumferential and axial tangent values merged closer. Circumferential tangents for the 0.25 strain were greater for large and small bronchi compared with the trachea (P ϭ 0.001 and P Ͻ 0.001), whereas axial tangents for the 0.25 strain were found to be greater for the trachea than for small bronchi (P Ͻ 0.001). The greatest tangent moduli was exhibited by axial trachea samples at the 0.25 strain (see Fig.  5 and Table 1). Among circumferential samples, the small bronchi had the highest average tangent modulus for each strain value. The observed nonlinear response and anisotropic behavior were similar to those of previous studies (12,45) when subjected to the same analysis.
Viscoelastic response: stress relaxation and time constants. Figure 7 shows the relaxation response over 300 s normalized to peak stress (means Ϯ SE) for the six groups. No groups exhibited greater than 50% percent stress relaxation over the testing period. Axial samples appeared to approach a plateau within the 300 s, whereas circumferential samples did not, qualitatively indicating anisotropy in the viscoelastic behavior of airway tissue. Consistent with this qualitative observation, quantitative metrics of viscoelastic relaxation (percent relax- Fig. 7. Stress-relaxation response over a 300-s hold. Average stress-time response Ϯ SE shown for axial and circumferential tissue specimens from the trachea, large bronchi, and small bronchi. Curves exhibited relaxation typical of biological soft tissues (18,26), with relaxation no more than 0.5 of peak maximum stress during the 300-s hold. From these relaxation curves, the fractional relaxation and time constants 1 and 2 in a 2-term Prony series relaxation function were calculated. Viscoelastic parameters are shown in Fig.  8, and values are listed in Table 2. ation and time-constants) exhibited both heterogeneity and anisotropy, as shown in Fig. 8 and listed in Table 2. Figure 8A demonstrates that percent relaxation was significantly anisotropic; axial samples exhibited significantly smaller percent relaxation (23.1 Ϯ 1.5%) than their circumferential counterparts (38.0 Ϯ 1.4%) (P Ͻ 0.001). Trachea axial and circumferential samples had significantly higher percent relaxations (34.5 Ϯ 1.7%) than the large (29.9 Ϯ 1.4%, P ϭ 0.005) or small bronchi (27.3 Ϯ 1.2%, P Ͻ 0.001).
The second (longer) time constant 2 (Fig. 8C) was nearly an order of magnitude larger than 1 , and had significant anisotropy for only the trachea samples (154.1 Ϯ 9.0 s axial compared with 99.1 Ϯ 5.8 s circumferential (P Ͻ 0.001). Heterogeneity was noted in the circumferential direction only, between the trachea and distal regions: 2 for the large bronchi was 136.0 Ϯ 6.7 s (P ϭ 0.002) and 137.7 Ϯ 5.3 s for the small bronchi (P Ͻ 0.001).

DISCUSSION
This research provides an intra-study comparison of the pseudoelastic and viscoelastic properties of proximal and distal bronchial tissues. Prior research has demonstrated that trachea tissue from humans and animal species is anisotropic, also confirmed by this study; however, exploration of the the mechanical responses of tissue from the distal airways embedded within the parenchyma underscoring heterogeneity and anisotropy, as additionally reflected in the viscoelastic response, has been limited. We have considered cartilaginous bronchi in this study, providing material property values needed in pulmonary airway mechanics research. The results of this work motivate consideration of tissue evolution in distal terminal airways, an important consideration given that they contribute nearly the entirety of total lung volume (62) and are the primary region of airway occlusion due to lack of cartilaginous reinforcement.
Pseudoelastic properties. Airway tissue exhibits the classical preconditioning, hysteresis, nonlinear elastic, and viscoelastic effects behavior of many well-established biological structures (10,18,26). Similar to other biological vessels, such as arteries and intestines (24,61), circumferential airway tissue samples were found to be more compliant than their axial counterparts. These findings are consistent with our visual observations of airway tissue striations (see Fig. 2) and with past studies citing elastic fiber tissue oriented about the axial direction (31,64) causing increased pseudoelastic linear modulus in the axial direction compared with the circumferential tissue orientation. The greater circumferential compliance, combined with a greater circumferential wall stress produced by internal air pressure, facilitates an increase in airway volume during inspiration that will be limited by the nonlinear stiffening at higher strains. This would be expected to reduce airway resistance but potentially could increase airway dead space, and future modeling studies incorporating these measured mechanical properties may provide relevant functional insights. Previous studies reported a relaxed modulus (expected to be lower than the pseudoelastic modulus reported here) from bovine or ovine airways with the innermost soft Fig. 8. Viscoelastic parameters describing relaxation of axial and circumferential tissue samples from the trachea, large bronchi, and small bronchi depicted anisotropy and heterogeneity. A: fractional relaxation defined as the ratio of stress reduction (peak minus final) value to peak stress exhibited significant anisotropy, Axial samples relaxed less and had smaller fractional relaxation values than their circumferential counterparts. Across both axial and circumferential tissue samples; samples from trachea relaxed more than those from large or small bronchi. B: the 1st time constant 1 exhibited significant anisotropy for only trachea samples and was significantly different for trachea samples than for those from large or small bronchi. C: the 2d time constant 2 also exhibited significant anisotropy in only the trachea samples and was significantly different for trachea samples than for those from large or small bronchi for only circumferential samples. Viscoelastic parameter values are listed in Table 2; *, **, and *** denote P values Ͻ 0.05, 0.01, and 0.001, respectively. connective tissue excised from the airway wall, precluding direct numerical comparison with moduli from this study; nevertheless, the notable anisotropy and increasing tangent moduli with increasing strain indicates a general consistency with prior findings (12,45).
Tissue from the three regions (selected from the extraparenchymal trachea to the furthest-most intraparenchymal region where specimen dimensions were still feasible for testing) and two orthogonal directions exhibited distinct pseudoelastic behaviors not previously reported. Noble et al. (45) tested circumferential deformation between the trachea and the immediately adjacent large bronchi and found no difference in the circumferential stress-strain response, in agreement with this study. We found that the pseudoelastic linear modulus of axial tissue specimens did not significantly differ across regions, whereas circumferential tissue specimens were heterogeneous, specifically for the small bronchi. This is similar to the phenomena observed in aortas and intestines, where distal locations differ more drastically in their circumferential elastic behavior but do not exhibit as much variation in axial behaviors (24,61).
In contrast, the circumferential pseudoelastic linear modulus of the small bronchi was notably increased more than twofold compared with proximal airways. It should be noted that the aspect ratio of test specimens varied across regions, as samples were prepared with similar length and width but the thickness of the airway wall varied both among and within regions. Across all samples, no correlations were found between modulus and the measured thickness (nor width, length, and crosssectional area), supporting the view that differences in material properties are related to intrinsic compositional and structural variations in bronchial tree tissue. In the future, use of histological and biochemical analyses will allow exploration of the microscopic and layered changes in the tissue properties between the trachea, large bronchi, and small bronchi to better understand collagen and elastin roles in the tensile variation.
The airway tissue material property analysis is fairly consistent with strain measurements conducted using rat lung imaging, where axial and circumferential deformations of the small and large airways were measured for large tidal volumes and various levels of positive end-expiratory pressure (58). Axial strains in the rat airways were found to be smaller than circumferential strains, yet no significant difference was observed in the axial strain measurements between the large and small airways. In contrast, heterogeneity was found for circumferential strains, with greater circumferential strain observed in the small airways of the rat. Our results are not fully consistent with this observation nor with the long-standing assumption that distal airways are more compliant (34,58). Past studies have found no significant differences in inflation and deflation specific compliances of porcine stem bronchus (41). Whereas the absolute porcine airway wall thickness decreased in the distal bronchial tree, the ratio of thickness to cross-sectional radius increased, as in past findings (30), suggesting that for a comparable pressure the average wall stress is lower in the small bronchi than in the larger airways. This combination of lower stress and higher circumferential pseudoelastic modulus leads to the expectation that the relative expansion of the small bronchi would be less than that of the larger airways. Additionally, it should be noted that the in vivo response (as in Ref. (58) is not purely a function of airway mechanical properties but is also influenced by other factors involved in the overall system of respiration such as surfactant-filled mucus on top of the airways, which varies the surface tension and results in reduced work of breathing (50). Future studies are needed to resolve this apparent discrepancy and more closely relate the local material properties to overall lung function.
The nonlinear transition point of representative axial samples in Fig. 5 demarks a fairly abrupt increase in tensile modulus. While soft tissues in tension typically display strainstiffening manifested as increased tangent moduli, in the context of proximal airways, this transition has been described as the "strain limit" and was surmised to be a protective mechanism (58,60). While the transition strain varies considerably among regions, it is interesting to note that (due to differences in the linear modulus) this transition occurred at similar stress levels across regions, suggesting a possible adaptation to developed airway stress. As axial and circumferential samples were subjected to shared testing protocols, stress-strain curves of the more compliant circumferential samples did not generally extend fully into the higher stiffness region. Nevertheless, qualitative examination of stress-strain curves suggested that the circumferential samples began this transition at stress levels roughly one-half of those seen in axial samples. Pressurevessel theory predicts that the average wall stress produced by internal pressurization will be twice as large in the circumferential direction as in the axial direction (20), suggesting that, for a given airway size, the strain limit will be reached first in the circumferential direction. While the absolute airway wall thickness decreased in the distal bronchial tree, the ratio of thickness to cross-sectional radius increased as in past findings (30), suggesting that for a comparable pressure the average wall stress is lower in the small bronchi than in larger airways. If the circumferential strain limit is also comparable across airway levels, this suggests that large airways would reach the strain limit first, strongly limiting further radial expansion, whereas small airways would continue to expand with increasing pressure. Although additional experiments across a wider strain range would be required to support this conjecture, such a phenomenon could help to explain the mechanism of small airway damage due to excessive forced expansion in ventilatorinduced lung injury and similar traumatic lung diseases (3). Viscoelastic properties. Minimal attention has been given to the region-and orientation-dependent viscoelastic response of the bronchial tree. Findings in this study are similar to the responses of many biological tissues: the tissue is sensitive to preconditioning, though not as strongly as in some other organs (10). Similarly, the tissue exhibits stress relaxation, with reduction by up to one-half of the peak stress over the hold period. Axial tissue samples reach the relaxed stress value faster than circumferential tissue samples, which did not clearly approach steady-state values during the testing period. This anisotropy is also evident in the percent relaxation of the tissue. Trachea samples exhibited greater percent stress relaxation than the large or small bronchi for both axial and circumferential tissues, possibly because of its denser smooth muscle layer. Although the stress relaxation test is not highly physiological, it is convenient for material characterization and is influenced by the same mechanisms producing creep (progressive strain under maintained stress) and hysteresis under cyclic loading. Further testing in those modes might reveal insights more directly relevant to physiological airway response during respiration.
The only point of reference for the viscoelastic measurements performed in this study is a report of the ex vivo stress-relaxation response of various layers of the human trachea (56). Although the species and sample preparation differed, the viscoelastic time constants 1 and 2 and the percent stress relaxation derived in our study are within the same range and order of magnitude as those reported by Safshekan et al (56). It is important to note that this comparison (as well as our comparisons between groups) is appropriate and valid, because the two studies fitted the same relaxation function (a two-term Prony series) to data describing an identical 300-s relaxation. Unlike our study, Safshekan et al. examined stress relaxation at different strains, leading to the conclusion that trachea tissue exhibits fully nonlinear viscoelasticity, as time constants depended on strain level. Since we examined relaxation at only a single strain, we are unable to make inferences regarding nonlinearity and for convenience assumed quasi-linear viscoelasticity to facilitate comparison of percent relaxation. While additional studies are required to explore thoroughly the regional viscoelastic responses of airway tissues, a primary motivation for the present study was an examination of the elastic response for the purpose of computational modeling. In that context, the limited stress relaxation (especially for the axial tissues) over short periods suggests a predominantly elastic behavior [as inferred for the trachea (56)] over time scales relevant to respiration, indicating that a nonlinear elastic material model is a suitable approximation to the actual tissue response. Should a viscoelastic model be employed, however, results of this study provide some insights to inform modeling choices.
The time constants did not display the sample categorical anisotropy as fractional relaxation values did. Only axial trachea tissue had higher time constants than circumferential. Large and small airways did not display significant anisotropy. Additionally, heterogeneity was observed between trachea and other regions, but no significance was found between large and small airways. Potentially, assuming uniform response of the bronchial tree's intraparenchymal regions may be sufficient for modeling purposes. The same claim cannot be made for the pseudoelastic and fractional relaxation response: anisotropy and heterogeneity should be accounted for.
Limitations and future directions. Although the use of porcine tissue allows characterization of healthy tissue behavior, some caution must be exerted in extrapolating material values to describe the behavior of human airway tissue. Nevertheless, as this study considered nonbranching regions, the monopodial and bipodial airway structure of porcine vs. human airways is irrelevant, and the similarity between human and animal lung structure allows that qualitative inferences may be made (42,47). The onset of disease will likely cause variation in measured material properties: prior works have measured increased airway resistance and rigidity due to airway remodeling in chronic diseases (5,9,17,39). Additionally, although this study provides pseudoelastic and viscoelastic material values, and facilitates new insights into regional and directional differences in airway tissue behavior, the in vivo loading environment is substantially more complex, and the tissue response may exhibit coupling between loading directions. Complementary testing methods such as tube inflation or planar biaxial testing can be utilized in the future to examine interdependence between loading directions. Activation of the smooth muscle surrounding the airway wall likely will lead to differences in observed anisotropic and heterogeneous stiffness responses compared with those seen in passive uniaxial tensile tests (31). Although such methods would extend this investigation by more closely reproducing physiological pulmonary loading states, this study provides new insights into airway tissue mechanics and allows interpretation of airway mechanical properties in the context of a large body of literature involving uniaxial testing of soft tissues.
The uniaxial tensile tests performed in this study reproducibly measured the pseudoelastic and viscoelastic mechanical properties in the axial and circumferential directions of the trachea, large bronchi, and small bronchi; however, measurements ex vivo eliminate tissue prestretch and boundary conditions. The extraparenchymal trachea is enclosed adjacent to the neck muscles and the spine, while the transitioning large bronchi are behind the heart, and small bronchi are intraparenchymal. Isolating tissue from supporting structures may be the cause of discrepancy with in vivo imaging studies (28). Residual stresses are unlikely to be a source of discrepancy, as human and porcine airways have not been found to display opening angles (in contrast to sheep and rabbit airways) (42).
Additionally, the parenchyma is known for its role in tethering open the airways (36,49); nevertheless, it is peeled off in ex vivo airway measurements, which can potentially impact the in situ prestretch the airways may be subjected to. It has been noted that static length measurements on ex vivo axial tissues do not suggest recoil compared with in situ measurements (44), consistent with the general agreement between this study and conclusions of in situ axial strain measurements in rat lung imaging (58). Moreover, no such observations of existing prestretch for the circumferential orientation have been explored, thus restricting the explanation for why small airways strain imaging results disagree with our material property measurements. Stiffer ex vivo small airways are not implausible, as detachment and removal of parenchymal teth-ering may decrease overall compliance with the loss of the soft parenchymal substrate contribution (33). A recent in situ preliminary study measured the elastic response of rabbit trachea by utilizing optical coherence tomography (OCT), and another used image registration to estimate heterogeneous Young's modulus of lung lobes (54,57); but until OCT can be extended to region-and direction-dependent exploration, and to avoid assumptions made by indirect and inferred measurements of stiffness from simplified pressure-vessel theory and fluid-structure interaction simulations, ex vivo measurements remain an invasive yet direct measurement (46,54,58).
The observations from this experimental study directly impact predictive computational tools. Stress-strain behaviors will inform future constitutive laws, resulting in more accurate airway-specific simulations instead of generic isotropic biological constitutive models, and informing the degree of anisotropy and material property gradient from the proximal to distal bronchial tree (27). Additionally, the existence of cartilaginous rings from the proximal to distal airways has been verified (37), and upon examining the collapse morphology in these porcine tissues, this study's finding conclude that frequent circumferential airway collapse requires inclusion of rigid cartilage rings in the previous bilayered models of airways (14,35,65,68). These rigid supports will drastically decrease the bifurcation pattern, and the resulting folds will obstruct the airways depending on cartilage ring failure instead of relying solely on the soft connective tissue. These experiments motivate revision of airway mechanics modeling.
The lungs can be controlled voluntarily, unlike many other organs. Therefore, time-dependent behaviors, such as breath duration, frequency, and speed, suggest a potential metric for disease categorization and progression through documentation of airway viscoelasticity. This study's anisotropic and heterogeneous findings regarding percent stress relaxation and time constants initiate this exploratory direction. Paired with findings of effective pseudoelastic moduli, the representation of airway mechanics can potentially be simplified as rheological models (i.e., standard linear solid), and these findings can be integrated into future patient-specific, physiologically relevant geometric and material simulations.
This exploration of tissue variation in the bronchial tree facilitates insight into clinical phenomena regarding observations of asthma hyperresponsiveness, airway stent design, and ventilated induced lung injury (3,21,52). In the case of the last, excessive ventilation, while necessary to stop alveoli collapse, can be damaging if small airways mechanotransduct higher stresses to the body. Additionally, airway stents require tissue material properties for a sufficient and integrative design. Last, airway wall heterogeneity has been credited with the hyperresponsiveness observed in asthma. Therefore, by capturing the heterogeneous and anisotropic material behavior of the bronchial tree and airways, we are one step closer to comprehending the underlying governing physiology of the lungs.

CONCLUSION
Advancing predictive computational tools for airway obstruction through tissue deformation and fluid-structure interaction models relies heavily on the characterization of bronchial mechanics. A series of uniaxial tensile tests of various regions and orientations were performed to supply the pulmonary mechanics community with essential airway tissue material property values. Given the inherent structural differences along the bronchial tree, past assumptions of airway uniformity are insufficient for predictive measures or computational models; results conclude both directional and regional variations. Anisotropy is consistent throughout the bronchial tree, and heterogeneity becomes more pronounced in distal regions. The gradient of evolving material behavior motivates explanations that agree with imaging strain deformation measures, and the collection of viscoelasticity measures within the bronchial tree. The measured pseudoelastic linear modulus, stress relaxation, and time constants emphasize noteworthy airway mechanics variations. Our results address a void in lung mechanics research, further characterizing distal bronchi, the site of diseases such as asthma and bronchitis, to answer questions surrounding pulmonary medicine.