Journal of Applied Physiology
Vol. 82, No. 6, pp. 2036-2044, June 1997
GAS EXCHANGE, MECHANICS, AND AIRWAYS

A finite-element model of oxygen diffusion in the pulmonary capillaries

Andreas O. Frank, C. J. Charles Chuong, and Robert L. Johnson

Biomedical Engineering Program, University of Texas at Arlington, Arlington 76019; and Department of Internal Medicine, University of Texas Southwestern Medical Center at Dallas, Dallas, Texas 75235

ABSTRACT

INTRODUCTION

METHODS

RESULTS

DISCUSSION

ACKNOWLEDGEMENTS

FOOTNOTES

REFERENCES

ABSTRACT

Frank, Andreas O., C. J. Charles Chuong, and Robert L. Johnson. A finite-element model of oxygen diffusion in the pulmonary capillaries. J. Appl. Physiol. 82(6): 2036-2044, 1997.We determined the overall pulmonary diffusing capacity (DL) and the diffusing capacities of the alveolar membrane (Dm) and the red blood cell (RBC) segments (De) of the diffusional pathway for O₂ by using a two-dimensional finite-element model developed to represent the sheet-flow characteristics of pulmonary capillaries. An axisymmetric model was also considered to assess the effect of geometric configuration. Results showed the membrane segment contributing the major resistance, with the RBC segment resistance increasing as O₂ saturation (SO₂) rises during the RBC transit: RBC contributed 7% of the total resistance at the capillary inlet (SO₂ = 75%) and 30% toward the capillary end (SO₂ = 95%) for a 45% hematocrit (Hct). Both Dm and DL increased as the Hct increased but began approaching a plateau near an Hct of 35%, due to competition between RBCs for O₂ influx. Both Dm and DL were found to be relatively insensitive (2~4%) to changes in plasma protein concentration (28~45%). Axisymmetric results showed similar trends for all Hct and protein concentrations but consistently overestimated the diffusing capacities (~2.2 times), primarily because of an exaggerated air-tissue barrier surface area. The two-dimensional model correlated reasonably well with experimental data and can better represent the O₂ uptake of the pulmonary capillary bed.

finite-element method modeling; hematocrit; plasma protein concentration

INTRODUCTION

OXYGEN UPTAKE WITHIN THE LUNGS, from the alveolar air space to a heme-binding site on a hemoglobin molecule inside of a red blood cell (RBC), is dictated by the diffusion characteristics of this pathway and the chemical reactions within the RBC. Conceptually, the total resistance for the oxygen uptake in this pathway (1/DL, where DL is lung diffusing capacity) can be expressed as the algebraic sum of that, due to a membrane segment (1/Dm, where Dm is membrane diffusing capacity) and an RBC segment [1/( · Vc), where  is the specific rate of gas uptake by RBCs and Vc is the pulmonary capillary blood volume]. This was originally defined by the Roughton and Forster equation (18) written as

(1)

where the apparent DL is in ml · min¹ · mmHg¹, Dm is in ml · min¹ · mmHg¹,  in blood with a normal hematocrit (Hct) is in ml · min¹ · mmHg¹ · ml blood¹, and Vc is in ml. The product of  and Vc is referred to as the RBC diffusing capacity (De). The diffusive transport across the membrane segment, which accounts for both the blood-gas tissue barrier and the plasma fluid, can be mathematically described as a simple passive diffusion process. However, oxygen transport within the RBC segment involves passive diffusion of oxygen, as well as binding of oxygen to hemoglobin and diffusion of oxyhemoglobin, i.e., facilitated diffusion of oxygen.

The lumped-parameter representation of Eq. 1 was a conceptual milestone, which allowed Dm and Vc to be quantified from experimental determination of DL for CO at different oxygen tensions (18). Such data could be translated into oxygen diffusing capacities. Equation 1, however, does not address the spatially distributed nature of the oxygen gas transport process. Federspiel (8) applied a finite-difference numerical method to calculate the DL by assuming spherical RBCs uniformly spaced inside a cylindrical capillary tube. However, capillaries in the lung are sandwiched between sheets of alveolar membrane as described by sheet flow (9); thus the capillary tube configuration may not accurately reflect the air-tissue barrier across which gas transport must occur in the lungs and may lead to an overestimation of diffusing capacities. Wang and Popel (22) studied the oxygen release in the systemic microcirculation using a finite-element model considering various RBC shapes. Both studies considered the transient diffusion of oxygen and oxyhemoglobin within the RBC coupled with the chemical reaction between oxygen and hemoglobin.

