Surfactant-Mediated Airway and Acinar Interactions in a Multi-Scale Model of a Healthy Lung

We present a computational multi-scale model of an adult human lung that combines dynamic surfactant physicochemical interactions and parenchymal tethering between ~16 generations of airways and subtended acini. This model simulates the healthy lung by modeling nonlinear stress distributions from airway/alveolar interdependency. In concert with multi-component surfactant transport processes, this serves to stabilize highly compliant interacting structures. This computational model, with ~10 k degrees of freedom, demonstrates physiological processes in the normal lung such as multi-layer surfactant transport and pressure–volume hysteresis behavior. Furthermore, this model predicts non-equilibrium stress distributions due to compliance mismatches between airway and alveolar structures. This computational model provides a baseline for the exploration of multi-scale interactions of pathological conditions that can further our understanding of disease processes and guide the development of protective ventilation strategies for the treatment of acute respiratory distress syndrome (ARDS).


INTRODUCTION
The lung is an extraordinary example of a physiological organ whose stability and function depends critically upon multiscale geometric interactions and processes that interlink tissue and biofluid mechanics. Airflow during inspiration and expiration changes the lung volume, causing the geometric change of highly compliant pulmonary tissues. During this volume change, tissue-level interactions like parenchymal tethering and molecular level interactions like surfactant re-distribution take place simultaneously in order to stabilize the lung. The lung's large surface area to volume ratio suggests the great importance of surfactant-mediated surface tension to lung stability. Meanwhile, the mechanical interdependency of lung units provides parenchymal tethering support for compliant airways that is instrumental to stabilizing the lung by maintaining adequate airflow to foam-like acinar structures that, in turn, support compliant airways. The large liquid-covered surface of the pulmonary tissue requires surfactant transport to dynamically modulate the surface tension during volume cycling in order to prevent fluid-structure instabilities, as it reduces the work of breathing (West, 1977(West, , 2012Notter, 2000). In this paper, we describe the development of a multi-scale model of a healthy lung from which future studies of pathophysiological systems can be explored.
The computational investigations of lung mechanics have been conducted on different scales. For example, at the organ level, our model is based upon an algorithm developed by Tawhai et al. who presented a computational method for generating the conducting airway network (Tawhai et al., 2000). At the tissue level, Lambert et al. (1982) and Fujioka et al. (2013) have developed computational models to describe the non-linear properties of the acinus and conductive airways, respectively. Micro-mechanical models of parenchyma support are based upon the recent work by Ryans et al. (2019), that is based on the foundational work of Wilson (1972), Lai-Fook et al. (1977, 1978, and Anafi and Wilson (2001). Surfactant transport properties are based on multilayer transport processes described by Krueger and Gaver (2000), which models surfactant behavior originally described by Clements (1957), Horn and Davis (1975), Smith and Stamenovic (1986), and Schurch (1995).
While these and many other studies have provided useful information and analysis about pulmonary mechanics on separate scales, the interactions between all of these interlinking processes have not been modeled directly. Ryans et al. demonstrated a reduced-dimension modeling approach for investigating the recruitment/de-recruitment of the airways without full fluid-structure interactions or surfactant transport processes (Ryans et al., 2016). Wall et al. (2010) developed a 3-D lung model with CT-based geometries up to a maximum of approximately seven generations, including not only airway wall deformability but also the influence of surrounding lung tissue. Filoche et al. (2015) presented a mathematical model of surfactant replacement therapy in a 3-D lung structure to investigate whether the instilled surfactant mixture actually reaches the adult alveoli/acinus in therapeutic amounts. Interlinking the multiscale processes is likely to introduce forcing and response characteristics that are not investigated in isolated models of isolated scales. For example, the multilayer surfactant behavior has heretofore been modeled in an oscillating bubble system with a pre-defined bubble radius. In contrast, in the pulmonary system the radius of airways and acini are driven by an oscillating pleural pressure on the surface of the lung. This, in turn, induces a heterogeneous distribution of mechanical stresses within the lung that dynamically modifies the lung sub-structure. In this process, the geometry of the airways and alveoli are modified, which induces surfactant transport that modulates the surface tension of the lining fluid on the internal structures of the airways and alveoli and, once again, modulates the tethering pressure and geometry of the system (Perun and Gaver, 1995). These feedback processes, when fully modeled, may elucidate pulmonary behavior that cannot understood from the reductionist approach and may exhibit emergent behavior (Suki and Bates, 1985).
In this study, our goal is to build a computational model of an adult human lung that combines dynamic multiscale interactions between ∼16 generations of airways and subtended acini generated by a space filling algorithm, including airflow, multilayer surfactant transport and airway/alveolar interdependency. We explore this system under an oscillating pleural pressure that is applied to simulate the breathing of a healthy lung. Our goal is to develop this FIGURE 1 | Schematic of a "simple" lung unit with one terminal airway attached to an acinus. model as a useful tool to understand the multi-scale lung mechanics and a critical baseline for the future study of the pathological conditions.

Framework
We have developed a rational engineering design approach for elucidating the surfactant physicochemical interactions and the mechanical interdependency within a healthy lung. It is assumed that under healthy conditions, the pleural pressure (P PL ) is always negative and varies with time during the normal breathing cycle. The framework of this simulation approach is based upon three components: I. Tissue Components represent the series and parallel segments of airways and acini. The 3-D structure of airways and acini is determined using an anatomically based space-filling algorithm (Tawhai et al., 2000) in order to introduce inter-connections between airways and acini via parenchymal tethering, which is heterogeneously distributed. II. Air Components represent the volume of the air spaces, with temporal changes based upon the flow rate within each segment. Segmental flow rates are determined by regional pressure differences and the status of each of the compliant airways and acini as determined from the Tissue Component. III. Liquid Components represent the lining fluid and surfactant in the airways and acini, which dynamically modulate the surface tension. This, in turn, influences the Tissue Component stresses and thus affects flow rate.

Interactions Among Components
To illustrate the concept of interactions among the three components, consider the presence of a terminal airway surrounded by parenchyma and connected to an acinus at one end as seen in Figure 1, where P PL is the pleural pressure, P AC is internal acinar pressure, R AW is the radius of airway, Q is the airflow toward the acinus. The pressure at the open end of the airway is considered a constant P 0 to simplify this situation.

I. Tissue Components
These are defined by airways and acini, shown in simplified form in Figure 1. The components are defined by a. Airways: The dynamic radius of the airway is defined by the tube law (Lambert et al., 1982) as a function of the airway transmural pressure (P TM,AW ), described conceptually as: where P int,AW represents the internal pressure from the internal air pressure and the Laplace pressure drop from surface tension, and thus is modulated by Air and Liquid Components. P ext,AW represents the external stresses from the time-dependent pleural pressure (P PL ) and is modulated by the parenchymal tethering using a shear modulus (G eff ) representation of the acinar support as described in Ryans et al. (2019). b. Acinus: Similarly, from Fujioka et al. (2013), the volume of acinus is a function of acinar transmural pressure where P int,AC represents the internal pressure from the acinus air pressure and the Laplace pressure drop from surface tension within alveoli that create the acinus, and P ext,AC = P PL .
In addition to structural coupling that exists through P int,AW and tethering forces between the acini and the exterior of airways, internal coupling exists through internal pressures because P int,AC is a function of the acinus volume. This sets P AC as the end-condition of the Air Component system, which is necessary for determining the flow rates in the system, as described below.

II. Air Components
These are driven by airflow to and from airways and the subtending downstream acinus in the following manner: a. Airflow (Q): Determined by the pressure difference at the two ends of the airway (P AW ), and the conductance (C) of the airway, Q = CP AW . b. Air Conductance (C): We simplify conductance as a function of the pressure-dependent airway radius, viscosity and length, following a Poiseuille relationship. c. Conservation of volume dictates that the flow determines the change of volume of the downstream acinus (V AC ).
As can be seen from the above descriptions, the Air Components create an internal coupling to the Tissue components through P int,AW and P int,AC .

III. Liquid Components
These determine the distribution of fluid in the airways and alveoli. Two fundamental aspects are imposed: a. Conservation of mass: The volume of liquid is conserved within each airway and acinus. So, as an airway or acinus changes shape, the liquid film thickness changes accordingly. b. Surfactant physicochemical interactions: Conservation of mass of surfactant is enforced. As the airway and alveolar surface areas change, this triggers multilayer surfactant transport (Krueger and Gaver, 2000). Through a surface tension equation of state, a dynamic surface tension exists and leads to a realistic full-lung pressure-volume loop.
As can be seen from these interactions, the liquid components directly influence P int,AW and P int,AC in Tissue Components.

System Evolution
The key to the modeling approach is the evolution of the system as each component interacts with each other. To evolve the system, the time-dependent pleural pressure P PL (t) is prescribed and the rest of the system is updated with an adaptive time increment t from the current state. The new state after the time increment is solved by the method described in the Modeling Implementation section and the breathing cycle is tracked.

Modeling Overview
To investigate surfactant-surface tension interactions and the mechanical interdependency within a healthy lung, the algorithm shown in Figure 2 was used to simulate the whole lung respiration. Initially the time-dependent pleural pressure P PL (t) is prescribed at the beginning of each time step. This provides the driving pressure for ventilation, as exists for normal respiration. We define the fundamental solution vectors for the airways, acini where P is airway/acinus pressure, R is airway radius, V is volume of acinus, M is a vector consists of the surfactant masses in each region of the multilayer structure, including M p for the primary layer, M s for the secondary layer, M b for the bulk liquid. Each airway or acinus is represented by two solution vectors, one of which is AC i (t) or AW j (t) containing tissue component properties, and the other is M k (t) containing liquid component properties. AC i (t), AW j (t) and M k (t) are solved simultaneously based on the mechanical relationship described in Network Flow Dynamics using a differential-algebraic equation systems solver.

I. Airway Morphology
The domain of one half of the lung is generated using an anatomically based space-filling algorithm from Tawhai et al. (2000). It consists of a network of airways that has up to an equivalent Horsfield 15th generation (bottom-up with terminal airway defined as generation 1) that terminate into acinar regions with mechanical properties described in the "Acinus Component" section (Fujioka et al., 2013). The equivalent Weibel generation (top-down from trachea as generation 0) of terminal airways and the residual volume of their attached acini are shown in Figure 3. We note that structural differences exist among airways and acini, which will appear as important contributors to the heterogeneity in surface tension and parenchymal tethering described in the "Results and Discussion" section. Airways and acini are positioned in a 3-D space in order to build inter-connections between the airways and acini (Tawhai et al., 2000). The morphology of the pulmonary airways was based on Lambert et al. by providing the airway maximum radius (R AW,max ) and normalized cross-sectional area α at P TM,AW = 0, i.e., α 0 = A 0 A max (Lambert et al., 1982).

II. Compliant Airway Components
We used the airway morphological data to construct a tube law that describes the radius of a compliant airway as a function of the airway transmural pressure. For a healthy lung, we assume that the acinus is the most compliant region of the system, so that the parenchyma tethers the embedded airway open during inspiration, as described in the "Parenchymal Tethering" section, below. To do so, we defined a tube-law that follows a sigmoidal form that matches the airway-based morphological parameters. The general form of the tube-law is defined as where P TM,AW is the transmural pressure defined in Equation (1), R AW,max is the maximum radius of this airway. The parameters a and b are generation-based constants shown in Table 1. The tube law behavior through the Weibel 6th generation as well as experimental data for the first 3 generations are shown in Figure 4 (Lambert et al., 1982). The internal airway pressure P int,AW is subject to the airway pressure determined by the average of the upstream and downstream pressure as well as the Laplace pressure drop from the liquid lining of the airway. The Laplace pressure drop and airway pressure can be expressed respectively as where P AW is the airway pressure, γ is the surface tension of the liquid lining, R in,AW is the radius from the center line Frontiers in Physiology | www.frontiersin.org FIGURE 4 | Tube law through the Weibel 6th generation. Airways smaller than the 6th generation are assumed to have the same tube law as the 6th generation. Experimental data is shown for the first 3 generations.
of the airway to the interior surface of the liquid lining, P DN is the downstream pressure, P UP is the upstream pressure. For the airway shape calculation, we assumed the pressure is uniform within an airway and determined the average pressure by Equation (5). The external airway pressure P ext,AW is subject to the pleural pressure and the parenchymal tethering from surrounding acini. We modeled tidal breathing by forcing the pleural pressure to be prescribed by a sinusoidally varying waveform with a 1:2 inspiration-expiration ratio that varies between −5 cmH 2 0 and −8.2 cmH 2 0 with a period of 5 s (Swan et al., 2012). The parenchymal tethering pressure is determined based on the work from Ryans et al. detailed in Parenchymal Tethering (Ryans et al., 2019). The peri airway pressure can be expressed as where P PL is the pleural pressure, G eff is the effective parenchyma shear modulus, R R is the fractional radius change from the "hole" radius of parenchyma as described below in Equation (12). From Equations (4) to (6), Equation (1) can be rewritten as (7) where the transmural pressure varies with time due to the change in pleural pressure and radius change of compliant airways.
To determine the length of airways, it is assumed that the length change of airways is proportional to the size change of their surrounding parenchyma. The parenchyma model is described in Parenchymal Tethering.

III. Acinus Component
The morphological information of the representative acinar regions connected to the terminal airways of the lung is based on the octahedral models depicted in Fujioka et al. (2013). We represent the behavior of a block of alveoli (∼1,200 alveoli representing ∼1.7 × 10 −3 cc of tissue) by a pressurevolume relationship derived from simulation results of this block. We follow the form suggested by Venegas et al. (1985) that successfully models whole-lung mechanics. As such, the pressure-volume relationship of this model is governed by where V AC is the acinar volume, V RV is the acinar volume at residual volume (RV), P TM,AC is the acinar transmural pressure defined in Equation (2), a-d are constants that describe the mechanical characteristics of the parenchyma. Since we add the dynamic surface tension to the Laplace pressure drop in our form of P TM,AC in Equation (11), this structural component is based upon parameters describing the tissue behavior only (zero surface tension model) with a = 1.03, b = 4.95, c = 5.04, d = 1.02. The external and internal acinar pressure are P ext,AC = P PL , and where P AC is the acinar pressure, P PL is pleural pressure, γ is the surface tension of the liquid lining and R ALV represents the mean radius of the interior surface of the alveolar lining fluid inside the acinus. At this stage, alveoli are assumed to be spherical. From Equations (9) and (10), Equation (2) can be rewritten as

IV. Parenchymal Tethering
We utilized the parenchyma model from Ryans et al. (2019) to describe the tethering pressure. In that model, a cylindrical "hole" is surrounded by an annular region of parenchyma. An airway is laminated inside this parenchymal hole. Without the laminated airway, the equilibrium hole radius is defined as R H so that when the strain in the annular parenchyma is uniform, this yields P ext,AW = P PL . For non-uniform strain, following Lai-Fook et al. (1978), we define the effective parenchyma shear modulus G eff by where R AW is the airway radius, P ext,AW is the peri-airway pressure, P PL is the pleural pressure. From analyses performed by Ryans et al. (2019), we assume the equilibrium hole radius is proportional to the cubic root of the weighted average of acinus volume within a distance of 5R AW,0 from the airway axis, where R AW,0 is R AW at zero transmural pressure. This behavior is consistent with the finite-element model by Ma et al. (2013). We also assumed the equilibrium hole radius R H,FRC is equal to R AW,FRC when the surrounding acinus volume is at the Functional Residual Capacity, V FRC . Since the exact value of FRC is unknown before the simulation, R H,FRC and R AW,FRC are estimated from the isolated airway and acini model. We note that the equilibrium hole radius setpoint at FRC is an assumption that should be experimentally validated and is related to positive strain-deviation of terminal airways, as described in section Parenchymal Tethering and Strain Deviation. From computational simulations as described from Ryans et al., the effective parenchyma shear modulus G eff is empirically determined as a function of mean transpulmonary pressure (Ryans et al., 2019). The transpulmonary pressure and effective parenchyma shear modulus can be expressed as where P TP is average transpulmonary pressure, P ALV is the weighted average of surrounding acinar pressures within a distance of 5R AW,0 from the airway axis,

I. Air Conductance
We assume that all air flows are laminar and fully developed. Airflow throughout the lung model is described using a Poiseuille relationship at each airway. From Poiseuille's law, the conductance of a compliant airway is prescribed as a function of its radius and length: where C AW is the conductance of airway, R in,AW is the radius from the center line of the airway to the interior surface of the liquid lining, µ is the viscosity of air, L is the length of the airway. We acknowledge that this is inaccurate when Reynolds numbers are high, as in the central airways (trachea and mainstem bronchi) and first generations of the model of the lung. We also neglect the loss of pressure at bifurcations.

II. Network Flow Dynamics
In our model, the axial component of the pressure gradient inside each airway is assumed to be a constant. We also assume that all air flows through airways are fully developed and without turbulence. The flow rates at the airway level are determined by where Q AW,i is the flow rate of the airway, C AW,i is the conductance of the airway, P UP,i and P DN,i are the upstream and downstream pressure of the airway, respectively. At the terminal end of the system (acinus), the flow changes the acinus volume as where P AC,i is the acinus pressure, V AC,i is the acinus volume. We assume that there is no pressure loss at the bifurcations. Conservation of flow rate at each bifurcation is satisfied by an equation of the form where P stands for the parent airway, and D stands for the daughter airways. Equations (16)-(18) are applied to each acinus, airway and bifurcation to build the differential-algebraic equation system. We note that we neglect the loss of pressure at the airway bifurcation (Filoche et al., 2015).

I. Conservation of Liquid Phase
We assume the amount of liquid is conserved within each airway/alveolus. Thus the liquid film thickness has to varies with time due to changes in radius of airway/alveolus. A parameter which characterizes this change is defined as where R out is the airway/alveolus radius, R in is the radius from the center(line) to the liquid lining that determines the Laplace pressure drop. The change in the volume of the airways and alveolus therefore modulates the film thickness, which feeds back to the stress balance Equations (7) and (11). We note that the change of airway caliber as well as length is accounted for in the conservation of mass calculation for determining the dynamic film thickness.

II. Surfactant Physicochemical Interactions
The mechanism of surfactant transport was based on the multilayer surfactant behavior analysis done by Krueger and Gaver (2000). Soluble surfactant in the liquid lining of airways/acini exists in the bulk fluid and two surface layers: primary and secondary layers. The primary layer is in direct contact with the air-liquid interface, while the secondary layer is created by collapse of the primary layer and resides between the primary layer and the bulk fluid. Adsorption/desorption can only happen between the bulk liquid and the primary layer. Only the surfactant in the primary layer can reduce the surface tension. The surface tension of water γ 0 = 72 dyne cm . From Krueger and Gaver (2000), the equilibrium surface tension γ ∞ = 22 dyne cm at the equilibrium surfactant concentration Γ ∞ = 3×10 −4 mg cm 2 . The minimum surface tension γ min = 4 dyne cm occurs at the maximum surfactant concentration Γ max = 3.3×10 −4 mg cm 2 . We assumed the relationship between the primary layer surfactant concentration and surface tension is linear from (0, γ 0 ) to (Γ ∞ , γ ∞ ) and from (Γ ∞ , γ ∞ ) to (Γ max , γ min ) with a slope of m and m ′ , respectively. The surface tension equation of state can be expressed as where m = −1.7 × 10 5 dyne · cm mg , m ′ = −5.3 × 10 5 dyne · cm mg . Figure 5 demonstrates an example hysteresis behavior of alveolar unit as determined from Equations (21) to (25) below. In Figure 5, Γ max = 4.2 × 10 −4 mg/cm 2 is the maximum surfactant concentration, γ min = 2 dyne/cm is the surface tension when Γ 1 = Γ max , Γ ∞ = 3 × 10 −4 mg/cm 2 is the surfactant concentration at equilibrium state, γ ∞ = 22 dyne/cm is the surface tension when Γ 1 = Γ ∞ .
As described by Figure 5, the alveolar unit starts contracting from the maximum volume point (top right corner). During the contraction, the surface tension decreases as the monolayer surfactant concentration increases with simultaneous adsorption/desorption occurring between the surface and bulk phases. When the monolayer surfactant concentration exceeds Γ ∞ , the adsorption from the bulk stops. At this time, the monolayer cannot sustain further compression without collapsing to create a multi-layer with surfactant molecules exuded from the primary to the secondary layer as the alveolar unit continues contracting. Desorption between the primary layer and bulk fluid continues as Γ 1 exceeds Γ ∞ and simultaneously the primary layer (Γ 1 ) collapses to transport mass to the secondary layer (Γ 2 ). At the bottom right corner, the primary layer concentration reaches Γ max . The speed of collapse increases by an order of magnitude when Γ 1 > Γ max . Therefore Γ 1 ∼ Γ max as the contraction continues, highly enriching the secondary layer. At the bottom left corner, the alveolar unit starts to expand. The surface tension increases during this expansion due to the decreasing concentration of the primary layer. When the concentration of the primary layer begins to fall below Γ ∞ , respreading starts from the secondary layer to the primary layer, and surfactant molecules are once again adsorbed from the bulk. This hysteresis loop is a key characteristic relating surfactant physicochemical behavior to the stabilizing properties associated with dynamic surface tension in the healthy lung (Notter, 2000).
The governing equations of the surfactant transport are where M 1 is the surfactant mass in the primary layer, M 2 is the surfactant mass in the secondary layer, M b is the surfactant mass in the bulk liquid, A is the surface area of the airliquid interface, j terms represent the transport between bulk and surface layers. Terms for collapse and respreading are defined as where Γ 1 is the primary layer surfactant concentration, Γ 2 is the secondary layer concentration, C collapse = 2 × 10 −4 ( mg s·cm 2 ), C respread = 5×10 −2 ( mg s·cm 2 ), Γ mls = 2.34 mg/cm 2 is the multilayer respreading concentration (corresponds to γ mls = 33 dyne/cm), serve as a switch to turn on the collapse/respreading when Γ 1 is above/below Γ ∞ , and describes the relationship between concentration and speed of transport together with the sigmoidal function The collapse term increases exponentially when Γ 1 is greater than Γ max based on the experimental evidence (Clements, 1957;Horn and Davis, 1975;Schurch, 1995;Krueger and Gaver, 2000). Note that we have modified the collapse and respreading terms from that originally introduced by Krueger and Gaver (2000) in order to stabilize the computations. In the present model we assume that collapse can be initiated (albeit slowly) as soon as the primary layer concentration exceeds Γ ∞ , and this collapse rate increases substantially as Γ 1 increases to Γ max . Likewise, respreading initiates when the primary layer concentration falls below Γ ∞ , and the respreading rate increases substantially near Γ mls .
Terms for adsorption and desorption are defined as (25) FIGURE 5 | Example hysteresis loop of an alveolus. Γ 1 is inversely proportional to the interfacial area. The primary layer begins to collapse when Γ 1 ≥ Γ ∞ , forming a secondary layer between the bulk and primary layer. When Γ 1 < Γ ∞ , the secondary layer dynamically respreads, which stabilizes the surface tension.

Computational Methods
High-performance computing simulations were conducted on Tulane's Cypress supercomputer. The code was constructed in C++ and utilized MPI library. The differential-algebraic equation systems solver used in the code was SUNDIALS IDAS solver with adaptive time stepping (Hindmarsh et al., 2005). The simulation utilizes 8 Intel Xeon E5-2680 v2 CPUs and 128 GB of RAM. Total CPU time for one simulation (180 s of breathing) is ∼5 h.

Overview
In this study, pleural pressure forcing was assumed to be periodic with −5 ≤ P PL ≤ −8.2 cmH 2 O and a frequency of 12 min −1 to simulate tidal breathing (Swan et al., 2012). The inhalationexhalation time ratio is assumed to be to be 1:2. As an initial condition, the dimensionless liquid film thickness ε, and the total amount of surfactant, M total = M 1 + M 2 + M b were prescribed at each airway as well as each acinus when the lung is at its residual volume (RV). In this study, we assumed the minimum alveolar radius equals 100 µm, and the alveolar liquid film thickness equals 0.1 µm at residual lung volume (Bastacky et al., 1985), so ε = 1 × 10 −3 at RV in the acinus components. For airway components, the initial ε is assumed to be 0.02 (Cassidy et al., 1985). We initially apply a uniform bulk surfactant concentration c aiway = 6 mg ml for airway components, and c alveolus = 35 mg ml for acinar components. During initialization, the equilibrium During the simulation, states of each component are recorded to track the evolution of the system. For visualization, a 3-D figure of the half-lung was generated. This figure reflects the space distribution of each airway and acinus as well as their real-time states. Our simulation was performed for several breathing cycles until stationary states were achieved. Figure 6 demonstrates the changes of surface tension in the lung (consisting of airways and acini) at four different time points (a-d) during a breathing cycle. Starting with end-exhalation (a), the surface tension in acini and small airways is as small as 8 dynes/cm, which is far below the equilibrium surface tension (22 dynes/cm) and can be achieved only due to the dynamic transport of surfactant. This low surface tension occurs because the air-liquid interface in each component contracts, increasing the primary layer interfacial surfactant concentration, and decreasing the surface tension (as surfactant also loads the secondary layer). As the lung volume increases during inspiration (b), the surface tension increases, implying a decrease of primary layer surfactant concentration. At the end of inspiration (c), the surface tension attains close to a maximum value in both acini and airways. During expiration (d), the volume of the lung decreases and thus the inner surface area of acini and airways are decreased, reducing the surface tension.
We note that at points a and c, the lung volumes and surface tensions are not respectively the minimum and maximum due to the airway network resistance and the Laplace pressure drop explained in the Pressure-Volume Curve section. Likewise, the lung volumes as well as the surface tensions at points b and d are not identical, even though these represent identical P PL values. The differentials in lung volume and surface tension are related to the phase-lag that occurs in the system, described below.

Pressure-Volume Behavior
The pressure-volume behavior described here reflects both lungs, that is twice the size of our computational domain. Our lung model has a full lung residual volume (RV) of 1.0 L. Assuming homogeneous acinar behavior, following Equation (8) the total lung capacity (TLC) is 6.0 L. The dynamic properties of the lung are determined by pressure cycling and Figure 7 presents the pressure-volume (PV) relationship for the entire lung. Under the given pleural pressure range investigated for normal breathing, the tidal volume (TV) equals 0.82 L, with a functional residual capacity (FRC) of ∼2.30 L. From these simulations of normal breathing, the maximum flow rate at the trachea is ∼0.92 L/s. The average compliance of the lung during tidal breathing is ∼0.25 L/cmH 2 O, while the accepted normal value is 0.2 L/cmH 2 O in human lungs (Harris, 2005). TV and compliance of our model is higher than the normal values because the simulated airway resistance is unrealistically low due to the assumptions of laminar, fully developed flow in "Modeling Implementation" as described in "Limitations". The total airway resistance excluding upper airways ranges from 0.6 to 0.9 cmH 2 O/L/s with lung volume oscillating between 2.30 L (38% of TLC) and 3.12 L (52% of TLC), which is less than the predicted flow resistance at low equivalent lung volumes predicted by Pedley et al. (1970).
Our model presents a typical PV hysteresis behavior, i.e., the compliance of the lung ( volume/pressure) is different between the inspiration and expiration phase. The PV phase lag that is demonstrated in Figure 6 is evident from points a-d in Figure 7. Figure 8 further elucidates the acinar contributions to PV hysteresis behavior-here we present the behavior in terms of local forcing ( P AC = P AC − P PL , represented by the blue curve), and global forcing (P PL , represented by the red curve). The global forcing demonstrates a rounded behavior similar to the full lung in Figure 7. However, the behavior of a specific acinus is best understood from the local forcing, P AC . This analysis demonstrates that the hysteretic behavior is reflected by two factors: a. The flow resistance of airways. As seen in Figure 8A, the hysteresis area of the P AC V curve is much smaller than P PL V, which results from the network airflow resistance. The time lag between the global and local pressure changes caused by the airway network resistance can be observed in Figure 8B, and causes the lung volume to lag the changes in P PL , as demonstrated by Figure 6. b. The Laplace pressure drop. The hysteresis area of the local PV curve ( P AC , represented by the blue curve) results from the behavior of surfactant transport as shown in Figure 5, as described in the next section.

Surfactant Transport
In order to understand the PV relationships in Figure 8, it is necessary to explore the dynamic surface tension behavior at the alveolar and airway levels (Figure 9). As seen in Figure 9A, the surface tension hysteresis behavior of acinus components   Figure 10A. (B) Representative surface-tension vs. normalized surface area of airways demonstrating variations in behavior. The acinar and airway surface tension hysteresis behavior is similar to the isolated model in Figure 5 and contributes to the PV hysteresis behavior of the entire lung. Surface tension variation of airways corresponds to observations in Figure 6.
is similar to the isolated interface model (Figure 5) except that a minimum surface tension plateau does not exist at the lowest acinar volumes. This deviation occurs because the surfactant concentration of the acinus model does not approach Γ max under the breathing conditions imposed by this simulation. In this model, "collapse" refers to the transport of surfactant from the primary to secondary layers. This transport occurs if Γ 1 > Γ ∞ , and occurs at a greater rate as Γ 1 increases, to a point where saturation will occur if Γ 1 ≥ Γ ∞ . In our full-lung model, the rate of interfacial contraction is not large enough to hyper-concentrate the alveolar primary layer to a level that leads to secondary layer saturation. Nevertheless, a significant secondary layer does exist in acinus components. We note that there is a sharp slope change in the local PV curve in Figure 8A. This abrupt change in compliance is due to the secondary layer respreading, which reduces surface tension and work of breathing. This sharp slope change corresponds to the slope change at point d in Figure 9A.
Airways ( Figure 9B) in our model have surface tension hysteresis behaviors similar to acini. Airways exhibit hysteresis areas that varies due to geometrical as well as tethering pressure difference. This corresponds to observations in Figure 6.
We also note that the surface tension increases slightly prior to alveolar and airway expansion (dotted circles), which is likewise observed in the acinar PV relationship (Figure 8A). We attribute this behavior to surfactant over-collapse. In our model, the direction of transport between the primary and secondary layer is only determined by the value of Γ 1 . During most of the phase a-b-c-d in Figure 9A, the primary layer concentration exceeds Γ ∞ , so surfactant transport from the primary to secondary layer occurs. However, as described above, this rate of transport is not large enough to saturate the secondary layer. Nevertheless, the secondary layer is able to respread to the primary layer during the phase d-a in Figure 9A. In phase aand c-d, the surface concentration is dominated by the surface area change, i.e., the surface concentration increases/decreases as the surface area decreases/increases. However, due to the sinusoidal pleural pressure, the change rate of the surface area in phase c-d decreases dramatically, while the collapse term in Equation (22) does not have any abrupt changes. Here, the surface concentration becomes transport-dominated, i.e., the surface concentration decreases as surfactant being exuded from the primary layer. Figure 6 demonstrated the surface tension in our model and shows spatial variation. This variation is due to physiologically based heterogeneity that is inherent in our model. For example, there are geometrical differences that exist from generation to generation based upon the space-filling algorithm used to construct the airway tree (Tawhai et al., 2000). Acini are likewise non-uniform because they fill different-sized regions of the lung depending upon the location of the terminal airways. The local surface tension depends upon the lining fluid/surfactant that is distributed within the alveolar components of the acini. While this is assumed to be a uniform concentration, the local strainfield is nonuniform. The physicochemical interactions from variable strain result in transport behavior that leads to nonuniform heterogeneous surface tensions in this dynamic model. Figure 10A demonstrates the average and standard deviations of the acinar surface tension. This figure demonstrates an upper limit of surface tension near γ ∞ with an insignificant variation. This occurs as the secondary layer respreads, and is indicated by point d in Figure 9A. In contrast, at low surface tension the variation is much greater.
The acinar minimum surface tension is found to be corelated to its residual volume, as shown in Figure 11A. We subdivide the data into two groups: (1) symmetrical acini that are created at the ends of two equivalent daughter terminal airways, and (2) an asymmetrical acinus that develops from a single terminal airway, with the other daughter airway non-terminal (see figure inset). We investigate the linear regressions using the form γ min = aV RV + b. On average, symmetrically produced acini have higher the surface tensions than asymmetric acini. We also explored the influence of central airway-to-acinus path length on the surface tension ( Figure 11B) and found no correlation between these values.
For airways, the surface tension heterogeneity also exists. Figure 10B shows the relationship between airway generation (Horsfield) and maximum/minimum airway surface tension. The points represent the average value of all the airways in that generation in one breathing cycle, and the bars represent the standard deviation. Similar to acini, the variation of the maximum surface tension is insignificant due to the surfactant transport mechanism, while a much greater variation of the minimum exists in each generation. The difference among generations can be found to correlate with the tube law defined by Equation (3) and Table 1.
To investigate the impact of a sigh to enhance surfactant surface-layer sorption, we introduced a sigh-like maneuver by temporarily tripling the magnitude of pleural pressure (from 3.2 to 9.6 cmH 2 O), and doubling the time period (from 5 to 10 s per breath). During this cycle, the maximum lung volume increased to 5.96 L. As shown in Figure 12, this high increment of lung volume temporarily depletes the secondary layer and reduces the surfactant concentration of the primary layer, and allows for greater sorption from the bulk to the interface. Upon exhalation, the primary layer is highly concentrated and significantly reduces the surface tension to very low values (<5 dynes/cm). Simultaneously, transfer to the secondary layer occurs, and this leads to a lower surface tension on the following inspiration compared to that observed during the stationary state of tidal breathing.

Parenchymal Tethering and Strain Deviation
As the red curve in Figure 13A shown, airways in our model has a hysteresis pressure-radius behavior similar to the acini. As discussed in the previous sections, this hysteresis behavior is a result from the flow resistance of the airway network and the surfactant behavior (same mechanics as in acini), and more importantly, the pressure-related tethering force. To describe the tethering effect of the parenchyma, we rewrite Equation (12) and define Strain as where R H is the hole radius of the parenchyma, R AW is the airway radius. If Strain is negative, it means R H < R AW , and the tethering force is compressing the airway. Likewise, if Strain is positive, R H > R AW , and the tethering force is pulling the airway open. In Figure 14, the strain deviation at four different time point of the breathing cycle are shown.
It can be observed in Figure 14 that Strain is positive almost all the time, which indicates that airways are pulled outward by parenchyma. To elucidate this stabilizing force, we investigate the parenchyma effect 2G eff R R and airways' own properties P AW − P PL (t) γ R in,AW (t) in Equation (7) separately. We simulate the system without tethering, i.e., forcing 2G eff R R = 0 in Equation (7). Figure 13A demonstrates the terminal airway radius with and without tethering (R AW, w P , and R AW, w o P , respectively) compared to the size of the hole radius, R H . This figure shows that R H is greater than the size of the terminal airways. This induces an outward tethering stress that increases the airway radius, so R AW, w P > R AW, w o P . This stent-like behavior depends upon the set-point for the equilibrium hole-radius (set at V FRC ) as described in "Parenchymal Tethering." This assumption leads to physiologically reasonable positive straindeviations during normal breathing for nearly all generations of airways, indicating that the parenchymal tethering supports the airway structure in this healthy lung model. We found that reducing the set-point of the equilibrium hole-radius to volumes below V FRC results in negative-strain deviation FIGURE 12 | Acinar (A) and representative airways, (B) surface tension after a deep sigh. The secondary layer respreading stabilizes the surface tension around γ ∞ . In airways, secondary layer depletion can be seen where the surface tension rises rapidly again after being stabilized around γ ∞ for a short period. Compared to the stationary tidal breathing state, the additional surfactant in the primary layer by respreading from the secondary layer and adsorption from the bulk during sigh leads to a lower surface tension at end-expiration. The representative airways are the same as those shown in Figure 9.
FIGURE 13 | (A) Predictions of terminal airway radius with and without parenchymal tethering. Parenchymal hole radius R H (green), airway radius with parenchymal tethering, R AW, w P (red), and predicted airway radius without parenchymal tethering, R AW, w o P (blue). This figure demonstrates the important feature that in this healthy lung model, R H > R AW, w o P . This implies that parenchyma provides an outwardly directed mechanical stress that stabilizes compliant airways. (B) The average maximum and minimum strain deviation of each Horsfield airway generation. Bars represent the standard deviation.
during portions of the ventilation cycle for mid-generation and terminal airways (data not shown). This behavior could be destabilizing because it would reduce the caliber of compliant airways and lead to an enhanced proclivity for airway obstruction. This parenchymal support behavior is consistent with earlier simulations conducted by Ryans et al. (2019). The tethering force is important to stabilize the open airway and maintain air flow, though surfactant deficiency can still lead to airway closure or meniscus obstruction in the unhealthy lung (Notter, 2000;West, 2012).  Figure 13B demonstrates the relationship between airway generation (Horsfield) and strain deviation Equation (26). The points represent the average value of all airways in that generation in one breathing cycle, and the bars represent the standard deviation. Small strain deviations exist in the most compliant airways, and this deviation is due to timedependence in the acinar response related to Equation (17), since the compliance of the airways is approximately equal to the acini. In contrast, for larger airways, significantly positive strain deviation exists due to several factors: (1) a compliance differences between the acinus and airway, (2) the pressure differences associated with large path-lengths between the large airways and surrounding acinus, and (3) the time-dependence in the acinar response.

Limitations
While this model provides a useful tool for simulating the multiscale interactions of the lung, there are significant limitations due to the assumed idealized conditions. For example, we neglect gravitational effects and assume a uniform pleural pressure distribution. As such, this model cannot yet simulate ventilation dependency that is linked to orientation. This model also neglects smooth muscle biomechanics that may be neuromodulated. Furthermore, we neglect important pressure losses due to turbulence and entrance-flow effects (Filoche et al., 2015). Our simulations of normal breathing induce a maximum flow rate of ∼0.92 L/s in the trachea. The Reynolds number is large (Re > 2,000) in airways of generations n ≤ 3 and the Womersley parameter indicates unsteady flow in generations n ≤ 4. Furthermore, we have assumed fully developed flow, which does not occur near bifurcations. Thus, we can expect that these regions will have greater pressure-drops than those predicted by this model. We also assume that each airway and acinus is initialized with the same dimensionless film thickness ε and a uniform surfactant concentration, and we do not allow for inter-airway or airway-to-acinar transport. Since this is a model of a healthy lung, we do not incorporate interfacial instabilities that can create plugs or regions of atelectasis. Despite these limitations, this multi-scale model a significant step in the development of a multi-scale model of the entire lung, and is capable of reproducing physiologically realistic behavior.
It should be noted that there are a large number of parameters in our system. Some of these are simple physical constants such as viscosity, equilibrium surface tension, adsorption/desorption rates, etc. that can be determined independently. Other parameters include physiological model constants related to the tube law, the set-point for tethering between parenchyma and airways, airway parent-daughter branch angles, etc. Some of these are better known than others, and could be patient specific. These might be considered "tunable" in the sense that they can be subject to revision when more information is available. With present technology, we cannot hope to create patientspecific geometric models of all branching airways and alveoli, lining fluid distributions and regionally dependent surfactant physico-chemical behavior. Despite the inability to define this type of specific information, when fully developed for healthy and pathological states, we will seek to determine the sensitivity of physiological behavior to key measurable parameters. This could provide insight into pathophysiological behavior and system interdependency that is not available from simpler systems. With this understanding, the step toward patient-specific analyses will become clearer.

CONCLUSIONS
We have developed and investigated a computational 3-D lung model that simulates multi-scale pulmonary interactions under healthy conditions. This model includes the airflow between airways and acini, surfactant transport in the liquid lining, and parenchymal tethering between the acini that surrounds airways. As shown above, these simulations demonstrate physiologically realistic macroscale PV relationships and predict microscale strain distributions that deviate from the uniform strain state (Ma et al., 2013). This deviation occurs due to ventilationinduced airflow pressures and non-equilibrium Laplace pressures resulting from surfactant physicochemical interactions (West, 1977(West, , 2012Krueger and Gaver, 2000;Notter, 2000). In the future, this model will provide a baseline for the study of pathological conditions; for example, airway closure and liquid plug movement that may be important contributors to the etiology of ARDS and ventilator-induced lung injury. Simulations of pathophysiological conditions may foster the development of improved methods of ventilation (Kollisch-Singule et al., 2014;Gaver et al., 2020).

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

AUTHOR CONTRIBUTIONS
HM, HF, and DG conceived, designed research, analyzed data, interpreted results of experiments, prepared figures, and drafted manuscript. HM and HF performed experiments. HM, HF, DH, and DG edited, revised manuscript, and approved final version of manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by NSF CBET 1706801 and NIH 1R01HL142702. Portions of this research were conducted on Tulane's Cypress supercomputer.