Simultaneous Blood–Tissue Exchange of Oxygen, Carbon Dioxide, Bicarbonate, and Hydrogen Ion

Blood–Tissue Exchange Unit and Nomenclature

The configuration of the blood–tissue exchange (BTEX) unit for the present study of oxygen (O₂), carbon dioxide (CO₂), bicarbonate ( ${\text{HCO}}_3^{-}$), and hydrogen ion (H⁺) transport and exchange in the heart is similar to the geometry of a classical Krogh cylinder model, e.g., see Li et al.⁵¹ as illustrated in Fig. 1. The BTEX unit consists of four concentric regions: The red blood cells (RBC), plasma, interstitial fluid (ISF), and parenchymal cells. Convective flow is F_r (ml g⁻¹ min⁻¹), and is one-dimensional, not accounting directly for intravascular velocity profiles. These regions are separated by the membranes of RBC, capillary, and parenchymal cells, each of which is characterized by its permeability–surface area product, PS. Because the capillary–tissue exchange length, L, that of the BTEX unit, is almost 1 mm,⁶ the model uses partial differential equations^7,9,11 to account for concentration gradients along the length. The anatomical volumes of these regions are denoted by V_rbc, V_pl, V_isf, and V_pc (units ml g⁻¹), their fractional water spaces by W_rbc, W_pl, W_isf, and W_pc (unitless), and corresponding volumes of distribution by $V_{\text{rbc}}^{\prime },V_{\text{pl}}^{\prime },V_{\text{isf}}^{\prime }$, and $V_{\text{pc}}^{\prime }(\text{units}\text{ml}{\mathrm{g}}^{-1})$, which are the anatomical volumes times their fractional water spaces, e.g. $V_{\text{isf}}^{\prime }=V_{\text{isf}}\times W_{\text{isf}}$.

The radial diffusion (perpendicular to the capillary length) of a species (O₂, CO₂, ${\text{HCO}}_3^{-}$, or H⁺) within a region (RBC, plasma, ISF, or parenchymal cells) is treated as if instantaneous because the radial diffusion distances are only a few microns (maximum of about 8 microns for half intercapillary distances in the heart) and the relaxation times for radial gradients are very short, a few milliseconds, compared to the average capillary transit times of a second or so.⁹ This assumption is further justified, for the heart, by there usually being four capillaries around each myocyte. Additional evidence was found in experiments by Deussen and Bassingthwaighte²⁷ in which oxygen was observed to be flow-limited in its exchange with tissue when oxygen consumption was blocked by cyanide in RBC-perfused rabbit hearts, thus illustrating that neither the set of membranes from RBC to mitochondria nor radial diffusion caused any hindrance to exchange (not reported in the Deussen and Bassingthwaighte,²⁷ paper). Given the rapid dissipation of intracapillary radial concentration gradients, the RBC can be well represented as a central column of fluid surrounded by plasma, a situation mathematically similar to Hb being dispersed throughout the plasma.^20,51 (This assumption fails at low hematocrits, Hct, as shown by Hellums³⁶). The RBC and plasma move by plug flow with different velocities; RBC move faster than the plasma, being in the higher velocity central stream. Plasma has zero velocity at the capillary wall, but because the radial diffusion is so fast the velocity profile is considered uniform at the mean velocity of the fluid. This contrasts with arterial velocity profiles because the peak velocity is about 1.6 times the average velocity.⁵ The small vessel (capillary) hematocrits average perhaps three-fourth of large vessel hematocrits, a Fahraeus–Lindquist effect³² with RBC velocities being 30–50% higher than plasma velocities. The blood flow (F_bl) in the capillary is assumed to be steady, so the RBC and plasma flows (F_rbc and F_pl) can be considered as constants. F_bl, F_rbc, and F_pl are in units of ml min⁻¹ (g of tissue)⁻¹.

In the blood space (RBC and plasma), there is both convection and diffusion (axial dispersion). The axial dispersion within a region is characterized by a pseudo diffusion or dispersion coefficient D (units cm² s⁻¹). D_pl characterizes the overall axial dispersion in the plasma by the combined processes of molecular diffusion, axial velocity distribution, flow disturbances, eddy currents, and mixing accompanying RBC rotation. Thus, the value of D_pl is much higher than the coefficient of molecular diffusion in the plasma. It is about 1 × 10⁻⁴ cm² s⁻¹ for O₂.51 We consider it to be approximately the same for CO₂, ${\text{HCO}}_3^{-}$, and H⁺. D_rbc characterizes the overall axial dispersion of RBC (as RBC are suspended particles in the plasma) via all the above processes, except that axial molecular diffusion of Hb is negligible. Thus, D_rbc is about 5 to 10 × 10⁻⁵ cm² s⁻¹.⁵¹ The dispersion coefficient of a solute species within RBC is considered to have the same dispersion coefficient as RBC. The axial transport (along the capillary length, L) of a species in the extra-vascular region (ISF and parenchymal cells) is due to axial molecular diffusion dissipating any concentration gradients, and its effect on the shapes of concentration–time curves is augmented by irregularity in the alignments of the capillaries in sets of neighboring BTEX units.

Offset positioning of starting and ending points of neighboring capillaries results not only in heterogeneity of axial diffusion distances, reducing the effective lengths and reducing the axial concentration gradients,^16,17 but also allows tracer solutes to take a short path from near the entrance of one capillary to the exit region of another, having quantitatively the same effect as enhanced extra-vascular axial diffusion. Accordingly, tracer water dilution curves show a small amount of diffusional shunting from inflow to outflow.⁸ To account for the heterogeneity in D/L, values of D_isf and D_pc are higher than the normal intra-myocardial diffusion coefficients.^56,64 Following the method of Safford et al.⁵⁷ for serial/parallel diffusion through and around cells and accounting for the heterogeneity in D/L, we set the values for dispersion coefficients D_isf and D_pc for all small solutes to 1 × 10⁻⁴ cm² s⁻¹.^27,51 The dispersion smooths axial gradients in each region but does not influence estimation of O₂ consumption rates.

