Pulmonary Oxygen Exchange in a Rhythmically Expanding–Contracting Alveolus–Capillary Model

: Pulmonary gas exchanges are vital to human health, and disruptions to this process have been associated with many respiratory diseases. Previous gas exchange studies have predominately relied on whole-body testing and theoretical analysis with 1D or static models. However, pulmonary gas exchanges are inherently a dynamic process in 3D spaces with instantaneous interactions between air, blood, and tissue. This study aimed to develop a computational model for oxygen exchange that considered all factors mentioned above. Therefore, an integrated alveolus–membrane–capillary geometry was developed with prescribed rhythmic expansion/contraction. Airﬂow ventilation, blood perfusion, and oxygen diffusion were simulated using COMSOL. The temporal and spatial distribution of blood ﬂow and oxygen within the capillaries were simulated under varying breathing depths and cardiac outputs. The results showed highly nonuniform blood ﬂow distributions in the capillary network, while the rhythmic oscillation further increased this nonuniformity, leading to stagnant blood ﬂow in the distal vessels. A static alveolus–capillary geometry underestimated perfusion by 11% for normal respirations, and the deviation grew with breathing depth. The rhythmic motion caused a phase lag in the blood ﬂow. The blood PO 2 reached equilibrium with the alveolar air after traveling 1/5–1/3 of the capillary network. The time to reach this equilibrium was signiﬁcantly inﬂuenced by the air–blood barrier diffusivity, while it was only slightly affected by the perfusion rate. The computational platform in this study could be instrumental in obtaining reﬁned knowledge of pulmonary O 2 exchanges.