In this paper, we present a model based on the finite-element method (FEM) that describes the transient oxygen transport in the pulmonary capillaries from the alveolar air space to the RBCs. We considered a modified two-dimensional (2D) geometry with parachute-shaped RBCs in a parallel-sided channel, which incorporates the available surface area at the air-tissue barrier and the RBC wall. This modified 2D geometry was used to represent the sheet-flow configuration of the pulmonary capillaries (9). To assess the effects of the capillary geometry, an axisymmetric model was also considered. The contributions to the diffusional resistance imposed by the membrane and RBC segments were determined and expressed as their respective diffusing capacities; from these, the total diffusing capacity was calculated. Effects of varying Hct and varying plasma protein concentrations on diffusive transport were also examined.

METHODS

Geometric model of a typical capillary segment. A modified 2D model with parachute-shaped RBCs in a parallel-sided channel was used to represent the sheet-flow characteristics (9) of pulmonary capillary blood flow. The model incorporated the available surface area for gas transport at the air-tissue barrier and the RBC wall. The model geometry consists of a cross section through the longitudinal axis of a typical pulmonary capillary segment (100 µm length) containing a variable number of equally spaced RBCs depending on the segmental Hct (Fig. 1A).

The parachute shape of the RBC was digitized from a photograph of Skalak and Branemark (20), fitted with cubic splines, and then used as the cross section for all the RBCs in the capillary segment (Fig. 1A). Human RBCs are known to have a mean volume of ~97 µm³ and a mean surface area of ~137 µm² (10), and it is important to incorporate these values into the model, since they can greatly affect the amount of transient gas transport. The 2D planar model was thus modified to have an effective depth of 4.92 µm, resulting in an RBC volume of 103 µm³ and an effective surface area of 125 µm². Note that the effective surface area considers only that which is used for gas transport at the RBC perimeter along the peripheral surface and does not include the area of the front or back faces of the RBC.

Assuming that the RBCs are equally distributed within a pulmonary capillary bed, we can study the oxygen diffusion characteristics by examining a typical unit segment, i.e., one RBC in its surrounding tissue-plasma barrier and the alveolar air. Additionally, we have taken advantage of the symmetry with respect to the x-axis and only modeled one-half of such a typical unit segment (Fig. 1A). Therefore, the model consists of three different regions (Fig. 1B): 1) the blood-gas tissue barrier, 2) the plasma fluid within the capillary, and 3) the RBC, each with respective gas-diffusive properties. The effect of varying Hct was studied by adjusting the "unit length" of the capillary segment, which effectively changes the volume of plasma fluid but does not alter the volume of the RBC within the capillary segment (Fig. 1B). In this study, the Hct was varied from 10 to 50% by increments of 5%.

Axisymmetric representation. By using the same shape as the 2D model (Fig. 1B), an axisymmetric model was also constructed so that the effects of capillary geometry could be assessed. Specifically, the differences in diffusing-capacity parameters could be evaluated, resulting from differences in the surface area available for gas transport and geometric shape factors. For the axisymmetric model, the y-axis is now interpreted as the radial direction (Fig. 1, A and B). As with the 2D planar configuration, the unit length of the diffusion space was adjusted to simulate the effect of varying Hct. All of the radial dimensions were maintained to be the same as those in the 2D model to allow for a direct comparison of the geometric configurations. The resulting volume and the surface area of the RBCs for the axisymmetric case were 103 µm³ and 133 µm², respectively.

Passive diffusion in the membrane segment. Due to the low Peclet number in pulmonary capillary blood flow, we have neglected the convective transport (1, 8, 13). When a reference frame moving with the RBC is used and the differences in the plasma fluid and RBC traveling speed are neglected, the passive diffusive transport in the membrane segment (blood-gas barrier and plasma fluid) can be described by a simple diffusion process written as

(2)

where t is time, ² is the Laplace operator (= ²/x² + ²/y² + ²/z²),  is the solubility coefficient for oxygen (assumed to be constant), and d_O2 is the diffusion coefficient for oxygen.