The transmembrane transport (permeation) of a species across the membranes is assumed to be passive and is characterized by the permeability surface area product parameters PS (units ml min⁻¹ g⁻¹), i.e., the product of barrier permeability, P (cm min⁻¹), and surface area, S (cm² g⁻¹). The PS is the barrier conductance for the solute: Permeability P is a diffusion coefficient in the barrier material D_m divided by the barrier thickness d_m (i.e., P = D_m/d_m). Solute species dissolving in the phospholipid bilayer have high Ps. The permeation of ${\text{HCO}}_3^{-}$ occurs through the band-3 anion transporter and that of H⁺ through the proton channels, thus their PSs are about an order of magnitude lower than those for dissolved gases such as O₂ and CO₂ (see Table 1 for values), and yet are still so high that the barrier has almost no influence on the kinetics of tracer transport. (However, the transport of CO₂ is slowed when carbonic anhydrase is inhibited, not because of a change in membrane PS, but because of slow transformation from bicarbonate to CO₂, the form in which most of the flux occurs.)

The consumption of O₂ and production of CO₂ occurs only in parenchymal cells, ignoring oxygenases elsewhere, and is characterized by the “gulosity” parameter ^∘G_pc (units ml min⁻¹ g⁻¹); ^∘G_pc is dependent on intracellular oxygen concentration, assuming Michaelis–Menten kinetics, which is effectively zero-order when P_O2 is high compared to its K_m (which is the case usually) and first-order when O₂ concentration is low. The CO₂ is produced inside the parenchymal cells from oxidative metabolism of glucose, proteins, and fatty acids. The ratio of the rate of CO₂ production to the rate of O₂ consumption is the respiratory quotient RQ: it is 1.0 for glucose metabolism, 0.8 for protein metabolism, and 0.7 for fat metabolism. In myocardium, RQ is normally around 0.8, indicating more fat usage than glucose.

The formation of ${\text{HCO}}_3^{-}$ is governed by the CO₂ hydration reaction which, in aqueous solution, is slow. However, in RBC, ISF, and parenchymal cells this reaction is facilitated by the enzyme carbonic anhydrase in cytosol and on the membranes. The membranes of RBC, capillary, and parenchymal cells are fairly permeable to ${\text{HCO}}_3^{-}$ (through the band-3 transporter), allowing dissipation of gradients. Charge neutrality is attained by the diffusion of chloride ions (Cl⁻) from a region of higher negative charge to a region of higher positive charge (the chloride shift) until the Gibbs–Donnan electrochemical equilibrium conditions are reached^39–41,58:

where R_rbc, R_cap, and R_pc (unitless) are the Gibbs–Donnan ratios across the membranes of RBC, capillary, and parenchymal cells. Under standard physiological conditions, pH in plasma, RBC, ISF, and parenchymal cells is about 7.4, 7.24, 7.2, and 7.1, respectively. Thus, R_rbc = 0.69, R_cap = 0.63, and R_pc = 0.79. The transport of Cl⁻ is not calculated in the model, but its regional equilibrium concentrations can be calculated from these equations.

Dynamic Mass Balance Equations for O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺

The governing equations for each of the four species (O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺) in each of the four regions (RBC, plasma, ISF, and parenchymal cells) of the BTEX unit are formulated by considering the detailed biophysical and biochemical phenomena discussed above with the constraint of mass conservation. The molar concentrations C(x, t), are functions of time (t) and axial distance (x) along the capillary of length (L, cm).

The nomenclature is based on that approved by the American Physiological Society^{10: s}C_r and ^sD_r denote the concentration and axial dispersion coefficient of a species “s” in a region “r”; ^sP S_m denotes the permeability surface area product for a species, s, across a membrane, m; s is one of the species O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺ (pre-superscripted as o, c, b, and h), r is one of the regions RBC, plasma, ISF, and parenchymal cells (subscripted as rbc, pl, isf, and pc), and m is one of the membranes RBC, capillary, and parenchymal cells (subscripted as rbc, cap, and pc). Table 1 lists parameters, representative values, and units.

Equations for Oxygen

Oxygen is transported in the plasma and ISF only in a dissolved form. However, in RBC and parenchymal cells it is transported in the dissolved form as well as in combination with Hb and myoglobin (Mb) as oxyhemoglobin (HbO₂) and oxymyoglobin (MbO₂). The dissolved O₂ in RBC and parenchymal cells is assumed to be in constant equilibrium with the HbO₂ and MbO₂ assuming that the kinetics of O₂ binding to Hb and Mb are fast compared with those of the transport processes. See Dash and Bassingthwaighte²⁶ for detailed reaction schema and kinetic parameters. The equations which follow use the equilibrium constants for the binding. The regional transport equations for convection, diffusion, permeation, and reaction are

where C_HbO2 and C_MbO2 are the concentrations of hemoglobin- and myoglobin-bound O₂ in RBC and parenchymal cells. The concentration of dissolved O₂ in a region r, ^∘C_r, is related to the partial pressure of O₂ in that region, P_O2,r, by Henry’s law: ^∘C_r = α_O2 × P_O2,r, where α_O2 is the solubility of O₂ in water.

Equations for C_HbO2 and C_MbO2

One mole of Hb can bind with four moles of O₂ and one mole of Mb can bind with one mole of O₂, so the concentrations of Hb- and Mb-bound O₂(C_HbO2 and C_MbO2) can be calculated from the O₂ saturations of Hb and Mb (S_HbO2 and S_MbO2) using

where C_Hbrbc and C_Mbpc are the concentrations of Hb and Mb in the water spaces of RBC and parenchymal cells. With an Hb content of 0.333 g ml⁻¹ RBC (0.15 g ml⁻¹ blood, Hct = 0.45), an RBC water fraction of 0.72 ml ml⁻¹ RBC, and MW_Hb = 64, 458 for the tetramer, then C_Hbrbc = 7.2 mM. Myoglobin concentrations in myocardial cell vary widely among species, but this example model used those in rabbit hearts, an Mb content of 0.0067 g ml⁻¹ cell (0.005 g ml⁻¹ tissue and cell volume fraction in tissue of 0.75 ml ml⁻¹ tissue), a myocyte water fraction of 0.8 ml ml⁻¹ cell, and MW_Mb = 16, 800, gives C_Mbpc = 0.5 mM.