Introduction
Normal oxygen exchanges in the lungs are critical to human health. Interruptions to pulmonary gas exchanges have been frequently observed in patients with respiratory distress, leading to ventilator usage and high rates of mortalities. The oxygen (O 2 ) is inhaled into the pulmonary alveoli, enters the blood to bind with the hemoglobin in red blood cells, and is then transported to body tissues, where metabolism occurs. Within alveoli, O 2 moves across thin respiratory membranes, which consist of alveolar squamous cells, capillary endothelial cells, and basement membranes [1]. Oxygen exchange occurs mainly through diffusion driven by the gradient of concentration or partial pressure (PO 2 ) [2]. Ambient air is a gas mixture, with each gas contributing independently to the total pressure. Note that only a fraction of air in the alveoli can be refreshed with each breath. This makes the air composition in the alveoli significantly different from that of inhaled air. Typically, the partial pressure in the alveolar was 569 mmHg for N 2 , 100-104 mmHg for O 2 , 50 mmHg for CO 2 , and 47 mmHg for water vapor (i.e., Henry's Law). By comparison, blood enters the lung capillaries with a partial pressure of 40 mmHg for O 2 and 46 mmHg for CO 2 [3].
The oxygen exchange between air and blood takes place across the alveoli-capillary interface. Because the blood volume in pulmonary capillaries at any moment is much lower than the alveolar air volume (i.e., 100 mL blood vs. 3500 mL air), the gas exchange essentially brings the blood oxygen level to the same level as the alveolar air. It is important that the composition of the alveolar air is monitored and adjusted to maintain amenable values [4]. For instance, if the O 2 level drops, deeper breathing will be taken automatically to bring PO 2 back to normal, or vice versa. Factors affecting O 2 exchange include partial pressure gradient, membrane thickness, and air-blood contact area. Respiratory diseases such as pulmonary edema or pneumonia thicken the alveolar membrane and hinder diffusion. Likewise, emphysema reduces the alveolar surface, decreases oxygen exchange, and leads to lower blood oxygen levels [5].
Previous pulmonary gas exchange studies predominately relied on whole-body lung diffusion capacity testing or theoretical analysis with 1D or static models [6][7][8][9]. However, the gas exchange between the alveoli and microvascular is inherently a dynamic process, which is featured by rhythmic motions of the geometry and instantaneous interactions between airflow ventilation, gas diffusion, and blood perfusion. The ventilation-perfusion ratio (V/Q ratio) that we often referred to was based on constant values of air and blood flows in/out of the alveoli and capillaries. However, both the airflow ventilation and capillary perfusion are highly variant in time and space, leading to a time-varying V/Q ratio, which can also significantly vary at different locations of the capillary network. A prominent example of the significance of transient O 2 exchange is the Nobel Prize in Physiology or Medicine 2019, awarded jointly to Kaelin, Ratcliffe, and Semenza "for their discoveries of how cells sense and adapt to oxygen availability" [10]. Moreover, the accumulated blood flow rate can also be influenced by the capillary's oscillating motion, leading to a modified V/Q ratio and subsequent oxygen availability [11].
The objective of this study was to numerically simulate the oxygen exchange in a dynamic alveolus-capillary model. Specific aims include: (1) Develop a computational platform for O 2 exchange in an expanding/contracting geometry that consists of an alveolus, a capillary network, and a membrane between them; (2) Determine the diffusivity of the air-blood interface (membrane) in comparison to empirical O 2 exchange data; (3) Quantify the effects of capillary pressure drop and breathing depth on the temporal blood flow rate; (4) Study the membrane diffusivity and perfusion rate on the temporal-spatial distribution of O 2 in the capillary network.

Alveolus-Capillary Model
A computer model was developed that consisted of one alveolus, a capillary network wrapping the alveolus, and the alveolar membrane between them, as illustrated in Figure 1a. The diameter of the alveolar duct was 70 µm [12,13], and the diameter of the alveolus (i.e., sphere) was 200 µm [14][15][16]. Under normal breathing conditions, the breathing depth, or tidal volume expansion ratio (η = ∆V/V), was 23% [15,17,18], which was equivalent to a radius expansion ratio of 7.15% (i.e., 3 (1 + 0.23) − 1). The capillary had one inlet and outlet, with a network of interconnected vessels between them (Figure 1a). The network consisted of four tori (circular vessels) along the vertical direction, seven half-tori in the x-direction, and one torus in the y-direction. The last torus aligned with the inlet and outlet and was considered the major vessel of the capillary network. Each capillary vessel had a C-ring chape with a hydraulic diameter of 8 µm (Figure 1a

Numerical Methods
COMSOL Multiphysics was used to develop the alveolus-membrane-capillary model geometry. Three materials (i.e., air, tissue, and blood) were applied to the above three domains, respectively (Figure 1a). The oxygen diffusivity in the air and blood was 0.242 × 10 −4 m 2 /s and 2.18 × 10 −9 m 2 /s, respectively [19,20]. No measured oxygen diffusivity through the alveolar membrane (Dm) has been found in the literature [21][22][23]. Note that the traditional O2 diffusion capacity (mL/min) was measured in the whole-body capacity, which cannot be applied to the alveolus-capillary model without introducing additional uncertainties [24][25][26]. In order to fill this gap, an inverse approach was used to identify the Dm value that gave an O2-equilibrium time of 0.25 s (i.e., the time that the blood O2 level reached the alveolar level), which was a well-accepted value for a normal lung under resting conditions [27,28]. Note that the lack of O2 diffusivity through the alveolus-capillary barrier is mainly due to the complex underlying mechanisms. First, O2 dissolves in and diffuses through the surfactant film, the alveolar epithelium, the interstitium, and the capillary endothelium; it then diffuses through the plasma, enters the erythrocyte, combines with hemoglobin, and is transported out of the capillary [29]. Figure 1c shows the computational mesh of the three-compartment models generated in COMSOL. The physics-controlled mesh was used for the alveolus and capillary with normal internal element size and fine prism element in the near-wall region, which were essential to capture near-wall transitions [30,31]. The extra fine mesh was implemented in the membrane (Figure 1c). The simulations were conducted in an AMD Ryzen 3960X 24-Core workstation with 3.79 GHz processors and 256 RAM.
A moving mesh was implemented to the three domains with a prescribed mesh displacement, which first expanded, then contracted, following a sinusoidal waveform (Figure 1b). The period of the waveform was 5 s, and the oscillation amplitude was determined by the breathing depths, as described in 2.1. Four breathing depths were investigated, 11.5%, 23%, 34.5%, and 46%, which gave the radius expansion ratios of 3.70%, 7.15%, 10.38%, and 13.44%, respectively.
Laminar flow was used to simulate airflow ventilation in and out of the expandingcontracting alveolus. An open boundary condition (BC) was specified at the alveolar inlet. The blood flow was also simulated as laminar, considering the mean velocity of approximately 1 mm/s [32,33]. Varying pressure drops (50-300 Pa) were prescribed at the

Numerical Methods
COMSOL Multiphysics was used to develop the alveolus-membrane-capillary model geometry. Three materials (i.e., air, tissue, and blood) were applied to the above three domains, respectively (Figure 1a). The oxygen diffusivity in the air and blood was 0.242 × 10 −4 m 2 /s and 2.18 × 10 −9 m 2 /s, respectively [19,20]. No measured oxygen diffusivity through the alveolar membrane (D m ) has been found in the literature [21][22][23]. Note that the traditional O 2 diffusion capacity (mL/min) was measured in the whole-body capacity, which cannot be applied to the alveolus-capillary model without introducing additional uncertainties [24][25][26]. In order to fill this gap, an inverse approach was used to identify the D m value that gave an O 2 -equilibrium time of 0.25 s (i.e., the time that the blood O 2 level reached the alveolar level), which was a well-accepted value for a normal lung under resting conditions [27,28]. Note that the lack of O 2 diffusivity through the alveolus-capillary barrier is mainly due to the complex underlying mechanisms. First, O 2 dissolves in and diffuses through the surfactant film, the alveolar epithelium, the interstitium, and the capillary endothelium; it then diffuses through the plasma, enters the erythrocyte, combines with hemoglobin, and is transported out of the capillary [29]. Figure 1c shows the computational mesh of the three-compartment models generated in COMSOL. The physics-controlled mesh was used for the alveolus and capillary with normal internal element size and fine prism element in the near-wall region, which were essential to capture near-wall transitions [30,31]. The extra fine mesh was implemented in the membrane (Figure 1c). The simulations were conducted in an AMD Ryzen 3960X 24-Core workstation with 3.79 GHz processors and 256 RAM.
A moving mesh was implemented to the three domains with a prescribed mesh displacement, which first expanded, then contracted, following a sinusoidal waveform (Figure 1b). The period of the waveform was 5 s, and the oscillation amplitude was determined by the breathing depths, as described in 2.1. Four breathing depths were investigated, 11.5%, 23%, 34.5%, and 46%, which gave the radius expansion ratios of 3.70%, 7.15%, 10.38%, and 13.44%, respectively.
Laminar flow was used to simulate airflow ventilation in and out of the expandingcontracting alveolus. An open boundary condition (BC) was specified at the alveolar inlet. The blood flow was also simulated as laminar, considering the mean velocity of approximately 1 mm/s [32,33]. Varying pressure drops (50-300 Pa) were prescribed at the capillary inlet and outlet. The mass transfer was simulated using the diluted species transport model. The oxygen partial pressure (PO 2 ) was 159 mmHg (8.557 mol/m 3 ) at the alveolar inlet (flux BC) and 40 mmHg (2.153 mol/m 3 ) at the capillary inlet. The initial PO 2 in the alveolus was 100 mmHg (5.382 mol/m 3 ). The inhaled oxygen was transported to the alveolar wall through both convection and diffusion, where it continued to diffuse through the membrane and enter the capillary vessels. The blood flow conveyed the O 2 towards the capillary outlet and further distributed it to different regions of the body.

Airflow
Airflow in the alveolus and blood flow in the capillary vessels are shown in Figure 2 at varying instants within one expansion and contraction cycle. The tidal expansion ratio was 23% (i.e., amplitude = 7.15% R 0 ), and the capillary pressure was 100 Pa [15,34]. The airflow streamlines appeared regular in distribution, which did not change much during the entire expansion phase. A similar streamline pattern was also observed during exhalation, only with the flow direction being reversed (lower panel, Figure 2). The mid-plane contour shows the velocity magnitude variations during one respiration cycle. The peak airflow velocity was found at 1.25 s during inhalation and 3.75 s during exhalation.
2, 2, FOR PEER REVIEW 4 capillary inlet and outlet. The mass transfer was simulated using the diluted species transport model. The oxygen partial pressure (PO2) was 159 mmHg (8.557 mol/m 3 ) at the alveolar inlet (flux BC) and 40 mmHg (2.153 mol/m 3 ) at the capillary inlet. The initial PO2 in the alveolus was 100 mmHg (5.382 mol/m 3 ). The inhaled oxygen was transported to the alveolar wall through both convection and diffusion, where it continued to diffuse through the membrane and enter the capillary vessels. The blood flow conveyed the O2 towards the capillary outlet and further distributed it to different regions of the body.

Airflow
Airflow in the alveolus and blood flow in the capillary vessels are shown in Figure 2 at varying instants within one expansion and contraction cycle. The tidal expansion ratio was 23% (i.e., amplitude = 7.15% R0), and the capillary pressure was 100 Pa [15,34]. The airflow streamlines appeared regular in distribution, which did not change much during the entire expansion phase. A similar streamline pattern was also observed during exhalation, only with the flow direction being reversed (lower panel, Figure 2). The mid-plane contour shows the velocity magnitude variations during one respiration cycle. The peak airflow velocity was found at 1.25 s during inhalation and 3.75 s during exhalation.

Capillary Blood Flow
One salient feature of the capillary flow was the oscillations of the capillary volumetric flow rate due to the rhythmic motion ( Figure 2). During expansion, the increased volume of the capillary network drew more blood into the vessels from both openings (inlet and outlet, upper panel, Figure 2). The storage of blood in the stretched capillary network slowed down the blood flow within it. At the center of the capillary network, the compliance-incurred speed decrease could become so significant that the blood flow becomes stagnant. This was especially obvious at t2 = 1.25 s when the expansion speed was the highest. After t2 = 1.25 s, the expansion slowed down, and the stagnant capillaries were gradually recruited. At the end of the expansion (t4 = 2.5 s), all capillaries were recruited and perfused, even though their flow rates could be highly heterogeneous.

Capillary Blood Flow
One salient feature of the capillary flow was the oscillations of the capillary volumetric flow rate due to the rhythmic motion ( Figure 2). During expansion, the increased volume of the capillary network drew more blood into the vessels from both openings (inlet and outlet, upper panel, Figure 2). The storage of blood in the stretched capillary network slowed down the blood flow within it. At the center of the capillary network, the complianceincurred speed decrease could become so significant that the blood flow becomes stagnant. This was especially obvious at t 2 = 1.25 s when the expansion speed was the highest. After t 2 = 1.25 s, the expansion slowed down, and the stagnant capillaries were gradually recruited. At the end of the expansion (t 4 = 2.5 s), all capillaries were recruited and perfused, even though their flow rates could be highly heterogeneous.
The blood flows were further examined in Figure 3 by comparing the perfusion distributions among capillaries within one breathing cycle. The white color in the capillaries means stagnant flow. Again, highly heterogeneous perfusions were observed throughout the entire respiration cycle, and locally stagnant flows occurred most of the time except at the beginning/end of expansion/contraction (t = 2.5, 5.0 s, Figure 3). The two most perfused regions were either close to the inlet or outlet. However, these two regions were not symmetric, with the inlet-proximity region more perfused during expansion and the outlet-proximity region more perfused during contraction ( Figure 3). Note that only the perfused capillaries could transport the ventilated O 2 out of the pulmonary region; without an efflux, a stagnant capillary could only absorb a limited amount of O 2 (i.e., till equilibrium). On the other hand, constant O 2 absorption occurred into an efflux. The absorption rate would increase with the efflux rate (perfusion limited). If the efflux was significantly large to carry away all available O 2 , the absorption rate could change to ventilation limited.
, FOR PEER REVIEW 5 The blood flows were further examined in Figure 3 by comparing the perfusion distributions among capillaries within one breathing cycle. The white color in the capillaries means stagnant flow. Again, highly heterogeneous perfusions were observed throughout the entire respiration cycle, and locally stagnant flows occurred most of the time except at the beginning/end of expansion/contraction (t = 2.5, 5.0 s, Figure 3). The two most perfused regions were either close to the inlet or outlet. However, these two regions were not symmetric, with the inlet-proximity region more perfused during expansion and the outletproximity region more perfused during contraction ( Figure 3). Note that only the perfused capillaries could transport the ventilated O2 out of the pulmonary region; without an efflux, a stagnant capillary could only absorb a limited amount of O2 (i.e., till equilibrium). On the other hand, constant O2 absorption occurred into an efflux. The absorption rate would increase with the efflux rate (perfusion limited). If the efflux was significantly large to carry away all available O2, the absorption rate could change to ventilation limited. Quantitation of perfusion distributions among capillaries within one breathing cycle is presented in Figure 4. In comparison to a perfectly sinusoidal waveform with no phase shift (dashed line), the inlet flow (blue color) had an apparent phase lag, while the outlet flow (red color) had a phase forward shift of approximately 0.55 s (Figure 4a). The velocity at 2.5 s (peak expansion) was higher than that at t = 0.0 s and 5.0 s (minimal volume) because the expansion drew in more blood flow through the capillary inlet. Figure 4b shows the blood flow speeds along the major capillary at nine points at t = 5.0 s. It was clear that the closer to the inlet or outlet, the higher the flow speeds. The perfusion distribution was approximately symmetric right vs. left; however, subtle differences were still discernable. For instance, using the grid line as a reference, it was observed that the flow speeds at points 1 and 3 were slightly higher than that at the corresponding points 9 and 7 ( Figure 4b). Such differences were presumably attributed to the inherent compliance elicited from the expanding-contracting model geometry. Figure 4c shows the flow distribution within the capillary network. By starting from the inlet, there was a cascade of bifurcations of the flow from the upper vessels to the lower vessels, as illustrated by the red clockwise curved arrows. The flows diverted from the upper vessels subsequently merged (but did not mix) with the flow in the lower vessel, which further bifurcated into two streams (zoomed view, Figure 4c). After the half capillary network, a reverse pattern existed for the flows towards the outlet, where a cascade of the flow-merging process occurred from the bottom to the top (counter-clockwise curved arrows, Figure 4c). Quantitation of perfusion distributions among capillaries within one breathing cycle is presented in Figure 4. In comparison to a perfectly sinusoidal waveform with no phase shift (dashed line), the inlet flow (blue color) had an apparent phase lag, while the outlet flow (red color) had a phase forward shift of approximately 0.55 s (Figure 4a). The velocity at 2.5 s (peak expansion) was higher than that at t = 0.0 s and 5.0 s (minimal volume) because the expansion drew in more blood flow through the capillary inlet. Figure 4b shows the blood flow speeds along the major capillary at nine points at t = 5.0 s. It was clear that the closer to the inlet or outlet, the higher the flow speeds. The perfusion distribution was approximately symmetric right vs. left; however, subtle differences were still discernable. For instance, using the grid line as a reference, it was observed that the flow speeds at points 1 and 3 were slightly higher than that at the corresponding points 9 and 7 ( Figure 4b). Such differences were presumably attributed to the inherent compliance elicited from the expanding-contracting model geometry.

Capillary Blood Flow Rate
The effects of the capillary pressure drop on the blood flow rate are presented in Figure 5. Overall, the capillary blood flow rate is proportional to the pressure drop ( Figure  5a). However, the phase shift increases with increasing pressure drop, which can be seen more clearly when replotted in Figure 5b by subtracting the instantaneous velocities from their averages. Here the solid lines represent the inlet velocities, and the dashed lines represent the outlet velocities. Note the approximate symmetry of the flow responses between the right and left vs. the apparent asymmetry between the upper and lower relative to the average (Figure 5b). The latter is because of the net positive flow within the capillary despite the flow velocity oscillations. The difference in the phase shift is caused by the interaction of the flow inertia and compliance, where a longer delay reaching the peak is observed for a larger pressure drop (Figure 5b). Figure 5c shows the normalized blood velocities by dividing the pressure drop. Therefore, the velocities under varying pressures successfully collapsed at the start/end of expansion/contraction (i.e., 0.0, 2.5, and 5.0 s). By contrast, the oscillating amplitudes of the normalized velocities were nonlinear, with a descending order from lower pressures to higher pressures (i.e., from 50 to 300 Pa, Figure  5c).
The effects of the tidal volume (i.e., expansion ratio, η = ΔV/V) on the capillary flows are shown in Figure 6. For a given pressure drop of 100 Pa (Figure 6a), the capillary flow speed remained approximately unchanged at 0, 2.5, and 5.0 s (the end of expansion or contraction). However, the flow oscillation increased with the tidal expansion ratio. Interestingly, all flow profiles under varying tidal volumes met at the same point, i.e., 3.1 s (Figure 6a). The data in Figure 6a were replotted in Figure 6b in terms of the flow fluctuations divided by the tidal expansion ratio η. Normalizing the flows using η nicely collapsed four profiles, and all profiles crossed one point at 2.5 s. Slight differences still existed at the flow peaks (Figure 6b). The larger the tidal volume, the higher the normalized peak and the more phase shift (arrows in Figure 6b).  Figure 4c shows the flow distribution within the capillary network. By starting from the inlet, there was a cascade of bifurcations of the flow from the upper vessels to the lower vessels, as illustrated by the red clockwise curved arrows. The flows diverted from the upper vessels subsequently merged (but did not mix) with the flow in the lower vessel, which further bifurcated into two streams (zoomed view, Figure 4c). After the half capillary network, a reverse pattern existed for the flows towards the outlet, where a cascade of the flow-merging process occurred from the bottom to the top (counter-clockwise curved arrows, Figure 4c).

Capillary Blood Flow Rate
The effects of the capillary pressure drop on the blood flow rate are presented in Figure 5. Overall, the capillary blood flow rate is proportional to the pressure drop ( Figure 5a). However, the phase shift increases with increasing pressure drop, which can be seen more clearly when replotted in Figure 5b by subtracting the instantaneous velocities from their averages. Here the solid lines represent the inlet velocities, and the dashed lines represent the outlet velocities. Note the approximate symmetry of the flow responses between the right and left vs. the apparent asymmetry between the upper and lower relative to the average (Figure 5b). The latter is because of the net positive flow within the capillary despite the flow velocity oscillations. The difference in the phase shift is caused by the interaction of the flow inertia and compliance, where a longer delay reaching the peak is observed for a larger pressure drop (Figure 5b). Figure 5c shows the normalized blood velocities by dividing the pressure drop. Therefore, the velocities under varying pressures successfully collapsed at the start/end of expansion/contraction (i.e., 0.0, 2.5, and 5.0 s). By contrast, the oscillating amplitudes of the normalized velocities were nonlinear, with a descending order from lower pressures to higher pressures (i.e., from 50 to 300 Pa, Figure 5c).  The combined influences from the tidal volume and pressure on capillary flows are shown in Figure 6c. Results under two more pressure drops (50 and 300 Pa) were added to Figure 6c, with long-dashed lines for 50 Pa and dotted lines for 300 Pa (Figure 6c). The flow fluctuations were sensitive to both the pressures and tidal volumes considered in this study, as flow profiles under categories covered the entire velocity range (Figure 6c). As The effects of the tidal volume (i.e., expansion ratio, η = ∆V/V) on the capillary flows are shown in Figure 6. For a given pressure drop of 100 Pa (Figure 6a), the capillary flow speed remained approximately unchanged at 0, 2.5, and 5.0 s (the end of expansion or contraction). However, the flow oscillation increased with the tidal expansion ratio. Interestingly, all flow profiles under varying tidal volumes met at the same point, i.e., 3.1 s (Figure 6a). The data in Figure 6a were replotted in Figure 6b in terms of the flow fluctuations divided by the tidal expansion ratio η. Normalizing the flows using η nicely collapsed four profiles, and all profiles crossed one point at 2.5 s. Slight differences still existed at the flow peaks (Figure 6b). The larger the tidal volume, the higher the normalized peak and the more phase shift (arrows in Figure 6b).  The combined influences from the tidal volume and pressure on capillary flows are shown in Figure 6c. Results under two more pressure drops (50 and 300 Pa) were added to Figure 6c, with long-dashed lines for 50 Pa and dotted lines for 300 Pa (Figure 6c). The flow fluctuations were sensitive to both the pressures and tidal volumes considered in this study, as flow profiles under categories covered the entire velocity range (Figure 6c). As The combined influences from the tidal volume and pressure on capillary flows are shown in Figure 6c. Results under two more pressure drops (50 and 300 Pa) were added to Figure 6c, with long-dashed lines for 50 Pa and dotted lines for 300 Pa (Figure 6c). The flow fluctuations were sensitive to both the pressures and tidal volumes considered in this study, as flow profiles under categories covered the entire velocity range (Figure 6c). As in Figure 6d, normalizing the flow fluctuations by η regressed well into three groups. Moreover, the higher the pressure, the less scattered the η-normalized flow fluctuations.

Alveolar Perfusion and Ventilation
The perfusion volume to the alveolus per breath was calculated by integrating the capillary blood flow rate over one respiration cycle (5 s). Figure 7a shows the calculated alveolar perfusion at varying capillary pressures and breathing depths. The breathing depth of 0% represents a static alveolus-capillary geometry. As alluded to in Figures 5 and 6, the rhythmic oscillation increased the alveolar perfusion more than a static geometry, and the increasing magnitude was sensitive to the breathing depth (Figure 7a). Relative to the static case, the increase was 5.7%, 10.8%, 15.6%, and 20.2% for the breathing depth of 11.5%, 23%, 34.5%, and 46%, respectively (Figure 7b). A similar observation of increased pulmonary perfusion with inspiration was also reported by Hayano et al. [35]. The perfusion per breath was also compared to the alveolar ventilation (V), which was calculated as the product of the alveolus volume and breathing depth, as shown in Figure 7c. It was observed that at η = 11.5%, the ventilation and perfusion matched at 115 Pa (dashed and solid red lines), while at η = 23%, the ventilation-perfusion match was found at 220 Pa (dashed and solid blue lines). Note that the prefusion Q varied with the capillary pressure (x-axis), while the ventilation V was independent of the capillary pressure (Figure 7c). The larger the capillary pressure drop between the vessel inlet and outlet, the higher the perfusion, and the more airflow ventilation was needed to achieve the V-Q matching.
J. Respir. 2022, 2, FOR PEER REVIEW 8 in Figure 6d, normalizing the flow fluctuations by η regressed well into three groups. Moreover, the higher the pressure, the less scattered the η-normalized flow fluctuations.

Alveolar Perfusion and Ventilation
The perfusion volume to the alveolus per breath was calculated by integrating the capillary blood flow rate over one respiration cycle (5 s). Figure 7a shows the calculated alveolar perfusion at varying capillary pressures and breathing depths. The breathing depth of 0% represents a static alveolus-capillary geometry. As alluded to in Figures 5  and 6, the rhythmic oscillation increased the alveolar perfusion more than a static geometry, and the increasing magnitude was sensitive to the breathing depth (Figure 7a). Relative to the static case, the increase was 5.7%, 10.8%, 15.6%, and 20.2% for the breathing depth of 11.5%, 23%, 34.5%, and 46%, respectively (Figure 7b). A similar observation of increased pulmonary perfusion with inspiration was also reported by Hayano et al. [35]. The perfusion per breath was also compared to the alveolar ventilation (V), which was calculated as the product of the alveolus volume and breathing depth, as shown in Figure  7c. It was observed that at η = 11.5%, the ventilation and perfusion matched at 115 Pa (dashed and solid red lines), while at η = 23%, the ventilation-perfusion match was found at 220 Pa (dashed and solid blue lines). Note that the prefusion Q varied with the capillary pressure (x-axis), while the ventilation V was independent of the capillary pressure (Figure 7c). The larger the capillary pressure drop between the vessel inlet and outlet, the higher the perfusion, and the more airflow ventilation was needed to achieve the V-Q matching.

Oxygen Exchange and Transport
3.3.1. Control Case (η = 0.23 and Dm = 2.0 × 10 −10 m 2 /s) Figure 8 shows the oxygen (O2) exchange in the alveolus-capillary model geometry, which involved both convection (alveolar airflow ventilation and blood-flow perfusion) and diffusion (within the alveolar airspace and blood, as well as across the alveolar membrane). In this example (control case), the tidal expansion ratio η was 0.23, and the diffusivity in the air, blood, and membrane were 2.42 × 10 −5 m 2 /s, 2.18 × 10 −9 m 2 /s, and 2.00 × 10 −10 m 2 /s, respectively. The initial O2 partial pressure (PO2) was 159 mmHg at the alveolar inlet, 100 mmHg inside the alveolus, and 40 mmHg at the capillary vessel inlet. Figure 8a shows the O2 temporal variation with a tidal volume of 235 and capillary pressure of 100 Pa at seven points (upper: 3, lower: 4). As expected, the O2 partial pressure at the capillary vessel inlet remained at 40 mmHg throughout the cycle (upper panel, Figure 8a). At the capillary outlet, the PO2 increased with time and reached equilibrium around 0.25 s. Meanwhile, the PO2 inside the alveolus decreased slightly from 100 mmHg to 97 mmHg   Figure 8 shows the oxygen (O 2 ) exchange in the alveolus-capillary model geometry, which involved both convection (alveolar airflow ventilation and blood-flow perfusion) and diffusion (within the alveolar airspace and blood, as well as across the alveolar membrane). In this example (control case), the tidal expansion ratio η was 0.23, and the diffusivity in the air, blood, and membrane were 2.42 × 10 −5 m 2 /s, 2.18 × 10 −9 m 2 /s, and 2.00 × 10 −10 m 2 /s, respectively. The initial O 2 partial pressure (PO 2 ) was 159 mmHg at the alveolar inlet, 100 mmHg inside the alveolus, and 40 mmHg at the capillary vessel inlet. Figure 8a shows the O 2 temporal variation with a tidal volume of 235 and capillary pressure of 100 Pa at seven points (upper: 3, lower: 4). As expected, the O 2 partial pressure at the capillary vessel inlet remained at 40 mmHg throughout the cycle (upper panel, Figure 8a). At the capillary outlet, the PO 2 increased with time and reached equilibrium around 0.25 s. Meanwhile, the PO 2 inside the alveolus decreased slightly from 100 mmHg to 97 mmHg during 0-0.14 s due to the O 2 transport from the alveolus to blood and recovered back to 98.2 mmHg hereafter.

Oxygen Exchange and Transport
J. Respir. 2022, 2, FOR PEER REVIEW 9 during 0-0.14 s due to the O2 transport from the alveolus to blood and recovered back to 98.2 mmHg hereafter.   Figure 8b) and increased quickly with alveolar expansion. Note the drastic color change from blue to red during the first 0.10 s. After 0.10 s, the capillary O2 increased at an ever-decreasing speed and reached full equilibrium around 0.25 s. The PO2 streamlines appeared uniform during 0-0.10 s but became apparently skewed to the vessel inlet after that. This was reasonable considering the more intensified mass transfer near the inlet than in other regions. The lower panel of Figure 8a shows the O2 partial pressure at P4-P1 with descending magnitudes, denoting an increasing O2 diffusion from air to blood. The location of the sampling points P1-P4 is illustrated in Figure 8b with t = 0.1 s.

Effects of Membrane Diffusivity (Dm)
Effects of the membrane diffusivity (Dm) on the oxygen exchange and transport in the alveolus-capillary model are shown in Figure 9. In the diffusivities considered herein (i.e., 2.0 × 10 −10 − 0.125 × 10 −10 m 2 /s), decreasing Dm notably prolonged the time that the oxygen reached an equilibrium between the alveolus and capillary vessels (i.e., the O2equilibrium time, or tO2E, Figure 9a). In comparison to 0.25 s under 2.0 × 10 −10 m 2 /s, the O2equilibrium time was around 0.75 s at 0.5 × 10 −10 m 2 /s (Figure 9a vs. Figure 9b), would be even longer at 0.25 × 10 −10 and 0.125 × 10 −10 m 2 /s (Figure 9c,d). The PO2 at the capillary vessels close to the inlet also decreased with decreasing membrane diffusivities. This resulted in a larger PO2 gradient between the alveolus and vessel that would partially compensate for the reduced O2 exchange from a lower Dm; however, the net O2 flux could still be lower than the capacity that the blood flow perfused, leading to a diffusion-limited O2 exchange (Figure 9c,d). Note that the difference between the red and green lines was due to the time difference in the blood flow from P4 to the vessel outlet.   Figure 8b) and increased quickly with alveolar expansion. Note the drastic color change from blue to red during the first 0.10 s. After 0.10 s, the capillary O 2 increased at an ever-decreasing speed and reached full equilibrium around 0.25 s. The PO 2 streamlines appeared uniform during 0-0.10 s but became apparently skewed to the vessel inlet after that. This was reasonable considering the more intensified mass transfer near the inlet than in other regions. The lower panel of Figure 8a shows the O 2 partial pressure at P4-P1 with descending magnitudes, denoting an increasing O 2 diffusion from air to blood. The location of the sampling points P1-P4 is illustrated in Figure 8b with t = 0.1 s.