Facilitated diffusion within the RBC segment. The diffusive transport inside the RBC is described by a more complicated facilitated diffusion due to the oxygen and hemoglobin interaction. Under typical conditions in the microcirculation, the facilitated diffusion can increase the rate of oxygen transport almost twofold (13). Thus the passive diffusion of oxygen, the passive diffusion of oxyhemoglobin, and the reaction between oxygen and hemoglobin must be considered, which can be written as

(3)

(4)

where SO₂ denotes the oxygen saturation, [Hb] the total hemoglobin concentration, and d_Hb the diffusion coefficient for hemoglobin. Both Eqs. 3 and 4 include generation terms accounting for the transient changes of oxygen and oxyhemoglobin due to the chemical reaction. Assuming instantaneous equilibrium for the oxyhemoglobin reaction (23), we can combine Eqs. 3 and 4 into one equation written as

(5)

with PO₂ and SO₂ governed by the Hill equation. Equivalently, the diffusive transport within the RBC can be written as

(6)

where dSO₂ /dPO₂ is the differential form of the Hill equation written as

(7)

with P₅₀ denoting the PO₂ at 50% oxyhemoglobin saturation. In the present work, considering normal human blood at 37°C and pH = 7.4, the empirically derived constants of P₅₀ (= 26 mmHg) and n (= 2.7) were used. Note that the diffusion of oxygen within the RBC is now represented by a single partial differential equation with nonlinear coefficients and one primary unknown (Eq. 6). Unlike the membrane segment, which is governed by a linear equation (Eq. 2), the RBC segment is described with nonlinear "diffusive properties," which need to be updated at different time steps to account for the effect of oxyhemoglobin interaction.

The assumption of instantaneous equilibrium within the RBC is valid as long as the speed of the chemical reaction greatly exceeds that for the diffusive transport. This assumption, however, must break down just inside of the RBC membrane, since the hemoglobin molecule is impermeable to the membrane, whereas oxygen is permeable. The region where deviation from equilibrium occurs has been shown to be only a thin layer inside of the RBC membrane, with the layer thickness varying depending on both spatial position along the membrane and temporal position during the RBC transit (22); i.e., Wang and Popel reported extreme cases for deviations in SO₂ of 10 and 4% at radial positions from the RBC membrane of 4 and 8% of the capillary radius, respectively. Because this thin layer only constitutes a small fraction of the RBC volume and becomes smaller as the RBC transit time decreases (22), the effects on the volume-weighted RBC PO₂ averaged over the entire transient have been neglected in the present work.

Initial and boundary conditions. The initial conditions require that the oxygen distribution (PO₂) for the entire domain be prescribed. This was established by a preliminary steady-state analysis (at t = 0) for the entire domain with the following boundary conditions Thus the initial PO₂ inside the RBC entering the pulmonary capillary segment was assumed to be uniform at 40 mmHg (23). The boundary conditions were maintained for the entire transient at i.e., the alveolar oxygen tension was assumed to be constant at 100 mmHg throughout the entire transient, considering that there is adequate stirring in the alveolar air space. Additionally, periodic boundary conditions were imposed by forcing the PO₂ distribution to match at the left and right edges of the unit length (Fig. 1, A and B) as It should be noted that the PO₂ distribution is continuous across the RBC-plasma interface. For the axisymmetric case, the coordinate pairs {X, Y } are understood to be {X, R} in the above initial and boundary conditions.

Discretized FEM model and solution method. The diffusive transport for the membrane segment (blood-gas barrier and plasma regions) is governed by Eq. 2, whereas for the RBC segment it is governed by Eq. 6. The physical properties used for each different region of the model are summarized in Table 1. With isoparametric formulation, a typical discretized FEM model consists of ~500 bilinear elements with nodal PO₂ as the primary degree of freedom (Fig. 1C). For temporal discretization, a time step of 2 ms was used for the first 100 ms, followed by a time step of 4 ms for the next 200 ms of the analysis. The total transient analysis was run for 300 ms.