Equations for S_HbO2 and S_MbO2

Of the several mathematical formulas available for describing the relationship of P_O2 to the O₂ saturation of Hb and Mb, S_HbO2 and S_MbO2 (see Popel⁵⁴ for a review), the modified Hill equation for S_HbO2 of Dash and Bassingthwaighte²⁶ allows the effects of intra-RBC pH, P_CO2, 2,3-DPG and temperature as well as the nonlinear O₂–CO₂ interactions to be incorporated into the model for simultaneous O₂ and CO₂ transport and exchange during transients. The expressions for S_HbO2 and S_MbO2 are

where K_HbO2 (with units M⁻¹) is a function of P_O2, P_CO2, and pH, or equivalently ^∘C_rbc,^c C_rbc, and ^hC_rbc, and K_MbO2 (M⁻¹) are given by

and $K_4^{\prime }$ (with units M⁻¹) is given by

The equilibrium constants $K_2^{\prime },K_2^{″},K_3^{\prime },K_3^{″},K_4^{″},K_5^{″}$, and $K_6^{″}$ and the empirical exponents n₀, n₁, and n₂ are as defined in Dash and Bassingthwaighte²⁶ (see Table 1 for values). The subscript S refers to values in standard physiological conditions at 37°C. The concentration of 2,3-DPG is assumed fixed at 4.65 mM. $P_{50}^{\text{Mb}}$ is the value of PO₂ at which Mb is 50% saturated by O₂; $P_{50}^{\text{Mb}}=2.39\text{mmHg}$^59,60; the $P_{50}^{\text{Hb}}$ for human hemoglobin is about 26.8 mmHg.

Michaelis–Menten Kinetics for Oxygen Consumption

The consumption of oxygen is facilitated by the mitochondrial cytochrome-c oxidase reaction and regulated via the availability of substrates (ADP, Pi, NADH, and FADH) to mitochondrial electron transport chain. The oxygen consumption in Eq. (2d) is the product of a “gulosity” ^∘G_pc (with units of a clearance or flow) times the concentration, with ^∘G_pc for MM kinetics being

where V_max (mol min⁻¹ g⁻¹) is the maximal rate of O₂ consumption and K_m is the Michaelis constant. Under standard physiological conditions, V_max is about 2.5 × 10⁻⁶ mol min⁻¹ g⁻¹ and K_m about 7 × 10⁻⁷ M (which corresponds to about 0.5 mmHg in the units of partial pressure as α_O2 = 1.46 × 10⁻⁶ M mmHg⁻¹.^17,31,51 The estimate of K_m about 0.7 μM is a functional one that makes excellent sense relative to myoglobin $P_{50}^{\text{Mb}}$ about 2.39 mmHg. Furthermore, the K_m value from in vitro measurements on isolated mitochondria are generally found to be 0.5–1 Torr P_o2, about 0.7 μM [O₂]. A K_m value of 0.7 μM is also in agreement with the “apparent” K_m value of 0.8 μM estimated by Korzeniewski and Zoladz⁴⁹ for the cytochrome-c oxidase reaction.

Equations for Carbon Dioxide

CO₂ along with other metabolic by-products such as water (H₂O) has a positive net flux from tissue to effluent blood. It undergoes reactions during its transport and exchange process.^26,37 In all regions, CO₂ is transported in both dissolved form and as bicarbonate ions, ${\text{HCO}}_3^{-}$. In RBC, it is also in combination with hemoglobin as carbamino hemoglobin, HbCO₂, which is assumed to be in constant equilibrium with the dissolved CO₂. The partial differential equations for CO₂ in the four regions are

where K_rbc, K_pl, K_isf, and K_pc (unitless) are the catalytic factors for CO₂ hydration reaction due to the presence of enzyme carbonic anhydrase in the four regions of the BTEX unit; K_pl is small, perhaps less than 100, but K_rbc, K_isf, and K_pc are of the order of 10,000; also ${kf}_1=[{\mathrm{H}}_2\mathrm{O}]k{f}_1^{\prime }=0.12{\mathrm{s}}^{-1}$ and ${kb}_1={kb}_1^{\prime }/K_1^{″}=89{\mathrm{s}}^{-1}/5.5\times {10}^{-4}\mathrm{M}=162,000{\mathrm{M}}^{-1}{\mathrm{s}}^{-1}$. These parameters are obtained from Hill et al.,^39–41 Salathe et al.,⁵⁸ and Huang and Hellums.⁴⁴ The last term in each of the equations represents the net rate of transformation of CO₂ to or from bicarbonate. In Eq. (8d), the rate of production of CO₂ inside the parenchymal cells is expressed as the rate of consumption of O₂ times the respiratory quotient, RQ. The concentration of dissolved CO₂ in a region r, ^cC_r, is related to the partial pressure of CO₂ in that region, P_CO2,r , by Henry’s law: ^cC_r = α_CO2 × P_CO2,r , where α_CO2 is the solubility of CO₂ in water. At body temperature (T = 37°C), α_CO2 = 3.27 × 10⁻⁵ M mmHg⁻¹.

Equations for C_HbCO2 and S_HbCO2

As one mole of Hb binds with 4 moles of CO₂, the concentration of Hb-bound CO₂ (C_HbCO2) can be calculated from CO₂ saturation of hemoglobin (S_HbCO2) using the formula:

In present notations, the equation for S_HbCO2 from²⁶ is given by

where K_HbCO2 (a function of P_O2, P_CO2, and pH, or ^∘C_rbc, ^cC_rbc, and ^hC_rbc) is given by

Thus, Eq. (10) for carbon dioxide and Eq. (4) for oxygen can both be calculated directly from the ambient concentrations of the free unbound solutes at each instant in time, critical to obtaining rapid computation. Equation (6) defines $K_4^{\prime }$.

Equations for Bicarbonate Ions

The amount of CO₂ transported as bicarbonate ions ( ${\text{HCO}}_3^{-}$) is governed by the CO₂ hydration reaction. At equilibrium or steady state, the exchange of ${\text{HCO}}_3^{-}$ through the band-3 anion channel is governed by the Gibbs–Donnan electrochemical equilibrium conditions: Eqs. (1a)–(1c). The dynamics of transport and exchange of ${\text{HCO}}_3^{-}$, indicated by the left superscript b, in the four regions are given by

As for the CO₂ equations, the last terms in Eqs. (12a)–(12d) represent the net rate of formation of ${\text{HCO}}_3^{-}$, and are the same as the last terms in Eqs. (8a)–8d.

Equations for Hydrogen Ions and Acid–Base Balance

Because the reactions of O₂ and CO₂ with Hb, Mb, and ${\text{HCO}}_3^{-}$ are influenced greatly by variation in pH,²⁶ the transport and exchange of hydrogen ions (H⁺) must be incorporated into any model for O₂, CO₂, and ${\text{HCO}}_3^{-}$. Variation in pH due to the CO₂ hydration reaction is controlled by the buffering action of bicarbonate ( ${\text{HCO}}_3^{-}$) in both blood and tissue. It is characterized by the buffering capacity parameter, β, which is the amount of base added, or equivalently acid removed, per unit pH change: β = +Δ[base]/ΔpH = −Δ [acid]/ΔpH. Here, the acid is H₂CO₃ (carbonic acid) and the base is ${\text{HCO}}_3^{-}$. As for ${\text{HCO}}_3^{-}$, the exchange of H⁺ through the band-3 proton channels at equilibrium or steady state is governed by the Gibbs–Donnan electrochemical equilibrium conditions: Eqs. (1a)–(1c). The four-region equations are

where β_rbc, β_pl, β_isf, and β_pc (with units M pH⁻¹ unit) are the buffering capacities in the four regions of the BTEX unit; β_pl is small, but β_rbc, β_isf, and β_pc are high; we choose β_pl = 6 mM pH⁻¹ unit, β_rbc = 54 mM pH⁻¹ unit, β_isf = 24 mM pH⁻¹ unit, and β_pc = 45 mM pH⁻¹ unit from Hill et al.^39–41 Salathe et al.,⁵⁸ and Huang and Hellums.⁴⁴

Equations for Capillary Hematocrits and RBC and Plasma Volumes

The capillary hematocrits (Hct_cap) are lower than the large-vessel hematocrits (Hct) by the Fahraeus–Lindquist effect^3,32,43 and consequently the RBC velocity (ν_rbc) is higher than the plasma velocity (ν_pl), and the erythrocytes’ mean transit times through the coronary system is shorter than that of plasma.⁹ The result is a “red cell carriage effect”³³ whereby solutes carried within RBCs travel faster than plasma proteins. In this model, Hct, blood flow (F_bl), capillary volume (V_cap), and the ratio of RBC to plasma velocity (R_vel) are known. The standard values for myocardial blood flow are Hct = 0.45, F_bl = 1 ml min⁻¹g⁻¹, V_cap = 0.07 ml g⁻¹, and R_vel = 1.6. The Hct_cap, RBC and plasma flows (F_rbc and F_pl), their volumes (V_rbc and V_pl), and their velocities (ν_rbc and ν_pl) are calculated through the following equations,⁵¹ the values of which are presented in Table 1 for standard physiological conditions:

The water fractions of various spaces are given by W_pl = 0.94, W_rbc = 0.72, W_isf = 0.92, W_pc = 0.8.51 The anatomical volumes of ISF and parenchymal cells are given by V_isf = 0.2 ml g⁻¹ and V_pc = 0.7 ml g⁻¹, so the water spaces or volumes of distribution of these two regions are about 0.18 and 0.56 ml g^−1.9

Steady-State Extraction and the Consumption Rate of Oxygen

The fractional steady-state extraction of O₂(E_ss) and its consumption rate (MR_O2) can be obtained from total O₂ content in the arterial and venous blood (^∘C_{Tot Art}): and ^∘C_{Tot Ven}):