Effects of Membrane Diffusivity (D m )
Effects of the membrane diffusivity (D m ) on the oxygen exchange and transport in the alveolus-capillary model are shown in Figure 9. In the diffusivities considered herein (i.e., 2.0 × 10 −10 − 0.125 × 10 −10 m 2 /s), decreasing D m notably prolonged the time that the oxygen reached an equilibrium between the alveolus and capillary vessels (i.e., the O 2 -equilibrium time, or t O2E , Figure 9a). In comparison to 0.25 s under 2.0 × 10 −10 m 2 /s, the O 2 -equilibrium time was around 0.75 s at 0.5 × 10 −10 m 2 /s (Figure 9a vs. Figure 9b), would be even longer at 0.25 × 10 −10 and 0.125 × 10 −10 m 2 /s (Figure 9c,d). The PO 2 at the capillary vessels close to the inlet also decreased with decreasing membrane diffusivities. This resulted in a larger PO 2 gradient between the alveolus and vessel that would partially compensate for the reduced O 2 exchange from a lower D m ; however, the net O 2 flux could still be lower than the capacity that the blood flow perfused, leading to a diffusion-limited O 2 exchange (Figure 9c,d). Note that the difference between the red and green lines was due to the time difference in the blood flow from P4 to the vessel outlet.  Figure 10 shows the effects of the perfusion rate on temporal-spatial oxygen exchange by varying the capillary pressure ∆P (50-200 Pa) with membrane diffusivity Dm being 2.0 × 10 −10 m 2 /s. Relative to the control case (∆P = 100 Pa), decreasing the perfusion rate increased the time difference between the red and green lines (∆t) and increased the PO2 in the vessels close to the inlet (i.e., P1-4, Figure 10a vs. Figure 10b). Conversely, this time difference ∆t was very small at a high perfusion rate (e.g., ∆P = 200 Pa, Figure 10d), reflective of the faster perfusion through the capillary network. The PO2 in the inlet-proximity vessels was also much lower under ∆P = 200 Pa, leading to a larger air-blood gradient in PO2 and a higher oxygen flux from the alveolus to capillary vessels, which was further transported out of the capillary by the blood flow. Note that 1 Pa = 0.0075 mmHg.   Figure 10 shows the effects of the perfusion rate on temporal-spatial oxygen exchange by varying the capillary pressure ∆P (50-200 Pa) with membrane diffusivity D m being 2.0 × 10 −10 m 2 /s. Relative to the control case (∆P = 100 Pa), decreasing the perfusion rate increased the time difference between the red and green lines (∆t) and increased the PO 2 in the vessels close to the inlet (i.e., P1-4, Figure 10a vs. Figure 10b). Conversely, this time difference ∆t was very small at a high perfusion rate (e.g., ∆P = 200 Pa, Figure 10d), reflective of the faster perfusion through the capillary network. The PO 2 in the inletproximity vessels was also much lower under ∆P = 200 Pa, leading to a larger air-blood gradient in PO 2 and a higher oxygen flux from the alveolus to capillary vessels, which was further transported out of the capillary by the blood flow. Note that 1 Pa = 0.0075 mmHg.  Figure 10 shows the effects of the perfusion rate on temporal-spatial oxygen exchange by varying the capillary pressure ∆P (50-200 Pa) with membrane diffusivity Dm being 2.0 × 10 −10 m 2 /s. Relative to the control case (∆P = 100 Pa), decreasing the perfusion rate increased the time difference between the red and green lines (∆t) and increased the PO2 in the vessels close to the inlet (i.e., P1-4, Figure 10a vs. Figure 10b). Conversely, this time difference ∆t was very small at a high perfusion rate (e.g., ∆P = 200 Pa, Figure 10d), reflective of the faster perfusion through the capillary network. The PO2 in the inlet-proximity vessels was also much lower under ∆P = 200 Pa, leading to a larger air-blood gradient in PO2 and a higher oxygen flux from the alveolus to capillary vessels, which was further transported out of the capillary by the blood flow. Note that 1 Pa = 0.0075 mmHg.  The t O2E (time for O 2 to reach equilibrium with the alveolar level) at P4 (red line) appeared to be insensitive to the perfusion rate, which was around 0.2 s for the four perfusion rates considered (Figure 10a-d). This suggested that for a D m of 2.0 × 10 −10 m 2 /s, the O 2 exchange was still perfusion limited for the four perfusion rates, as illustrated by the limited membrane area for mass diffusion (cyan dashed rectangle) at 0.2 s. By contrast, a large portion of the membrane still existed to be recruited for extra O 2 demand (Figure 10b).