Table  1.   Physical properties used in the FEM model Parameter Value Units Source (Ref. No.) Blood-gas barrier   O2 1.4  nmol · cm3 · mmHg1 8   dO2 2.4 × 105 cm2/s 8 Plasma   O2 1.4  nmol · cm3 · mmHg1 8   dO2 for protein level     Low (5 g/100 ml) 2.49 × 105 cm2/s 24     Average       (6.9 g/100 ml) 2.40 × 105 cm2/s 24     High (10 g/100 ml) 2.25 × 105 cm2/s 24 Red blood cell   O2 1.40 nmol · cm3 · mmHg1 8   dO2 2.4 × 105 cm2/s 8   [Hb] 2.0 × 104 nmol/cm3 8   dHb 1.4 × 107 cm2/s 8 FEM, finite-element model; O2, solubility coefficient for oxygen; dO2, diffusion coefficient for oxygen; [Hb], hemoglobin concentration; dHb, diffusion coefficient for hemoglobin.

Governing Eqs. 2 and 6 were solved by using the Newton-Raphson iteration technique, with the time transient portion solved by applying the general form of the trapezoid rule (2). The analysis is to find, at each time step, the PO₂ distribution at all nodal points that satisfy both the initial and boundary conditions. The FEM software ANSYS 5.0A (Swanson Analysis System, Houston, PA) running on a DECstation 5000 was used for the analysis.

(8)

(9)

(10)

(11)

RESULTS

11

1

1

11

1

1

11

1

1

11

1

1

11

11

1

1

11

11

1

1

11

11

1

1

11

1

1

DISCUSSION

Variability in transient Dm and De. The progressive importance of the RBC segment in total DL can be seen in the distribution of oxygen flux through the RBC wall membrane (Fig. 5). In an early stage of RBC transit (SO₂ = 75%), with the low SO₂ in the RBC, oxygen diffuses across the RBC membrane with little resistance. Oxygen enters the RBC most rapidly at its wall membrane closest to the capillary surface, since it is the path offering the least resistance (i.e., maximal PO₂ gradient). Much less flux is seen at the concave part of the RBC, which is relatively hidden from the alveolar surface. This suggests that the plasma fluid in this vicinity does not contribute to Dm and the utilization of the plasma fluid (membrane segment) is not uniform. At 85% SO₂, due to the gradual increase in the relative resistance of the RBC segment with respect to the membrane segment, the distribution of oxygen flux over different regions of the RBC wall membrane becomes more uniform as compared with an SO₂ of 75%. Similarly, the spatial distribution of oxygen flux within the plasma becomes more uniform. Finally, at 95% SO₂, total flux into the RBC becomes so low that there is no regional preference in oxygen flux through the RBC wall membrane or plasma fluid.

Thus there is a gradual progression during the RBC transit toward a more uniform distribution in the oxygen flux across the RBC wall membrane and in the utilization of the plasma fluid with respect to oxygen transport. This results in the gradual decrease in Dm (Fig. 2), since, effectively, the diffusion of oxygen from the air-tissue barrier into the RBC takes a longer path because of the increased regional resistance of the RBC. Therefore, the regional changes in De occurring throughout the transient cause the utilization of the RBC wall membrane and plasma fluid to be altered, which affects Dm.

Comparison with Dm measurements. To compare with available experimental measurements, we calculated the mean Dm at the total lung volume by using the following extrapolation

(12)

where m,RBC is the mean Dm per RBC calculated by the FEM model throughout the transient, and VRBC is the FEM model volume of a RBC (103 µm³). Note that m,RBC is an implicit function of Hct. Using experimentally determined data for Vc of 100.3 ml at a resting state (cardiac index ~7 l/min) and 162.8 ml at an exercise state (cardiac index ~15.5 l/min) (15), we approximated m,total lung for the FEM model results according to Eq. 12. Values in m,total lung were obtained at resting and exercise states for both the 2D and axisymmetric model configurations (Fig. 6). These extrapolated Dm values for the total lung showed a linear increase up to a Hct of ~35% before gradually approaching a plateau. Measurements of Dm,CO, based on the rebreathing technique (15), were multiplied by 1.24 to yield estimated values for Dm,O₂, shown as constant values in Fig. 6. Shaded regions represent ranges from a mean resting state of 54.1 ml · min¹ · mmHg¹, to a moderate exercise state of 69.9 ml · min¹ · mmHg¹, and to a peak exercise state of 115.7 ml · min¹ · mmHg¹ (15). A visual comparison indicates that the 2D planar FEM model predicted total Dm,O₂ within a reasonable range, more so for the resting state.