where

The MR_O2 in Eq. (15) equals the integral of ^∘G_pc^∘C_pc over the axial length L in Eq. (2d).

Computational Methods

To obtain the solutions for concentration profiles, one must prescribe initial and boundary conditions for the governing partial differential equations (2a)–(2d), (8a)– (8d), (12a)–(12d), and (13a)–(13d). For equilibrium/steady-state solutions, the initial conditions are chosen arbitrarily, ^sC_r(x, t) = ^sC_r0(x). The boundary conditions at the entrance to the capillary are defined by Robin or Danckwerts conditions such that the convective inflow is balanced by the diffusive upstream efflux into the inflow:

where v_r is the velocity, F • L/V_r, and D is the diffusion coefficient, as applied to both RBC and plasma regions, and where ^sC_in is the concentration of the species ‘s’ (O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺) in the inflowing arterial blood, where plasma and RBC solutes are assumed to have equilibrated prior to entry into the exchange region at x = 0. These can be computed from standard arterial P_O2, P_CO2, and pH (see Table 1). The concentration of ${\text{HCO}}_3^{-}$ at the arterial end is obtained using the Henderson–Hasselbalch equation and Gibbs–Donnan electrochemical equilibrium condition.²⁶

At the downstream end of the capillary and at both the upstream and downstream ends in all the stagnant extravascular regions, the boundaries are reflecting, meaning that the derivatives are zero so that there is no flux across the boundary:

The model is coded and run in the JSim model simulation environment²³ and is solved by invoking either the LSFEA (Lagrangian Sliding Fluid Element Algorithm,¹³) or TOMS690 Chebychev Polynomial Algorithm¹⁹ or TOMS731 Moving-Grid Algorithm²¹ from our library of numerical PDE solvers. LSFEA is the fastest of these, even though accurate solutions at the inlet require a finer grid or a larger number of axial segments than do the TOMS algorithms.¹³ In contrast, TOMS731 is considerably faster than TOMS690, and the solutions are equally accurate. All the results presented in this paper were obtained using the TOMS731 solver with 51 axial grid points. The time step of integration is adaptively controlled by the solver. The solver computes the concentrations of total O₂ in RBC and parenchymal cells (^oC_rbc + C_HbO2 and ^oC_pc + C_MbO2) and total CO₂ in RBC (^cC_rbc + C_HbCO2) at each axial grid point at each time point. The concentrations of free O₂ in RBC and parenchymal cells (^oC_rbc and ^oC_pc) and free CO₂ in RBC (^cC_rbc) are then obtained iteratively from the coupled Eqs. (4), (5), (9), and (10).