Discussion
In this study, high sensitivity was observed in the temporal/spatial O 2 distribution to the alveolus-capillary membrane diffusivity D m . Surprisingly, direct measurements of D m through the membrane have not yet been found; even though the O 2 diffusing capacity has been measured in the whole-body capacity (e.g., in mL/min) [36], it cannot be applied in the microscale model in this study. This lack of O 2 air-blood diffusivity can be partially from the alveolus' tiny dimension that defies current quantification techniques and partially from the limited knowledge regarding the O 2 diffusion pathway that includes not only the alveolar wall but also the red blood cell wall and hemoglobin binding [37,38]. The latter is also demonstrated by the 2019 Nobel Prize in Physiology or Medicine, which recognizes the discoveries related to cells' response to oxygen availability [39]. In order to find an appropriate O 2 diffusivity to start with, preliminary simulations have been conducted. We started with a water-equivalent diffusivity (2.18 × 10 −9 m 2 /s) for the air-blood barrier but obtained an alveolar-capillary O 2 equilibrium time of 20 ms, which was too short in comparison to 0.25 s. After a series of trial-and-error, an O 2 diffusivity of 2.0 × 10 −10 m 2 /s was selected, which gave the alveolar-capillary equilibrium time at 0.25 s under resting conditions (Figure 8). Note that this D m value was merely a result from the attempt of filling the data-missing gap; future in vivo or in vitro measurements are needed to verify/improve the data accuracy.
The results in Figure 9 can have implications in patients with impaired diffusions of the lungs, such as alveolar fibrosis (thickening of the alveolar wall). Other conditions that decrease the diffusing capacity include interstitial edema, sarcoidosis, scleroderma, and pneumonia [40][41][42][43]. In these patients, the capillary PO 2 may not be able to reach the alveolar PO 2 within the breathing cycle, and the oxygen flux from the air to blood is limited by the membrane diffusivity (i.e., diffusion-limited). As a result, an insufficient amount of oxygen can reach the cells and satisfy the demand for cell metabolism [44,45]. By contrast, the results in Figures 8 and 10 represent healthy conditions. With a normal membrane diffusivity, the capillary PO 2 reaches the alveolar PO 2 when the blood is about one-fifth to one-third of the way through the capillary network (Figures 8 and 10). Within this region, diffusion drives the oxygen transfer occurred because of the PO 2 gradient, even though this gradient decreases progressively from the capillary inlet to the equilibrium point. Beyond that region, the oxygen transfer becomes perfusion limited. In circumstances when a larger cardiac output is needed, the diffusion region becomes larger (i.e., microvacuolar recruitment) to meet the increased demand, as illustrated by the gradually increasing PO 2 gradient at the four sampling points (P1-P4) in Figure 10a-d.
The computational model can serve as a platform to study gas exchange under healthy and abnormal breathing conditions. Ideally, the oxygen via ventilation (V) is just sufficient to fully saturate the blood perfusion (Q), as indicated by the cross points in Figure 7c. These two variables, V and Q, act as the major determinants of blood oxygen level [46,47]. In this regard, Figure 8 can be used as a guideline to understand the V-Q match or determine the V needed for a given physical activity (i.e.,~Q). For a typical adult, one-liter blood can hold approximately 200 mL of oxygen, while one-liter air contains around 210 mL of oxygen [48]. Thus, the ideal ventilation/perfusion ratio (V/Q) is around 0.95. On the other hand, a V/Q mismatch can elicit Type 1 respiratory failure [49,50]. A shunt has a zero V/Q, which has perfusion but no ventilation, while a dead space has an infinite V/Q, with ventilation but no perfusion [51]. Exercise leads to a higher demand for oxygen, which regulates the ventilation and perfusion rates, and eventually reaches a new equilibrium [52], a similar transition process as shown in Figure 10a-d.
Surprisingly long computational times were required to simulate O 2 exchange in the computational model developed in this study, even though simplified alveolar and capillary morphologies were used. With a normal mesh for fluid physics in the alveolus and capillary and a refined mesh in the membrane (Figure 1b), the computational time for O 2 exchange during 0-0.75 s was 39-64 COMSOL hours using an AMD Ryzen 3960X 24-Core workstation with 3.79 GHz processors and 256 RAM. The lower the membrane diffusivity, the longer the computation time was required. Several reasons may be responsible for this high demand: a large mesh for the multi-compartment system, multi-physics, multiphases, multi-scale, and fluid-structure interactions. The speeds for the wall motion, airflow ventilation, and blood perfusion were in the range of 0.1-1 mm/s, and the diffusion speeds through the membrane were even lower. It was also observed that COMSOL automatically refined the time step to improve convergence. Time step sizes of 1 × 10 −6 s were often used during simulations, even though a step size of 0.01 s was specified. In addition, this COMSOL-based model was highly CPU demanding; running two test cases would take nearly 100% CPU. In this regard, simulating O 2 exchange in more complex alveoli-capillary geometries can still be prohibitive.
Limitations of this numerical investigation included the lack of comparable measurements, idealized alveolus-capillary model geometry, simplified motion waveforms, constant capillary pressures, and oxygen only. Due to the submicron dimensions of both the alveoli and capillaries, direct in vivo measurements of the gas exchange parameters and associated physical properties were rare (if not existent at all). This lack of data is not surprising considering that the alveolar diameter is 0.2-0.4 mm and the capillary vessel is 8 µm or so, which falls below the resolutions of current characterization or visualization techniques [5,53,54]. The acinar region has complex morphologies with multiple alveoli ventilated via the alveolar duct or pores of Kohn [55]. Even more complicated is the pulmonary microvasculature that wraps the alveolar sac and has sophisticated architectures [56]. Adopting a simplified model as in this study circumvented the prohibitive computational resource that was otherwise required to resolve the above complexities while could still provide insights into the O 2 behaviors in this dynamic, multi-phase (gas, liquid, solid), and multi-physics (convection and diffusion) system. No pathological conditions were considered in this study, such as shunt, partial obstruction, membrane abnormality, etc. [57][58][59]. This study did not consider CO 2 removal. As a result, the interactions between O 2 and CO 2 in the blood, i.e., the Bohr and Haldane effects, were also excluded, which might limit full equilibration between air and blood in certain scenarios [60]. By considering the dead spaces in the respiratory tract, the assumption of the alveolar inlet PO 2 equivalent to the ambient condition should also be improved in future studies.