Our 2D FEM model has slightly overpredicted Dm,O₂ because it assumed an ideal condition of uniform RBC distribution and spacing within the capillaries, which offers the optimal use of the membrane segment, i.e., highest values in Dm (5). It is known that there is spatial and temporal fluctuation in both RBC distribution and spacing in the pulmonary capillary bed under physiological conditions. For the resting state, the corresponding Hct was in the range of 18-30% (Fig. 6), somewhat lower than the physiological values of 28~37% (3, 7) because of the overprediction in Dm,O₂. For the exercise state, this source of error is further exaggerated, and the deviation in the estimated Hct is probably even larger.

2D vs. axisymmetric model. For all cases with axisymmetric configuration, mean diffusing capacities of 2.0-2.4 times their 2D counterparts were obtained. Factors contributing to the difference in diffusing capacities between geometric configurations can be separated into two categories: 1) the available surface area for gas transport, and 2) the geometric shape of the diffusion space. For the axisymmetric case, the available surface area at the air-tissue barrier was 2.29 times the 2D case considering the same Hct; also, there was a slightly higher RBC surface area (1.06 times). Note that the unit length (Fig. 1, A and B) for the axisymmetric case was necessarily different from the 2D so that we could match their Hct values. The increased surface area, primarily that at the air-tissue barrier, was responsible for most of the exaggerated diffusing capacity parameters. In contrast, for the axisymmetric configuration, the converging effect of the radially inward flux of oxygen would counteract the transport because of the progressively reducing area associated with the geometric shape. The effects due to geometric shape between the two configurations can be quantified in terms of Nusselt number (see below).

The calculated mean Dm values for the total lung derived from the axisymmetric cases were found to be significantly higher than the measured values, whereas the 2D cases more closely corresponded to the measurements (Fig. 6). This suggests that the 2D model is more representative of the oxygen transport in the pulmonary capillary bed. Our comparisons demonstrate that the geometric representation of the diffusion space needs to be carefully incorporated for accurate assessment of the diffusing capacity parameters.

Comparison of Nusselt numbers with other models. To compare with other mathematical models, we calculated the Nusselt number (Nu), a measure of mass transport conductance, defined as

(13)

where RA denotes one-half of the alveolar septum thickness for the 2D case and the radius of the alveolar septum for the axisymmetric case. Calculated Nusselt numbers were added to Fig. 14 of Ref. 14 for comparison (Fig. 7A). It should be noted that all of these previous works on systemic microcirculation utilized cylindrically shaped diffusion spaces with an axisymmetric geometry and the values presented were for a Hct of 25% and hemoglobin saturation of 70%. They were derived to represent oxygen delivery to the tissues, as opposed to oxygen uptake by the lungs; thus the hemoglobin saturation is outside of the range represented by the present model. For comparison purposes, we calculated Nusselt numbers of the present models for a hemoglobin saturation of 75% and a Hct of 25%. The results for the present models compare reasonably with the previous works (Fig. 7A) (12, 16, 17, 19, 22).

The Nusselt number, written in terms of flux (flow/area), represents the ratio between the rates of actual mass transfer and the diffusive transfer. For the axisymmetric case, the Nusselt number was found to be 20% lower than its 2D counterpart at a Hct of 25%, suggesting a higher efficiency with the latter configuration, which is to mimic the sheet-flow configuration. Note that the greater diffusing capacity parameters (~2.2 times) for the axisymmetric configuration are primarily due to its exaggerated surface area at the air-tissue barrier. Effect of surface area is removed in the calculation of Nusselt number, since it is based on flux and not flow (flux × area). Thus the differences in Nusselt number between the axisymmetric and 2D cases are due to the converging effect with the former, i.e., oxygen molecules are competing for the progressively reducing surface area available.