This model’s capabilities extend those of the previous models,^20,27,51,55 and it is computationally fast in spite of its comprehensive nature, so this model is our preferred vehicle for data analysis for oxygen-related studies, for example for ¹⁵O-oxygen studies using positron emission tomography where many regions of interest are analyzed within an organ. The model is available for download and public use from the Physiome website at: http:/nsr.bioeng.washington.edu/jsim/models/webmodel/NSR/Exchange_O2_CO2_HCO3_and_H/

The Default Parameter Set

The model requires many parameters, and while there are many sources for partial information, it has not been possible to find a consistent set of values in the literature. Consequently, a self-consistent set of parameter values were chosen to represent realistic normal physiological conditions; as listed in Table 1 with references. Where the value of a parameter was not available, known relationships to other parameters (correlation) were used to constrain the missing parameter value. When applying the model to data analysis or simply making predictions from the model, the aim is to fix the values for as many parameters as possible, and thereby to minimize the degrees of freedom for data fitting or exploring behavior under varying physiological conditions. Known physicochemical parameters can be fixed. Anatomic parameters can often be fixed or markedly constrained using data coming from other experiments than the one being analyzed; parameters for which there is low sensitivity, e.g., cell volumes and water spaces, fit into this category.⁶⁵ In fitting the model to actual experimental data (e.g., tracer-dilution outflow curves following the injection of tracer-labeled ¹⁵O-oxygen in PET imaging studies) to estimate the physiological parameters of particular interest, oxygen consumption, the RQ, or response times for a change in metabolism, all but a few parameters need to be determined by optimization to fit data.^27,51

With Table 1’s parameter values, the effects of blood flow, hematocrits, arterial P_O2, and oxygen consumption on the concentration profiles can be studied by changing parameter values one at a time, a kind of sensitivity analysis. The results presented later in Figs. 2–6, give a portrayal of the model functions and the model’s utility in exploring physiological interventions.

Axial Profiles at Standard Physiological Conditions

The computed steady-state axial profiles for O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺ in the four regions of the blood–tissue exchange unit using the standard parameter values presented in Table 1 are shown in Fig. 2. The local ^∘G_pc at each x is defined by Eq. (7). Oxygen P_O2 (left upper) decays along the capillary, while the concentrations of metabolic products are all rising, P_CO2 (right upper), bicarbonate [ ${\text{HCO}}_3^{-}$] (left lower), and H⁺ ion (right lower), in all four regions (RBC, plasma, ISF, and parenchymal cells). The arterial inflow concentrations were constant at P_O2 = 100 mmHg, P_CO2 = 40 mmHg, pH_pl = 7.4, pH_rbc = 7.24, ${[{\text{HCO}}_3^{-}]}_{\text{pl}}=24.4\text{mM}$, and ${[{\text{HCO}}_3^{-}]}_{\text{rbc}}=16.8\text{mM}$. In the venous outflow, the computed solutions were P_O2 = 39 mmHg, P_CO2 = 46.5 mmHg, pH_pl = 7.36, pH_rbc = 7.2, ${[{\text{HCO}}_3^{-}]}_{\text{pl}}=26.4\text{mM}$, and ${[{\text{HCO}}_3^{-}]}_{\text{rbc}}=18.2\text{mM}$. These simulated results are in good agreement with the literature and standard physiological observations (e.g., see Huang and Hellums³⁷ and Hellums et al.⁴⁴).

Figure 2 shows steep gradients in O₂, more or less exponentially decaying, while the increasing concentrations in CO₂, ${\text{HCO}}_3^{-}$, and H⁺ show considerably smaller gradients because they are heavily buffered. In contrast to the exponential decay in dissolved free O₂ concentration, the Hb-bound O₂ concentration decays very nearly linearly from a saturation, S_HbO2 = 0.98 at the entrance to S_HbO2 = 0.74 at the exit. This is due to the nonlinear (sigmoid) nature of the oxyhemoglobin dissociation curve and nonlinear–nonlinear interaction of sigmoid–exponential curves. The trend in Hb-bound CO₂ concentration was opposite, though HbCO₂ increased only slightly along the capillary length. The steady-state O₂ extraction at the venous end was 25.4%. The average intracellular O₂ consumption, MR_O2, was under these conditions equal to V_max, 2.5 μ mol min⁻¹g⁻¹, so the intracellular O₂ level was sufficient to saturate the cytochrome-c and to put cellular respiration into the zero-order mode, nearly constant at all positions along the BTEX unit. The [O₂] in RBC and plasma as well as in ISF and parenchymal cells were close to each other, indicating the RBC and plasma as well as ISF and parenchymal cells regions were equilibrated with O₂. However, even with ^∘P S_cap at the rather high value of 200 ml g⁻¹ min⁻¹ there was an [O₂] gradient across the capillary wall, which no experimental evidence has either corroborated or contradicted, suggesting that this permeability is underestimated. In contrast, [CO₂] was in close equilibrium in all four regions of the blood–tissue exchange unit (except near the arterial end), as there was negligible gradient in [CO₂] across the membranes. The [ ${\text{HCO}}_3^{-}$] and [H⁺] were found to vary negligibly in axial direction (probably due to the buffering action of ${\text{HCO}}_3^{-}$) and were in equilibrium in different regions obeying the Gibbs–Donnan electrochemical equilibrium conditions given by Eqs. (1a)–(1c). The pH inside the parenchymal cells, driven by the production of CO₂ in this oversimplified model of metabolism, is too low compared to observed pHs in normal myocardium, simply because there is no representation of H⁺ ion consumption by the processes of oxidative metabolism in the model. This lack reduces all regional pHs compared to those for a fully-developed metabolic model.

Effect of Varying Blood Flow, Hematocrit, and Arterial P_O2 on Axial Profiles