Conclusions
No previous studies have been reported that considered the dynamic interactions among pulmonary ventilation, membrane diffusion, and capillary perfusion in moving, three-dimensional alveoli-capillary geometry. This study aimed to fill this gap and develop a numerical platform that could evaluate all aforementioned factors, albeit with simplified geometries, kinematics, and physical properties. Temporal-spatial distributions of the blood flow and O 2 were simulated under varying ventilation and perfusion conditions. Specific findings include: (1) Blood perfusion was nonuniform among the capillary vessels, and the geometry oscillation further increased the nonuniformity; (2) A static alveolus-capillary model underestimated the blood flow rate by 11% under resting conditions and increased linearly with the breathing depth for a given capillary pressure; (3) The blood flow had a phase lag (~0.55 s) than the alveolar motion under resting conditions; the phase lag increased with increasing cardiac output; (4) The blood oxygen level reached the alveolar level around 1/5-1/3 of the capillary; oxygen exchange was diffusion driven within this region and was perfusion limited beyond this region; (5) The time to reach the air-blood equilibrium in PO 2 was sensitive to the membrane diffusivity and was relatively insensitive to the blood flow rate; (6) Without measured alveolus-capillary barrier diffusivity for oxygen in the literature, a value of 2.0 × 10 −10 m 2 /s was proposed for normal conditions based on the match to empirical data.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data in this study are available from the corresponding author.