As with the diffusing capacity parameters, the Nusselt number from the present work is seen to be strongly dependent on Hct. Calculated mean Nusselt numbers in terms of Nu  Nu_Hct=25 were plotted as a function of Hct with data from previous works (12, 16) in Fig. 7B. They reveal the dependency of the Nusselt number on Hct for various-sized capillaries. The data from Groebe and Thews (12), based on oxygen delivery with axisymmetric capillaries of a 5.5-µm diameter, show a stronger dependency on Hct as compared with the present models (capillary diameter of 8 µm). On the other hand, data from Nair (16) represent larger vessel diameters (20-100 µm), showing a much smaller dependency on Hct. These comparisons suggest that the Nusselt number (oxygen mass transfer conductance) becomes more strongly dependent on Hct as the capillary vessel size decreases, and for small capillary vessels the Hct can have a significant effect on oxygen transport.

Comparison with RBC conductance . Using the mean e per RBC throughout the entire transient, we calculated an equivalent  value of 5.9 ml · min¹ · mmHg¹ · ml whole blood¹ for an Hct of 45%. This was obtained on the assumption of infinite reaction velocity between oxygen and hemoglobin. To calculate it, we have considered a transient analysis from 0 to 80% SO₂ to simulate the experimental conditions of rapid-reaction techniques (21, 24), which indicated nearly constant  values for this SO₂ range. Above 80% saturation, the RBC conductance  progressively falls as the saturation rises due to the progressive decrease in the number of binding sites provided by the last heme on the hemoglobin molecule (21). Measured values of  have been reported to be 2.8 (21) and 3.8 ml · min¹ · mmHg¹ · ml whole blood¹ (24). However, the rapid-reaction technique may never completely eliminate the unstirred plasma layer around the RBC and, hence, tends to underestimate the  value; whereas the present work tends to overestimate the  value due to the assumed infinite reaction velocity.

Neglect of convective diffusion. Our model neglected convective diffusion of oxygen within the pulmonary capillary. This assumption was made with the consideration that the effect of convective diffusion is small compared with that of conductive diffusion. Two possible mechanisms could contribute to the convective diffusion within the capillary bed: 1) relative motion between the RBC and the blood plasma, and 2) relative motion between the RBC membrane and its interior contents.

Convection due to the former mechanism would lead to increased oxygen conductance of the membrane segment, which could reduce the role of the membrane segment as the limiting resistance of the pulmonary diffusion process (13). Flow velocity at capillary is relatively low (~0.2 cm/s), with a low Reynolds number (~0.001) (4). The relative motion of plasma fluid with respect to the RBC, particularly that at the cell-endothelial gap, is known to be nonuniform with a complex circulating pattern associated with the propulsion of the RBCs (4). The plasma fluid moves at a speed lower than that of the RBC as a consequence of the no-slip velocity boundary condition at the capillary wall. With an effective Peclet number of <1 (8), the enhanced mixing due to the stirring motion in plasma should not appreciably affect the diffusive transport of dissolved gasses such as oxygen (1), especially for higher Hct values (13).

The second mechanism considers the relative motion between the RBC membrane and the cytoplasm within, known as "tank treading" (11). It occurs within the RBC segment of the diffusion space and would thus tend to increase the RBC conductance. An increase of the RBC conductance would tend to decrease its effect on either DL or Dm, due to its already relatively high value.

Summary. We determined diffusing capacities (DL, Dm, and De) for oxygen using a 2D FEM, developed to represent the sheet-flow configuration of pulmonary capillaries. Results showed the membrane segment contributing the major resistance, with the RBC segment resistance increasing as SO₂ rises during the RBC transit. Both Dm and DL increased as the Hct was increased but gradually approached a plateau as the Hct exceeded 35%. Both Dm and DL were found to be relatively insensitive to changes in plasma protein concentration. Axisymmetric results showed similar trends to the 2D model for all Hct and plasma protein concentrations but consistently overestimated the diffusing capacities by ~2.2 times, primarily due to the exaggerated air-tissue barrier surface area. The 2D model correlated reasonably well with experimental data and can better represent the oxygen uptake at the pulmonary capillary bed.

ACKNOWLEDGEMENTS

The authors thank the reviewers for very constructive comments and suggestions.

FOOTNOTES

Address for reprint requests: C. J. C. Chuong, Biomedical Engineering Program, PO Box 19138, 501 W 1st St., EL-220, Univ. of Texas at Arlington, Arlington TX 76019.

Received 22 March 1996; accepted in final form 14 February 1997.

REFERENCES