The delivery of O₂ to the tissue is determined mainly by three factors: Capillary blood flow, hematocrit, and arterial P_O2. These also affect the axial profiles of O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺. All these have similar qualitative effects and can be explained from the behavior of one of the factors. The effect of varying capillary blood flow (F_bl) on steady-state axial profiles of O₂, CO₂, and H⁺ in the four regions of the blood–tissue exchange unit are shown in Figs. 3 and 4. Other parameters, including the conditions at the arterial end, are fixed at the standard physiological values presented in Table 1. The steady-state axial profiles of ${\text{HCO}}_3^{-}$ (not shown) can be deduced from the steady-state axial profiles of P_CO2 and pH (using the Henderson–Hasselbalch equation). Figure 3 actually shows the steady-state axial variations of P_O2 and P_CO2 in RBC and parenchymal cells for F_bl = 1.5, 1.0, 0.5 ml min⁻¹ g⁻¹. The profiles of P_O2 and P_CO2 in plasma and ISF are not included in these plots because these profiles almost match with those in RBC and parenchymal cells. Figure 4 shows the steady-state axial profile of pH in RBC, plasma, ISF, and parenchymal cells for the same three values of F_bl.

It is seen that as the capillary blood flow (perfusion or O₂ supply) decreases with the O₂ consumption, while V_max (O₂ demand) remains constant, the concentration of O₂ in the blood–tissue exchange unit decreases. However, the P_O2 does not go to zero anywhere. With arterial P_O2 = 100 mmHg, V_max = 2.5 μ mol min⁻¹ g⁻¹, and F_bl = 1.5, 1.0, 0.5 ml min⁻¹ g⁻¹, we see that the venous P_O2 is 46, 39, 27 mmHg and O₂ extraction is 17.3, 25.4, 48.1%. A decrease in blood flow also results in an increase in the concentrations of CO₂, H⁺, and ${\text{HCO}}_3^{-}$ in the system. Thus, a severe blood flow decrease (ischemia) can result in acute hypoxia and acidosis (lower pH) inside the parenchymal cells and affect the normal cellular function.

Effect of Varying Oxygen Consumption on Axial Profiles

The effect of varying O₂ consumption V_max on steady-state axial profiles O₂, CO₂, and H⁺ in the four regions of the blood–tissue exchange unit is shown in Figs. 5 and 6. Other parameters, including the conditions at the arterial end, are fixed at the standard physiological values presented in Table 1. Figure 5 shows the steady-state axial variations of P_O2 and P_CO2 in RBC and parenchymal cells for V_max = 2.5, 5.0, 7.5 μmol min⁻¹ g⁻¹. Figure 6 shows the steady-state axial variation of pH in RBC, plasma, ISF, and parenchymal cells for the same three values of V_max. These results are consistent with the results in Figs. 3 and 4, i.e., the increasing O₂ consumption V_max (O₂ demand) has similar qualitative effects on the species concentration profiles as does decreasing capillary blood flow (perfusion or O₂ supply). Both cause decreases in the regional O₂ concentration and increases in the regional concentrations of CO₂, H⁺, and ${\text{HCO}}_3^{-}$ in the system. Quantitatively the differences between the regional concentration profiles are larger with varying V_max than with varying F_bl, using the same relative increments, comparing Fig. 6 with Fig. 4. With arterial P_O2 = 100 mmHg, F_bl = 1 ml min⁻¹ g⁻¹, and V_max = 2.5, 5.0, 7.5 μmol min⁻¹ g⁻¹, the venous P_O2 is seen to vary from 39 to 26 to 20 mmHg and O₂ extraction is seen to vary from 25.4 to 50.2 to 66.8%. With V_max = 5.0 μmol min⁻¹ g⁻¹ the downstream end of the parenchymal cells has a P_O2 of 10 mmHg and pH of 7.03. With V_max = 7.5 μmol min⁻¹ g⁻¹, the downstream half of the parenchymal cells has a P_O2 < 5 mmHg and pH of 7.01. Therefore, when there is an increasing energy demand and O₂ consumption (associated with exercise), unless the capillary blood flow and O₂ supply increase to meet the energy demand, the normal cellular function will be severely affected.

It may be noted here that in acute situations with very low blood flow (F_bl) or very high oxygen consumption (V_max) the blood P_O2 can become very low (below 20 mmHg) in the downstream end of the blood–tissue exchange region, and the HbO₂ saturation will fall below 30%. Our model for HbO₂ saturation²⁶ begins to lose accuracy below S_HbO = 30%, so the computed solutions are less accurate in this regime of low P_O2.

Potential Applications

This model is the start of, or a part of, whole-body metabolic modeling for systems physiology. While it encompasses only the simultaneous transport and exchanges of O₂, CO₂, ${\text{HCO}}_3^{-}$, and H⁺ in the microcirculation, it links these processes to the metabolism of O₂ and therefore to the whole realm of intracellular metabolism. Augmenting such a model by adding biochemical systems for intermediary metabolism, oxidative phosphorylation, NAD/NADH reactions, the tricarboxylic acid cycle, and hydrogen balance will improve the characterization of spatial profiles and representation of kinetics during transients. Through whole-body representation of multiple organs, such a model can be expanded into a model for glucose metabolism, for example.

The model itself is worthwhile as a basic vehicle for data analysis for whole-organ studies by multiple indicator dilution techniques, by positron emission tomographic or NMR image sequences. For application to BOLD (blood oxygen level dependent) NMR, one has only to pull out the signal for Hb rather than HbO [where HbO really represents Hb(O₂)₄]. As a component of a model further developed to include mitochondrial oxygen and carbon dioxide kinetics, and linked into a whole-organ model by adding the arterial and venous transport systems, one will be able to evaluate the basis for the delays in responses to changes in oxygen utilization with changes in cardiac contractile performance. This has been a puzzle: “Is it mitochondrial delay or is it merely a transport time?”

What is modeled here has been well known for years, but the complexity of the problem and the slowness of computation inhibited developing full sets of equations for practical applications. The model is free for general use and available for download. Investigators anywhere can explore the model behavior with or without delving into the code by going to the website for JSim.²³ This fits with current trends toward open-source software and widespread collaboration in research. Models including the O₂–CO₂ interactions were presented by Huang and Hellums⁴⁴ and Ye et al.,⁶⁷ but for different reasons neither of their approaches lends itself to the quantitative analysis of experimental data from transients in venous P_O2 or tracer ¹⁵O-oxygen in positron emission tomography (PET) imaging studies. The model of⁴⁴ accounts for radial diffusion, but is computationally slow, making it impractical for routine fitting of experimental time courses from PET imaging studies by automated iterative parameter adjustment optimization. The model of Ye et al.⁶⁷ was compartmental, thereby failing to account for axial gradients along the capillary, and therefore lacking the critical narrow impulse response required to fit high-resolution data of the sort analyzed by Deussen and Bassingthwaighte²⁷ and Li et al.⁵¹ The current model solves the first of these problems by an approximation (for radial diffusion) and by using highly efficient computational algorithms for Hb/HbO and for the blood–tissue exchange, and the second by using very efficient partial rather than ordinary differential equations.

This general model can readily be applied to arterioles and venules and to other organs (e.g., the lungs) for gaining insight into the processes of O₂ and CO₂ transport and exchange as well as acid–base regulation along the length of the microcirculatory exchange regions from arterioles to venules.

Improvements Over Prior Models

The current model improves our previous blood–tissue exchange models of O₂ transport and metabolism^20,27,51 for it uses the analytically invertible expressions of Dash and Bassingthwaighte²⁶ for oxyhemoglobin (HbO₂) and carbamino hemoglobin (HbCO₂) dissociation curves which explicitly account for the effects of P_O2, P_CO2, pH, and nonlinear O₂–CO₂ interactions in computing P_O2 and P_CO2 from the total O₂ and CO₂ concentrations (^∘C_rbc + C_HbO2 and ^cC_rbc + C_HbCO2) in RBC at each axial and time grid point. This reduces the computation time by about a factor of 20 compared to solving the Adair-type equations by iteration.²⁶

Computational efficiency is a major issue in analyzing the data from PET imaging studies of tracer-labeled ¹⁵O-oxygen because of the broad heterogeneity of flows and metabolism in the heart (e.g., see King et al.,⁴⁷ Caldwell et al.,²⁴ reviews by Bassingthwaighte et al.¹⁴ The steady-state spatial distribution of regional blood flows within organs is always broad, having coefficients of variation of 25% or more, depending on the tissue under study and the ambient humoral, nervous, and metabolic conditions.^9,47 To allow an approximate description of spatial heterogeneity at the level of spatial resolution of about 0.5–1% of the heart mass (the resolution used for PET analysis) and still account for the heterogeneity within the observed region-of-interest, our capillary–tissue exchange model can be incorporated into a multicapillary configuration in which the blood–tissue exchange regions are composed of a set of parallel, independent capillary–tissue exchange units such as those described.^27,48

The O₂ consumption in such small tissue regions can be measured directly using the tracer kinetics of ¹⁵O-oxygen in PET imaging⁵¹ and in isolated perfused rabbit hearts.²⁷ For standard PET imaging data analysis, one uses 30–100 regions of interest within the LV myocardium, so efficiency is very important. To extend our current nontracer model for the analysis of tracer ¹⁵O-oxygen from PET imaging studies, it needs to be integrated with an identical model for tracer ¹⁵O-oxygen and its product, tracer ¹⁵O-water, which is formed in the mitochondria. The steady-state solutions of the current model will provide the volumes of distribution (virtual volumes) for the tracer ¹⁵O-oxygen that accounts for the Hb binding space and water space. The current model differs from that of Li et al.⁵¹ in the sense that the virtual volumes will be functions of P_O2, P_CO2, and pH, and therefore will vary under varying physiological conditions.

Assessment of the Model’s Features

As reviewed in detail by Hellums et al.,³⁷ there are a number of issues to be addressed in any integrated model that contains simplifications or approximations of the real geometry and physiology. These issues include (1) the exchanges across the arterioles and venules walls and the complexity of diffusional interactions among the capillaries, arterioles, and venules, (2) the radial gradients in blood, particularly in arterioles and venules, and in tissues, (3) the increase in intracapillary diffusional resistance due to the spacings between RBC in the capillary, (4) the existence of a continuous disequilibrium between O₂ and Hb and between O₂ and Mb during steady states involving diffusion and consumption (and likewise for CO₂ with Hb), (5) the anion exchanges ( ${\text{HCO}}_3^{-}$ and Cl⁻) across the membranes, and (6) at the level of molecular interactions, the release of approximately 0.7 H⁺ and binding of 60 H₂O molecules²⁵ when O₂ binds to Hb.

The components of Issue 1, arteriolar and venular exchanges, can be handled approximately, using three or more BTEX units in series, as was done.²⁷ The first serves as the arteriole, the second as the capillary–tissue region and the third as the venule. The parameters for the volumes and permeabilities of the arteriolar and venular sections are adjusted to match reality. The resistances between blood and tissue are the sum of those of the vessel walls (with a minor degree of O₂ consumption by the smooth muscle and other cells) and of the loose aureolar connective tissue surrounding the vessels and allowing their expansion and contraction. Further, the arrangement in many tissues including the heart where each arteriole is accompanied by two venules (venae comitantes⁶) means that some oxygen diffuses from arterioles to interstitium to venules without passing by the muscle cells. While this is most obvious for heat exchange,⁹ it is a major factor in shaping the impulse responses for tracer water¹⁶ and oxygen¹⁷; this can be accounted for, again approximately but rather well, by increasing the D_isf, the dispersion coefficient in the interstitial fluid space, to allow graded fluxes between the upstream and downstream ends of the capillary–tissue exchange unit, which accounts approximately for the asynchrony or offsettedness of capillary entrances and exits. An improved treatment would include the arterio-venous exchanges directly in addition to the dispersion.

Compensating for the omission of radial gradients (Issue 2) is important. O₂ and CO₂ gradients will be poorly represented in arterioles and venules. Radial diffusion in 20–40-micron diameter vessels takes longer than in capillaries where RBCs rub along the capillary wall; decreasing P S_cap can give a good approximation to this enlarged resistance when the resistance and fluxes are calculated from a more detailed model such as that of Huang and Hellums.⁴⁴ Their model tested very well against the data on 1.5 mm diameter tubes from Dorson and Voorhees²⁸ on O₂ and CO₂ exchange. Fortunately the fluxes across walls of arterioles and venules are small,^44,67 so finding an exact correction is not critical. As for intra-tissue diffusion, with half-intercapillary distances of about 8 microns in the heart,⁶ the diffusion relaxation times are short enough (τ = R²/2D or about 30 ms maximum) that the radial gradients are almost dissipated within one time step for the numerical solution, allowing us to consider them negligible within the shorter distances across each of the radial regions.

Issue 3 relates to the hematocrit-dependency of the effective blood–tissue exchange, defined by Hellums³⁶ in a nicely principled analysis: Transfer of O₂ is faster when plasmatic gaps between RBC are smaller. Others^{29,30,34,37,44} discussed that even in steady-state modeling there was an influence of RBC spacing in the capillary on the efficiency of O₂ transfer from blood to the tissue when the hematocrit was lower than 20%. This is not directly addressed in the current model, but the reduction in oxygen transfer due to the decreased P_O2 within the plasmatic gaps can be well approximated, as far as delivery to the tissue is concerned, by an appropriate reduction in P S_cap.

Slow release of O₂ from Hb (Issue 4) retards exchange into peripheral tissues. Slow dissociation reduces the effective free O₂ concentration at the RBC membrane and reduces transmembrane flux,³⁷ exactly analogous to the situation for albumin-facilitated diffusion of fatty acid in plasma when there is endothelial uptake.⁴ Likewise, slow association retards RBC uptake of O₂ in the lung, contributing to any alveolar–arterial O₂ concentration difference. Knowing the rate constants and the dimensions involved allows translation into an equivalent increase in the membrane resistance. While the rate constants for CO₂ binding with Hb are known, those for O₂ binding are concentration-dependent and treated empirically, as described by Hellums et al.³⁷, their Eq. (7), where the “on” and “off” rates are defined as a function of HbO₂ saturation to account for the cooperativity.

Slow anion exchange (Issue 5) may be important at the RBC membrane because a delay in ${\text{HCO}}_3^{-}$ and Cl⁻ exchange influences the pH inside the RBC. This would slow the pH change along the length of the exchange region, retarding the rise of pH and carbamino formation in the capillary blood (Fig. 2) and slowing CO₂ loss in the lung. This is not accounted for in our present model.

Finally, Issue 6: There is the release of about 0.7 mole H⁺ per mole of O₂ binding to Hb, and 1.5–2 moles H⁺ per mole of CO₂ binding to Hb,⁴² that are not included in the present model, but could be added stoichiometrically to the equation for H⁺ fluxes. In going from a P_O2 of 40–100 mmHg, 1 mole of Hb goes from 75 to 98% saturated, a gain of 23%. With [Hb] = 7.2 mM (Table 1), the O₂ binding would displace about 0.7 times 0.23 times 7.2 mM, a change in [H⁺] of 1.16 mM, over 3 pH units. Since CO₂ binds less readily to HbO₂ than to Hb, and since in the lung the CO₂ is coming off when O₂ is going on, and vice versa in the systemic tissue, the net H⁺ exchange is always smaller than 1.16 mM. Because the amount of CO₂ bound to Hb is much less than O₂, there is still a net H⁺ balance dominated by the O₂ binding or release.

While it is important to address the preceding six issues in any model of O₂ and CO₂ exchange, the present model is comprehensive enough that the model parameters can be varied to fit many different physiological and pathological conditions. For example, one can look at the distributions of O₂, CO₂, ${\text{HCO}}_3^{-}$ , and H⁺ and the regulation of acid–base balance in tissue cells in normoxia, hypoxic hypoxia (low arterial P_O2), and anemic hypoxia (low hematocrits) conditions, and with varying arterial plasma or RBC pH. The results shown in this paper (Figs. 3–6) are the simulated concentration profiles with decreasing blood flow (F_bl) and increasing O₂ consumption (V_max). These cases correspond to the ischemic or low-perfusion (low F_bl) and exercise or higher energy demand (high V_max) conditions. The transient behavior of these profiles at the onset of exercise (step increase in V_max) from steady state can be computed and analyzed in understanding the dynamics of O₂–CO₂ exchange and acid–base balance.

The respiratory quotient RQ (the ratio of CO₂ production to O₂ consumption) was fixed at an average value of 0.8 (protein metabolism). While RQ can vary from 0.7 (fat metabolism) to 1.0 (carbohydrate metabolism), a model for a regulated change in CO₂ production would require the model to include the major metabolic pathways (glycolysis, glycogenolysis, and TCA cycle) and oxidative phosphorylation^15,49 along with the present model to complete the link between O₂ consumption and CO₂ and H⁺ production in the parenchymal cells.

There are 16 axial dispersion coefficients (^sD_r; 4 species × 4 regions) and 12 permeability surface area products (^sP S_m; 4 species × 3 membranes) in the current model. There is neither enough information nor a real necessity to set them all at distinct values, because the gas exchanges between regions are so rapid. As discussed above and in the “Mathematical Formulation of the Model” section, the Ds were assigned higher values than the pure molecular diffusion coefficients to account for the heterogeneities of spacing and capillary entrance regions. The radial diffusion is a separate issue. While the true membrane PSs are high, and must be for O₂ whose exchange is almost flow-limited in the microcirculation,⁵¹ they are reduced somewhat to account for the retarding influences of diffusion as discussed above. The PSs for O₂ and CO₂ are set at the values used for O₂ by Li et al.⁵¹ in their studies of PET data analysis. The PSs for ${\text{HCO}}_3^{-}$ and H⁺ are set one order of magnitude lower than those for O₂ and CO₂, because their molecular polarity reduces transmembrane transport compared to nonpolar solutes.

In conclusion, the present model describes physiologically realistic features of O₂ transport, exchange, and metabolism, including the simultaneous CO₂ transport and exchange and acid–base balance. It is rather general and is widely applicable to measuring O₂ consumption in the heart and other organs. The model parameters can be varied to simulate different physiological situations. The model is computationally efficient, so the tracer transients can be computed quickly, favoring this model as a preferred vehicle for routine data analysis (e.g., in PET imaging studies).