Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2868387

Formats

Article sections

Authors

Related links

Comput Methods Programs Biomed. Author manuscript; available in PMC 2010 May 12.

Published in final edited form as:

PMCID: PMC2868387

NIHMSID: NIHMS197430

James B. Bassingthwaighte: ude.notgnihsaw.gneoib@bb

The publisher's final edited version of this article is available at Comput Methods Programs Biomed

See other articles in PMC that cite the published article.

A computational model of oxygen transport from red blood cells to mitochondria with subsequent reaction to water is presented. This computational model consists of a five region convection-diffusion-reaction mathematical model which is solved using a standard numerical time-split method. The unique feature of this mathematical model is the treatment of the red blood cells and the plasma as two separate flows. The numerical method is second order accurate overall. This computational model is useful for analyzing residue data from positron emission tomography or data from multiple indicator dilution curves.

Many models for oxygen transport from blood into tissue have been investigated [1–7], but they are mainly concerned with steady state and do not allow for transient solutions. Also, the models are over-simplified and neglect certain critical aspects of the oxygen transport system. The models of Federspiel and Sarelius [3], Federspiel and Popel [4] and of Groebe and Thews [7] take into account that hemoglobin comes in discrete packages in the red cell, a factor which Hellums [8] identified in his 1977 paper as being critical at low hematocrits. However, even these studies like those before them give a first approximation only, since they are concerned with steady state situations. The study of steady state solutions does allow one to gain much understanding of complex situations, as reviewed by Popel [9], but there is a practical need for transient solutions for analyzing experimental data.

In contrast to those, Rose and Goresky [10] developed a model which assumes that hemoglobin is homogeneously distributed in the capillary. Their model was used for describing the transient change in ^{18}O-oxygen by analysis of coronary venous dilution curves following an arterial bolus injection of ^{18}O-labeled oxygen molecules. They concluded that the permeability for oxygen diffusion from capillary into tissue was low relative to that for diffusion across slices of tissue. However, it is possible that the inhomogeneous distribution of the oxygen carrier hemoglobin leads to an extra impediment to oxygen transport [11]. The impediment to oxygen transport does depend on the capillary hemocrit, which will determine the distance between red blood cells. Furthermore, the speed of dissociation of oxygen from the hemoglobin molecule, intraerythrocyte diffusion, membrane permeation, permeation through the endothelial cells and diffusion through the interstitial fluid may all together play a very significant role under normal circumstances.

A similar model by Deussen and Bassingthwaighte [12] advances the field by accounting for the tracer water produced by metabolism of the tracer oxygen. This model, like that of Rose and Goresky, uses linear conductance parameters for the membrane permeability-surface area products and intracellular transformation or clearance. Like theirs, it considers the blood as a “hemoglobin solution,” so that all of the sequential resistance's beginning with release from the binding site on hemoglobin, are lumped together as a single capillary wall conductance or permeability-surface area product, PS_{c}. It is a multicapillary model, allowing accounting for intraorgan flow heterogeneity. The Deussen–Bassingthwaighte model in addition allows axial diffusion parallel to the flow path. Its most important development is not the diffusion, but the transformation of ^{15}O-oxygen to ^{15}O-water, making it a dual species model; it is therefore useful for analyzing ^{15}O-water residue or ^{15}O-water or ^{3}H-water outflow curves, ^{15}O-oxygen, ^{18}O-water, ^{17}O-oxygen, or ^{11}C-carbon monoxide and other solutes that may or may not undergo a single transformation. This model has been applied to studies on the positron emitting tracers mentioned, but has a particularly strong potential for providing estimates of regional oxygen consumption from a sequence of positron emission tomographic (PET) images. The first role of the model proposed here is to extend this capability.

There is a need to calculate local flow and oxygen consumption in an organ from residue curves of ^{15}O-oxygen obtained with PET scans. Compartmental models have been used to analyze PET data on tracer oxygen residue curves for the brain [13,14]. Since oxygen is metabolized to water, and the oxygen label in water and molecular oxygen cannot be discerned, the washout of water must be modeled as well. This new axially distributed model accounts for oxygen transport and metabolism using linearized approximations for hemoglobin and myoglobin binding, but considers radial diffusion within each region to be instantaneous. It accounts for a separation between erythrocytes and plasma and for the more rapid velocity of erythrocytes than plasma, a feature whose importance can now be evaluated. The lumped mixing chamber tracer-kinetic models commonly used for PET detection of labeled water residue curves suffer a severe limitation under certain circumstances because of the compartmental assumption [15], and this distributed model being more anatomically realistic avoids this pitfall.

Thus, the second purpose is to develop a physiologically and anatomically reasonable model, with the inflow-to-capillary discontinuity in concentrations implicit in the compartmental modeling, and to capture its essence in a robust, accurately computable model.

We idealize the geometry of a capillary-tissue exchange unit and model the flowing regions and the tissue as axisymmetric cylindrical elements as shown in Fig. 1. The red blood cells are idealized as a separate homogeneous flowing region contained within the plasma region. Since we assume axisymmetric, we only have to solve five one-dimensional partial differential equations.

The five region red blood cell model (rbc model) consists of the following set of partial differential equations. The notation is consistent with Bassingthwaighte et al. [16]. The concentrations, *C*, for O_{2} are given by the solutions to

$$\begin{array}{l}\frac{\partial {C}_{\text{rbc}}}{\partial t}=-{u}_{\text{rbc}}\frac{\partial {C}_{\text{rbc}}}{\partial x}-\frac{{\text{PS}}_{\text{rbc}}}{V{\text{'}}_{\text{rbc}}}\left({c}_{\text{rbc}}-{C}_{\text{p}}\right)-\frac{{G}_{\text{rbc}}}{V{\text{'}}_{\text{rbc}}}{c}_{\text{rbc}}+{D}_{\text{rbc}}\frac{{\partial}^{2}{C}_{\text{rbc}}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{p}}}{\partial t}=-{u}_{\text{p}}\frac{\partial {C}_{\text{p}}}{\partial x}-\frac{{\text{PS}}_{\text{rbc}}}{{V}_{\text{p}}}\left({C}_{\text{p}}-{C}_{\text{rbc}}\right)-\frac{{\text{PS}}_{\text{g}}}{{V}_{\text{p}}}({C}_{\text{p}}-{C}_{\text{isf}})-\frac{{\text{PS}}_{\text{ec1}}}{{V}_{\text{p}}}({C}_{\text{p}}-{C}_{\text{ec}})-\frac{{G}_{\text{p}}}{{V}_{\text{p}}}{C}_{\text{p}}+{D}_{\text{p}}\frac{{\partial}^{2}{C}_{\text{p}}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{ec}}}{\partial t}=-\frac{{\text{PS}}_{\text{ecl}}}{V{\text{'}}_{\text{ec}}}({C}_{\text{ec}}-{C}_{\text{p}})-\frac{{\text{PS}}_{\text{eca}}}{{V}_{\text{ec}}^{\text{'}}}({C}_{\text{ec}}-{C}_{\text{isf}})-\frac{{G}_{\text{ec}}}{{V}_{\text{ec}}^{\text{'}}}{C}_{\text{ec}}+{D}_{\text{ec}}\frac{{\partial}^{2}{C}_{\text{ec}}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{isf}}}{\partial t}=-\frac{{\text{PS}}_{\text{g}}}{{V}_{\text{isf}}^{\text{'}}}({C}_{\text{isf}}-{C}_{\text{p}})-\frac{{\text{PS}}_{\text{eca}}}{{V}_{\text{isf}}^{\text{'}}}({C}_{\text{isf}}-{C}_{\text{ec}})-\frac{{\text{PS}}_{\text{pc}}}{{V}_{\text{isf}}^{\text{'}}}({C}_{\text{isf}}-{C}_{\text{pc}})-\frac{{G}_{\text{isf}}}{{V}_{\text{isf}}^{\text{'}}}{C}_{\text{isf}}+{D}_{\text{isf}}\frac{{\partial}^{2}{C}_{\text{isf}}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{pc}}}{\partial t}=-\frac{{\text{PS}}_{\text{pc}}}{{V}_{\text{pc}}^{\text{'}}}({C}_{\text{pc}}-{C}_{\text{isf}})-\frac{{G}_{\text{pc}}}{{V}_{\text{pc}}^{\text{'}}}{C}_{\text{pc}}+{D}_{\text{pc}}\frac{{\partial}^{2}{C}_{\text{pc}}}{\partial {x}^{2}},\end{array}$$

and the concentrations for H_{2}O are given by the solutions to an analogous set of equations but in which there is a positive production term *G _{i}C_{i}* from the flux of oxygen tracer being transformed to water tracer.

$$\begin{array}{l}\frac{\partial {C}_{\text{rbc}}^{\ast}}{\partial t}=-{u}_{\text{rbc}}\frac{\partial {C}_{\text{rbc}}^{\ast}}{\partial x}-\frac{{\text{PS}}_{\text{rbc}}^{\ast}}{{V}_{\text{rbc}}^{\ast}\text{'}}({C}_{\text{rbc}}^{\ast}-{C}_{\text{p}}^{\ast})+\frac{{G}_{\text{rbc}}}{{V}_{\text{rbc}}^{\ast}\text{'}}{C}_{\text{rbc}}+{D}_{\text{rbc}}^{\ast}\frac{{\partial}^{2}{C}_{\text{rbc}}^{\ast}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{p}}^{\ast}}{\partial t}=-{u}_{\text{p}}\frac{\partial {C}_{\text{p}}^{\ast}}{\partial x}-\frac{{\text{PS}}_{\text{rbc}}^{\ast}}{{V}_{\text{p}}}({C}_{\text{p}}^{\ast}-{C}_{\text{rbc}}^{\ast})-\frac{{\text{PS}}_{\text{g}}^{\ast}}{{V}_{\text{p}}}({C}_{\text{p}}^{\ast}-{C}_{\text{isf}}^{\ast})-\frac{{\text{PS}}_{\text{ecl}}^{\ast}}{{V}_{\text{p}}}({C}_{\text{p}}^{\ast}-{C}_{\text{ec}}^{\ast})+\frac{{G}_{\text{p}}}{{V}_{\text{p}}}{C}_{\text{p}}+{D}_{\text{p}}^{\ast}\frac{{\partial}^{2}{C}_{\text{p}}^{\ast}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{ec}}^{\ast}}{\partial t}=-\frac{{\text{PS}}_{\text{ecl}}^{\ast}}{{V}_{\text{ec}}^{\ast}\text{'}}({C}_{\text{ce}}^{\ast}-{C}_{p}^{\ast})-\frac{{\text{PS}}_{\text{eca}}^{\ast}}{{V}_{\text{ec}}^{\ast}\text{'}}({C}_{\text{ec}}^{\ast}-{C}_{\text{isf}}^{\ast})+\frac{{G}_{\text{ec}}}{{V}_{\text{ec}}^{\ast}\text{'}}{C}_{\text{ec}}+{D}_{\text{ec}}^{\ast}\frac{{\partial}^{2}{C}_{\text{ec}}^{\ast}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{isf}}^{\ast}}{\partial t}=-\frac{{\text{PS}}_{\text{g}}^{\ast}}{{V}_{\text{isf}}^{\ast}\text{'}}({C}_{\text{isf}}^{\ast}-{C}_{\text{p}}^{\ast})-\frac{{\text{PS}}_{\text{eca}}^{\ast}}{{V}_{\text{isf}}^{\ast}\text{'}}({C}_{\text{isf}}^{\ast}-{C}_{\text{ec}}^{\ast})-\frac{{\text{PS}}_{\text{pc}}^{\ast}}{{V}_{\text{isf}}^{\ast}\text{'}}({C}_{\text{isf}}^{\ast}-{C}_{\text{pc}}^{\ast})+\frac{{G}_{\text{isf}}}{{V}_{\text{isf}}^{\ast}\text{'}}{C}_{\text{isf}}+{D}_{\text{isf}}^{\ast}\frac{{\partial}^{2}{C}_{\text{isf}}^{\ast}}{\partial {x}^{2}},\\ \frac{\partial {C}_{\text{pc}}^{\ast}}{\partial t}=-\frac{{\text{PS}}_{\text{pc}}^{\ast}}{{V}_{\text{pc}}^{\ast}\text{'}}({C}_{\text{pc}}^{\ast}-{C}_{\text{isf}}^{\ast})+\frac{{G}_{\text{pc}}}{{V}_{\text{pc}}^{\ast}\text{'}}{C}_{\text{pc}}+{D}_{\text{pc}}^{\ast}\frac{{\partial}^{2}{C}_{\text{pc}}^{\ast}}{\partial {x}^{2}},\end{array}$$

where *C* refers to O_{2} and *C** refers to H_{2}O. The subscripts are rbc for red blood cells, p for plasma, ec for endothelial, isf for interstitial fluid, and pc for parenchymal cells.

For notational convenience, we define the following constants for the ratios of regional volumes of distribution to the plasma space:

$$\begin{array}{llll}{\mathit{\gamma}}_{\text{rbc}}=\frac{{V}_{\text{rbc}}^{\prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{ec}}=\frac{{V}_{\text{ec}}^{\prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{isf}}=\frac{{V}_{\text{isf}}^{\prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{pc}}=\frac{{V}_{\text{pc}}^{\prime}}{{V}_{\text{p}}},\\ {\mathit{\gamma}}_{\text{rbc}}^{\ast}=\frac{{V}_{\text{rbc}}^{\ast \prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{ec}}^{\ast}=\frac{{V}_{\text{ec}}^{\ast \prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{isf}}^{\ast}=\frac{{V}_{\text{isf}}^{\ast \prime}}{{V}_{\text{p}}},& {\mathit{\gamma}}_{\text{pc}}^{\ast}=\frac{{V}_{\text{pc}}^{\ast \prime}}{{V}_{\text{p}}}.& & & \end{array}$$

Also, the following rate constants are defined again relative to the plasma space:

$$\begin{array}{llll}{k}_{1}=\frac{{G}_{\text{rbc}}}{{V}_{\text{p}}},& {k}_{12}=\frac{{\text{PS}}_{\text{rbc}}}{{V}_{\text{p}}},& {k}_{12}^{\ast}=\frac{{\text{PS}}_{\text{rbc}}^{\ast}}{{V}_{\text{p}}},& \\ {k}_{2}=\frac{{G}_{\text{p}}}{{V}_{\text{p}}},& {k}_{23}=\frac{{\text{PS}}_{\text{ecl}}}{{V}_{\text{p}}},& {k}_{24}=\frac{{\text{PS}}_{\text{g}}}{{V}_{\text{p}}},& \\ {k}_{23}^{\ast}=\frac{{\text{PS}}_{\text{ec1}}^{\ast}}{{V}_{\text{p}}},& {k}_{24}^{\ast}=\frac{{\text{PS}}_{\text{g}}^{\ast}}{{V}_{\text{p}}}& & \\ {k}_{3}=\frac{{G}_{\text{ec}}}{{V}_{\text{p}}},& {k}_{34}=\frac{{\text{PS}}_{\text{eca}}}{{V}_{\text{p}}},& {k}_{34}^{\ast}=\frac{{\text{PS}}_{\text{eca}}^{\ast}}{{V}_{\text{p}}},& \\ {k}_{4}=\frac{{G}_{\text{isf}}}{{V}_{\text{p}}},& {k}_{45}=\frac{{\text{PS}}_{\text{pc}}}{{V}_{\text{p}}},& {k}_{45}^{\ast}=\frac{{\text{PS}}_{\text{pc}}^{\ast}}{{V}_{\text{p}}},& {k}_{5}=\frac{{G}_{\text{pc}}}{{V}_{\text{p}}},\end{array}$$

where the *k _{i}'s* refer to O

The known inputs, constant over time to the model are: whole blood flow *F*_{blood} in units of ml g^{−l} min^{−1}, length of capillary *L* in units of cm, the total volume *V*_{cap} of the capillary in units of ml, the red blood cell virtual volume ratio
${\mathit{\theta}}_{\text{rcb}}={V}_{\text{rcb}}^{\prime}/{V}_{\text{rcb}}$, various permeability-surface area products PS in units of ml g^{−l} min^{−1}, reaction rates *G* in units of ml g^{−l} min^{−1}, flow rates *F* in units of ml g^{−l} min^{−1}, and diffusion coefficients *D* in units of cm^{2} s^{−1}. The small vessel hematocrit is lower than the larger vessel hematocrit (and the rbc velocities higher than plasma velocities in small vessels), so we account for this at the capillary level only, using the hematocrit ratio Hct_{ratio} (Hct_{ratio} = Hct_{cap}/Hct_{LV}), where the large vessel hematocrit is Hct_{LV}. Normally, the hematocrit is physiologically less than 1 and rbc velocity greater than 1, so that the time for exchange between rbc and tissue is reduced compared to that for solutes carried only in the plasma. Other parameters that are needed are calculated from these inputs. These are:

$$\begin{array}{l}{V}_{\text{rbc}}={V}_{\text{cap}}{\text{Hct}}_{\text{LV}}{\text{Hct}}_{\text{ratio}},\\ {V}_{\text{p}}={V}_{\text{cap}}-{V}_{\text{rbc}},\\ {\text{Hct}}_{\text{cap}}=\frac{{V}_{\text{rbc}}}{{V}_{\text{cap}}},\\ {\text{Hct}}_{\text{LV}}=\frac{{V}_{\text{rbc}-\text{LV}}}{{V}_{\text{p}-\text{LV}}+{V}_{\text{rcb}-\text{LV}}},\\ {V}_{\text{rbc}}^{\prime}={\mathit{\theta}}_{\text{rbc}-{\text{O}}_{2}}{V}_{\text{rbc}},\\ {V}_{\text{rbc}}^{\ast \prime}={\mathit{\theta}}_{\text{rbc}-{\text{H}}_{2}\text{O}}{V}_{\text{rbc}},\\ {F}_{\text{rbc}}={F}_{\text{blood}}{\text{Hct}}_{\text{LV}},\\ {F}_{\text{p}}={F}_{\text{blood}}(1-{\text{Hct}}_{\text{LV}}),\\ {u}_{\text{p}}=\frac{{F}_{\text{p}}L}{{V}_{\text{p}}},\\ {u}_{\text{rbc}}=\frac{{F}_{\text{rbc}}L}{{V}_{\text{rbc}}},\\ {\mathit{\theta}}_{\text{rbc}-{O}_{2}}=\frac{{V}_{\text{rbc}}^{\prime}}{{V}_{\text{rbc}}}=[\text{Hgb}]f(p{\text{O}}_{\text{2}}),\\ {\mathit{\theta}}_{\text{rbc}-{\text{H}}_{\text{2}}\text{O}}=\frac{{V}_{\text{rbc}}^{\ast}\prime}{{V}_{\text{rbc}}},\end{array}$$

where [Hgb] is the concentration of hemoglobin, and *f*(*p*O_{2}) is the functional relationship between O_{2} saturation and O_{2} partial pressure. The functional relationship *f*(*p*O_{2}) for the hemoglobin-containing erythrocytes and the myoglobin-containing myocytes or parenchymal cells are described in detail by Deussen and Bassingthwaighte [12]. They give methods for calculating the diffusion terms. The rbc and plasma velocities are *u*_{rbc} and *u*_{p}.

In matrix notation, the rbc model is

$$\begin{array}{cc}{c}_{t}=\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}-\mathbf{\text{U}}{c}_{x}+{\mathbf{\text{K}}}_{c},& 0\le x\le L.\end{array}$$

(2.1)

Boundary conditions for regions rbc and p are either

$$\begin{array}{cc}c(0,t)={c}_{\text{in}}(t),& {c}_{x}(L,t)=0,\end{array}$$

(2.2)

or

$$\begin{array}{cc}\mathbf{\text{D}}{c}_{x}(0,t)=\mathbf{\text{U}}(c(0,t)-{c}_{\text{in}}(t)),& {c}_{x}(L,t)=0.\end{array}$$

(2.3)

Boundary conditions for regions ec, isf, and pc are based on reflections at *x*=0 and *x*=*L*

*c _{x}*(0,

Initial conditions for all five regions are

*c*(*x*,0) = 0.

The matrices **U, D**, and **K** are

**U** = diagonal(*u*_{rbc},*u*_{p},0,0,0,*u*_{rbc},*u*_{p},0,0,0),

**D** = diagonal(*D*_{rbc},*D*_{p},*D*_{ec},*D*_{isf},*D*_{pc},
${D}_{\text{rcb}}^{\ast},{D}_{\text{p}}^{\ast},{D}_{\text{ec}}^{\ast},{D}_{\text{isf}}^{\ast},{D}_{\text{pc}}^{\ast}$),

and

$$\mathbf{\text{K}}=\left[\begin{array}{ll}{\mathbf{\text{K}}}_{11}& {\mathbf{\text{K}}}_{12}\\ {\mathbf{\text{K}}}_{21}& {\mathbf{\text{K}}}_{22}\end{array}\right].$$

The first submatrix **K**_{11} is

$$\left[\begin{array}{ccccc}-({k}_{1}+{k}_{12})/{\mathit{\gamma}}_{\text{rbc}}& {k}_{12}/{\mathit{\gamma}}_{\text{rbc}}& 0& 0& 0\\ {k}_{12}& -{k}_{2}-{k}_{12}-{k}_{23}-{k}_{24}& {k}_{23}& {k}_{24}& 0\\ 0& {k}_{23}/{\mathit{\gamma}}_{\text{ec}}& -({k}_{3}+{k}_{23}+{k}_{34})/{\mathit{\gamma}}_{\text{ec}}& {k}_{34}/{\mathit{\gamma}}_{\text{ec}}& 0\\ 0& {k}_{24}/{\mathit{\gamma}}_{\text{isf}}& {k}_{34}/{\mathit{\gamma}}_{\text{isf}}& -({k}_{4}+{k}_{24}+{k}_{34}+{k}_{45})/{\mathit{\gamma}}_{\text{isf}}& {k}_{45}/{\mathit{\gamma}}_{\text{isf}}\\ 0& 0& 0& {k}_{45}/{\mathit{\gamma}}_{\text{pc}}& -({k}_{5}+{k}_{45})/{\mathit{\gamma}}_{\text{pc}}\end{array}\right].$$

The second submatrix is **K**_{12} = 0 and the third is

$${\mathbf{\text{K}}}_{21}=\text{diagonal}(-{k}_{1}/{\mathit{\gamma}}_{\text{rbc}}^{\ast},-{k}_{2},-{k}_{3}/{\mathit{\gamma}}_{\text{ec},}^{\ast}-{k}_{4}/{\mathit{\gamma}}_{\text{isf}}^{\ast},-{k}_{5}/{\mathit{\gamma}}_{\text{pc}}^{\ast}).$$

The fourth submatrix **K**_{22} is

$$[\begin{array}{ccccc}-{k}_{12}^{\ast}/{\mathit{\gamma}}_{\text{rbc}}^{\ast}& {k}_{12}^{\ast}/{\mathit{\gamma}}_{\text{rbc}}^{\ast}& 0& 0& 0\\ {k}_{12}^{\ast}& -{k}_{12}^{\ast}-{k}_{23}^{\ast}-{k}_{24}^{\ast}& {k}_{23}^{\ast}& {k}_{24}^{\ast}& 0\\ 0& {k}_{23}^{\ast}/{\mathit{\gamma}}_{\text{ec}}^{\ast}& -({k}_{23}^{\ast}-{k}_{34}^{\ast})/{\mathit{\gamma}}_{\text{ec}}^{\ast}& {k}_{34}^{\ast}/{\mathit{\gamma}}_{\text{ec}}^{\ast}& 0\\ 0& {k}_{24}^{\ast}/{\mathit{\gamma}}_{\text{isf}}^{\ast}& {k}_{34}^{\ast}/{\mathit{\gamma}}_{\text{isf}}^{\ast}& -({k}_{24}^{\ast}+{k}_{34}^{\ast}+{k}_{45}^{\ast})/{\mathit{\gamma}}_{\text{isf}}^{\ast}& {k}_{45}^{\ast}/{\mathit{\gamma}}_{\text{isf}}^{\ast}\\ 0& 0& 0& {k}_{45}^{\ast}/{\mathit{\gamma}}_{\text{pc}}^{\ast}& -{k}_{45}^{\ast}/{\mathit{\gamma}}_{\text{pc}}^{\ast}\end{array}].$$

To solve this system of equations, we split the main equation (2.1) into simpler subproblems, and solve each subproblem separately.

The following sections are a brief overview of a standard time-split method used in the red-blood cell five-region convection-diffusion-reaction model. Our aim is to give enough detail of the numerical methods to show our claim of second order accuracy.

We begin with a simple example that highlights the necessary error analysis required for an understanding of time-split methods. We follow with a section on the error analysis for our diffusion-convection-reaction equation. We discuss the boundary condition corrections that arise in our time-split method in the next section. Finally, we discuss the numerical methods used to solve each of the subproblems. More complete discussions of standard and widely used time-split methods can be found in LeVeque [17,18], Strikwerda [30], and Yanenko [31].

Consider the following Cauchy problem for the convection-reaction equation in one dimension:

$$\begin{array}{ll}{c}_{t}=-{\mathit{\text{Uc}}}_{x}+\mathit{\text{Kc}},& -\infty <x<\infty ,\\ c(x,{t}_{0})=f(x),\end{array}$$

(3.4)

where *U* and *K* are constant. This equation is of the general form

*c _{t}* =

where *A* depends on *c* and its spatial derivatives. To determine the solution at a later time, we use the exact solution operator

$$c({t}_{1})=S({t}_{1},{t}_{0})c({t}_{0}).$$

(3.5)

For the case we are looking at *t* does not appear explicitly, so the solution operator only depends on the time elapsed:

*S*(*t*_{1},*t*_{0}) = *S*(*t*_{1} − *t*_{0}).

If we use a two-level difference scheme, we have

$${c}^{n+1}=Q(k){C}^{n},$$

(3.6)

where *Q*(*k*) is an approximation to *S*(*k*), and *k* is the temporal step size. Eq. (3.6) is the difference analogue to

$$c({t}_{n+1})=S(k)c({t}_{n}),$$

(3.7)

The method (3.6) is accurate of order *p* if, for smooth functions *c*, the local truncation error (*Q*(*k*) − *S*(*k*))*c* is *O*(*k ^{p}*

We can derive the exact solution operator for the convection-reaction equation for the case when *c*^{r} and therefore **U,K** ^{r}^{×}^{r}

$${c}_{t}=-\mathbf{\text{U}}{c}_{x}+\mathbf{\text{K}}c.$$

(3.8)

We can expand the solution in a Taylor series about the point (*x,t*) to get

$$\begin{array}{l}c(x,t+k)=c(x,t)+{\mathit{\text{kc}}}_{t}(x,t)+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}(x,t)+\dots \\ =c+k(-\mathbf{\text{U}}{c}_{x}+\mathbf{\text{K}}c)+\frac{1}{2}{k}^{2}(-\mathbf{\text{U}}{c}_{\mathit{\text{xt}}}+\mathbf{\text{K}}{c}_{t})+\dots \\ =c+k(-\mathbf{\text{U}}{c}_{x}+\mathbf{\text{K}}c)+\frac{1}{2}{k}^{2}(-\mathbf{\text{U}}[-\mathbf{\text{U}}{c}_{\mathit{\text{xx}}}+\mathbf{\text{K}}{c}_{x}]+\mathbf{\text{K}}[-\mathbf{\text{U}}{c}_{x}+\mathbf{\text{K}}c])+\dots \\ =c+k(-\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+\mathbf{\text{K}})c+\frac{1}{2}{k}^{2}({\mathbf{\text{U}}}^{2}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}^{2}+[-\mathbf{\text{UK}}-\mathbf{\text{KU}}]\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+{\mathbf{\text{K}}}^{2})c+\dots \\ c(x,t+k)=\text{exp}(-k\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+k\mathbf{\text{K}})c(x,t).\end{array}$$

(3.9)

Thus we have for Eq. (3.8), the solution operator

$$S(k)=\text{exp}(-k\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+k\mathbf{\text{K}}).$$

(3.10)

Consider Eq. (3.8) again, but now split into two subproblems

$${c}_{t}=-\mathbf{\text{U}}{c}_{x},$$

(3.11)

$${c}_{t}=\mathbf{\text{K}}c.$$

(3.12)

Each of these subproblems has an exact solution operator. By a similar derivation for *S*(*k*) above, we can derive ones for (3.11) and (3.12) to get for (3.11), *S*_{1}(*k*) = exp(− *k***U** _{x}), and for (3.12), *S*_{2}(*k*) = exp(*k***K**).

The time split method is based on the fact that

$$S(k)\approx {S}_{2}(k){S}_{1}(k),$$

(3.13)

when *k* is small. In some cases, this splitting is exact. For example, for *c _{t}* =

To see the splitting error, we look at

$$\begin{array}{l}{S}_{2}(k){S}_{1}(k)-S(k)=\text{exp}(k\mathbf{\text{K}})\text{exp}(-k\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x})-\text{exp}(-k\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+k\mathbf{\text{K}})\\ =(\mathbf{\text{I}}+k\mathbf{\text{K}}+\frac{1}{2}{k}^{2}{\mathbf{\text{K}}}^{2}+\dots )\times (\mathbf{\text{I}}-k\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+\frac{1}{2}{k}^{2}{\mathbf{\text{U}}}^{2}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}^{2}+\dots )-\left(\mathbf{\text{I}}+k(-\mathbf{\text{U}}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+\mathbf{\text{K}})+\frac{1}{2}{k}^{2}({\mathbf{\text{U}}}^{2}\phantom{\rule{0.2em}{0ex}}{\partial}_{x}^{2}+[-\mathbf{\text{UK}}-\mathbf{\text{KU}}]\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+{\mathbf{\text{K}}}^{2})+\dots \right)\\ =\frac{1}{2}{k}^{2}(-\mathbf{\text{KU}}+\mathbf{\text{UK}})\phantom{\rule{0.2em}{0ex}}{\partial}_{x}+\text{O}({k}^{3}).\end{array}$$

(3.14)

It can be verified that the splitting is exact only if the matrices **K** and **U** commute.

For the four region model in Bassingthwaighte et al. [16]

$$\begin{array}{l}\mathbf{\text{U}}=\left[\begin{array}{llll}\mathit{\alpha}& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\\ \mathbf{\text{K}}=\left[\begin{array}{llll}{k}_{11}& {k}_{12}& {k}_{13}& 0\\ {k}_{21}& {k}_{22}& {k}_{23}& 0\\ {k}_{31}& {k}_{32}& {k}_{33}& {k}_{34}\\ 0& 0& {k}_{43}& {k}_{44}\end{array}\right],\end{array}$$

(3.15)

$$\begin{array}{l}\mathbf{\text{KU}}=\left[\begin{array}{llll}\mathit{\alpha}{k}_{11}& 0& 0& 0\\ \mathit{\alpha}{k}_{21}& 0& 0& 0\\ \mathit{\alpha}{k}_{31}& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right],\\ \mathbf{\text{UK}}=\left[\begin{array}{llll}\mathit{\alpha}{k}_{11}& \mathit{\alpha}{k}_{21}& \mathit{\alpha}{k}_{31}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right].\end{array}$$

(3.16)

Therefore, we see **KU** – **UK** ≠ 0. So the local error on smooth solutions is O(*k*^{2}) and hence, the splitting (3.11 – 3.12) is only first order accurate.

We can increase the order of the splitting to O(*k*^{3}) by using Strang splitting. Previously, to advance the solution one step in time

$${C}^{n+1}={Q}_{2}(k){Q}_{1}(k){C}^{n},$$

(3.17)

where *Q*_{1}(*k*) is the approximation to *S*_{1}(*k*), *Q*_{2}(*k*) is the approximation to *S*_{2}(*k*), and *C ^{n}* is the approximation to c at time

$${C}^{n+1}={Q}_{1}(k/2){Q}_{2}(k){Q}_{1}(k/2){C}^{n},$$

(3.18)

the O(*k*^{2}) terms drop out of Eq. (3.14) and we are left with a local truncation error that is O(*k*^{3}) which makes our splitting accurate to O(*k*^{2}). In Eq. (3.18), we show taking a half step with *Q*_{1} followed by a whole step with *Q*_{2}, followed by another half step with *Q*_{1}. When these steps are repeated over and over, the two half steps with *Q*_{1} can be combined into a single whole step to give

$$\begin{array}{l}{C}^{n+1}={Q}_{1}(k/2){Q}_{2}(k){Q}_{1}(k){Q}_{2}(k){Q}_{1}(k)\cdot \cdot \cdot {Q}_{2}(k)\times {Q}_{1}(k/2){C}^{0},\end{array}$$

(3.19)

which takes the solution from time *t*_{0} to *t _{n}*

Consider the constant coefficient diffusion-convection-reaction equation for five regions, the first two of which have convection.

$$\begin{array}{cc}{c}_{t}=\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}-\mathbf{\text{U}}{c}_{x}+\mathbf{\text{K}}c,& 0\le x\le L.\end{array}$$

(3.20)

Boundary conditions for regions 1–2

*c*(0,*t*) = *f*(*t*),

*c _{x}*(

Boundary conditions for regions 3–5

*c _{x}*(0,

*c _{x}*(

Initial conditions for regions 1–5

$$c(x,0)=0.$$

(3.21)

We can solve this problem using the following time-split scheme:

$${c}_{t}^{\u2020}=-\mathbf{\text{U}}{c}_{x}^{n},$$

(3.22)

$${c}_{t}^{\u2020\u2020}=\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}^{\u2020},$$

(3.23)

$${c}_{t}^{n+1}=\mathbf{\text{K}}{c}^{\u2020\u2020}.$$

(3.24)

We can then approximate Eqs. (3.22), (3.23) and (3.24) by finite differences to get

$${C}^{\u2020}={Q}_{1}{C}^{n},$$

(3.25)

$${C}^{\u2020\u2020}={Q}_{2}{C}^{\u2020},$$

(3.26)

$${C}^{n+1}={Q}_{3}{C}^{\u2020\u2020},$$

(3.27)

where *Q*_{1}, *Q*_{2}, and *Q*_{3} are the finite difference approximations to the operators in Eqs. (3.22), (3.23) and (3.24). To time step with Strang splitting as in Eq. (3.18), we use

$$\begin{array}{l}{C}^{n+1}={Q}_{1}(k/2){Q}_{2}(k/2){Q}_{3}(k)\times \underset{\text{repeats}}{\underset{\_}{{Q}_{2}(k/2){Q}_{1}(k){Q}_{2}(k/2){Q}_{3}(k)}}\cdot \cdot \cdot {Q}_{3}(k)\times {Q}_{2}(k/2){Q}_{1}(k/2){C}^{0}.\end{array}$$

(3.28)

For the four-region model without red blood cells [27] which is discussed further in the section on model utility, the error was

$${C}^{n+1}={Q}_{1}(k){Q}_{3}(k){Q}_{2}(k)\underset{\text{repeat}}{\underset{\_}{{Q}_{1}(k){Q}_{3}(k){Q}_{2}(k)}}\cdot \cdot \cdot {Q}_{1}(k)\times {Q}_{3}(k){Q}_{2}(k){C}^{0}.$$

(3.29)

The truncation error for the time split method of Eqs. (3.25), (3.26) and (3.27) can be represented as the sum of the trunctation errors in each subproblem plus the splitting error:

$${E}_{\text{total}}={E}_{\text{splitting}}+{E}_{\text{convection}}+{E}_{\text{diffusion}}+{E}_{\text{reaction}}.$$

(3.30)

Previously, this has been

$${E}_{\text{total}}=O(k)+0+O(k)+O({k}^{6})+O({h}^{2}).$$

(3.31)

For the proposed five-region rbc model, this error would be

$${E}_{\text{total}}=O({k}^{2})+O({k}^{p})+O({k}^{2})+O({k}^{q})+O({h}^{2}).$$

(3.32)

The methods chosen to solve the convection subproblem and the reaction subproblem would be such that *p* and *q* would be at least equal to 2. This would give an overall truncation error of O(*k*^{2} + *h*^{2}), which would be an improvement over the error for the four-region model of O(*k* + *h*^{2}).

For the splitting scheme given in Eqs. (3.22), (3.23) and (3.24), we solve exactly each of the subproblems

$${c}_{t}^{\u2020}=-\mathbf{\text{U}}{c}_{x}^{\u2020},$$

(3.33)

$${c}_{t}^{\u2020\u2020}=\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}^{\u2020\u2020},$$

(3.34)

$${c}_{t}^{\u2020\u2020\u2020}=\mathbf{\text{K}}{c}^{\u2020\u2020\u2020}.$$

(3.35)

However, we need to use modified boundary conditions for each subproblem, and not the boundary conditions for the original problem (3.20).

For the convection subproblem (3.33), we need to specify boundary conditions *c*^{†}(0,*t _{n}*

$${c}^{\u2020}(0,{t}_{n+1})={c}^{\u2020}(0,{t}_{n})+k{c}_{t}^{\u2020}(0,{t}_{n})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.36)

We use Eq. (3.33) to replace ${c}_{t}^{\u2020}(0,{t}_{n})$ to get

$${c}^{\u2020}(0,{t}_{n+1})={c}^{\u2020}(0,{t}_{n})-k\mathbf{\text{U}}{c}_{x}^{\u2020}(0,{t}_{n})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.37)

We also have for the initial conditions *c*^{†}(*x,t _{n}*) =

$${c}^{\u2020}(0,{t}_{n+1})=c(0,{t}_{n})-k\mathbf{\text{U}}{c}_{x}(0,{t}_{n})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.38)

Now we can use the original p.d.e. (3.20) to replace the **U***c _{x}* term to get

$${c}^{\u2020}(0,{t}_{n+1})=c(0,{t}_{n})+k({c}_{t}(0,{t}_{n})-\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})-\mathbf{\text{K}}{c}_{x}(0,{t}_{n}))+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.39)

Since

$$c(0,{t}_{n+1})=c(0,{t}_{n})+k{c}_{t}(0,{t}_{n})+\cdot \cdot \cdot ,$$

we can simplify Eq. (3.39) to get

$${c}^{\u2020}(0,{t}_{n+1})=c(0,{t}_{n+1})-k(\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})+{\mathbf{\text{K}}}_{c}(0,{t}_{n})).$$

(3.40)

We see from this expression that we have a modified form of the boundary condition which includes a correction for diffusion and reaction.

Now consider the diffusion equation subproblem (3.34). If we solve this equation by a stable implicit method like Crank–Nicolson, we need modified boundary conditions at time *t*_{n+1}. The initial conditions will be *c*^{††}(*x,t _{n}*) =

$${c}^{\u2020\u2020}(0,{t}_{n+1})={c}^{\u2020\u2020}(0,{t}_{n})+k{c}_{t}^{\u2020\u2020}(0,{t}_{n})+\frac{1}{2}{k}^{2}{c}_{u}^{\u2020\u2020}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.41)

Using the subproblem (3.34) to replace the first time derivative, we get

$${c}^{\u2020\u2020}(0,{t}_{n+1})={c}^{\u2020\u2020}(0,{t}_{n})+k\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}^{\u2020\u2020}(0,{t}_{n})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020\u2020}(0,{t}_{n})+\cdot \cdot \cdot .$$

(3.42)

Our initial conditions allow us to substitute for

$${c}^{\u2020\u2020}(0,{t}_{n+1})={c}^{\u2020}(0,{t}_{n+1})-k\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}^{\u2020}(0,{t}_{n+1})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n+1})+\cdots .$$

(3.43)

Now we can use Eq. (3.40) to replace *c*^{†}(0,*t _{n}*

$${c}^{\u2020\u2020}(0,{t}_{n+1})=c(0,{t}_{n+1})-k\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})-k\mathbf{\text{K}}c(0,{t}_{n})+k\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}^{\u2020}(0,{t}_{n+1})+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n+1})+\cdot \cdot \cdot .$$

(3.44)

From Eq. (3.40), we can differentiate twice with respect to *x* to get

$${c}_{\mathit{\text{xx}}}^{\u2020}(0,{t}_{n+1})={c}_{\mathit{\text{xx}}}(0,{t}_{n+1})-k(\mathbf{\text{D}}{c}_{\mathit{\text{xxxx}}}(0,{t}_{n})+\mathbf{\text{K}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})).$$

This expression can then be substituted for ${c}_{\mathit{\text{xx}}}^{\u2020}(0,{t}_{n+1})$ in Eq. (3.44) to give

$${c}^{\u2020\u2020}(0,{t}_{n+1})=c(0,{t}_{n+1})+k\mathbf{\text{D}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})-k\mathbf{\text{K}}c(0,{t}_{n})-k\mathbf{\text{D}}({c}_{\mathit{\text{xx}}}(0,{t}_{n+1})-k(\mathbf{\text{D}}{c}_{\mathit{\text{xxxx}}}(0,{t}_{n})+\mathbf{\text{K}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})))+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n+1})+\cdot \cdot \cdot .$$

This expression simplifies to

$${c}^{\u2020\u2020}(0,{t}_{n+1})=c(0,{t}_{n+1})-k\mathbf{\text{K}}c(0,{t}_{n})+k\mathbf{\text{D}}(-k(\mathbf{\text{D}}{c}_{\mathit{\text{xxxx}}}(0,{t}_{n})+\mathbf{\text{K}}{c}_{\mathit{\text{xx}}}(0,{t}_{n})))+\frac{1}{2}{k}^{2}{c}_{\mathit{\text{tt}}}^{\u2020}(0,{t}_{n+1})+\cdot \cdot \cdot .$$

(3.45)

With the modified boundary conditions (3.40) and (3.45) and the fact that the initial conditions for Eq. (3.35) are *c*^{†††}(*x,t _{n}*) =

By splitting the problem into three subproblems, we can use different numerical methods to solve each subproblem. We use a combination of exact and random choice for the convection subproblem, Crank–Nicolson for the diffusion subproblem, and an o.d.e solver that automatically chooses between Adams–Moulton and Gear's method for the reaction subproblem.

Since the convection subproblem (3.33) is two separate one-way wave equations, we can select the Δ*t* and Δ*x* such that we can use the exact solution for one of the regions. We choose the normally slower region, that of the plasma, for the exact solution, and use Δ*x* = *L/N*_{seg} and Δ*t* = Δ*x/u*_{p}. To update the solution in the plasma region, we just do a shift

$$\begin{array}{ll}{C}_{j}^{n+1}={C}_{j-1,}^{n}& j=1,2,\dots ,N.\end{array}$$

For *j* = 0, we use the modified boundary condition as given by Eq. (3.40). In the rbc region, which is the higher velocity region, we use the random choice method as described by Colella [19]. To update the solution at the next time step *t _{n}*

$${C}_{j}^{n+1}={C}_{\text{exact}}^{n}((j-\frac{1}{2}+{a}^{n+1})\Delta x,{t}_{n+1}),j=1,2\dots ,N,$$

where
${C}_{\text{exact}}^{n}$ is the solution to the Riemann problem and a^{n+1} is a van der Corput sampling sequence (see Ref. [19] for complete details). Since we assume that the fast rbc velocity *u*_{rb}* _{c}*, is in the range

*u*_{rbc} = (*u* _{p} + *u*_{p} + *u*_{p} + 2*u*_{p})/4 = 125*u*_{p}.

Further modifications have to be made to the boundary conditions. When the updating is
${C}_{j}^{n+1}={C}_{j-1}^{n}$, we only need to specify
${C}_{0}^{n+1}$ which is given by Eq. (3.40). However, when the updating is
${C}_{j}^{n+1}={C}_{j-2}^{n}$, we need to specify both
${C}_{0}^{n+1}$ and
${C}_{1}^{n+1}$. For
${C}_{0}^{n+1}$, we would use Eq. (3.40) as before. For
${C}_{1}^{n+1}$, we to correct the time step in Eq. (3.40), given by *k*, which is too large. The actual time increment is between 0 and *k*/2 depending on the speed *u*_{rbc}. If *u*_{rb}* _{c}* =

We solve the diffusion subproblem by using a standard Crank-Nicolson approach. The details and the code may be requested from the authors. This allows us to use standard linear algebra libraries [29] to solve the resulting linear system of equations.

To solve the reaction subproblem, we use the ordinary differential equation solver SDRIV2 from Kahaner et al. [20]. This solver allows for the dynamic switching between Adams–Moulton and Gear's method depending on the stiffness of the problem. The accuracy of the solution is chosen so that the errors from solving the reaction subproblem are not greater than the other errors in Eq. (3.30).

After the subproblems are solved, we have
${C}_{j}^{n+1}$ for *j* = 2,…,*N* − 1. For
${C}_{0}^{n+1}$, we just use either Eq. (2.2) or Eq. (2.3). If *u*_{rbc} > *u*_{p}, then we also need to correct
${C}_{1}^{n+1}$. Whatever material is flowing in from the boundary undergoes not only convection, but also reaction and diffusion.
${C}_{1}^{n+1}$ needs to be corrected for one time step's worth of reaction and diffusion, using Eq. (3.40) and Eq. (3.45), otherwise conservation of mass will be violated. Without this correction, the solution ends up having too much material flowing in from the boundary during those updating steps when both
${C}_{0}^{n+1}$ and
${C}_{1}^{n+1}$ are updated from the boundary.

Figs. 2 and and33 show the results from the rbc model with the splitting scheme given in Eqs. (3.22), (3.23) and (3.24). The numerical values for the parameters are listed in the legend of Fig. 2. These values were choosen to test the numerical method. Fitting the model to experimental data is given in the next section.

Concentration–time profiles in the input (left panel insert), outflow (right panel), and concentration–space profile in the rbc region at *t* = 18 s (left panel). The following values were used for the parameters: *F*_{blood} = 2.5 ml g^{−l} **...**

Fig. 2 shows in the left panel the position of the tracer oxygen and tracer water concentration distance profiles at *t* = 18 s, just 3 s after the input function drops to zero. Because of the unrealistically large value of 1.0 ml g^{−1} chosen for *V*_{cap}, the capillary transit time is *V*_{cap}/*F*_{blood} = 24 s, so the bolus front is only 3/4 of the way along the capillary. The front of the pulse at *x/L* = 0.75 is blunted, more so than is the trailing edge, as expected. The outflow curves reveal the separation in concentrations of erythrocytes and plasma, due to the resistance at the red cell membrane and the loss of oxygen from plasma into tissue.

Fig. 3, using the same parameters, illustrates the profiles in the five regions at 12 s, just as the pulse front reaches the capillary midpoint. As was also evident in Fig. 2, there are finite concentrations of oxygen beyond the mean bolus front because of the axial diffusion. The relative magnitude of the the oxygen concentration in the plasma as compared to the concentration in red blood cells is due to the artificially high values for PS_{rbc} used for testing. The concentrations of tracer water, produced by oxygen consumption, are highest in the interstitium because the production by transformation was occurring in endothelial cells and interstitium, consuming most of the substrate oxygen before it reached the parenchymal cells.

To calculate an overall order of accuracy of the rbc model, we look at the error in output concentrations for several values of *N*_{seg}. Fig. 4 shows the output concentration profiles in both the red blood cell region and the plasma region for the cases *N*_{seg} = 5, 10, 20, 40, 80, 160, 320. The 2-norm of the output concentrations is shown in Table 1. We can see that as *N*_{seg} increases, the error becomes progressively smaller. If we assume asymptotic error

Concentration profiles of O_{2} and H_{2}O in the outflow of the red blood cell region and the plasma region for *N*_{seg} = 10, 40, 160, 320.

$${\Vert {C}_{\text{out}}^{{N}_{\text{seg}}}-{C}_{\text{out}}^{{N}_{\text{seg}}=320}\Vert}_{2}=\text{constant}{(1/{N}_{\text{seg}})}^{p}$$

is valid then the overall order *p*, from a least squares fit of a straight line to the log-transformed 2-norm of the error is 1.85.

Fig. 5 shows the results of fitting the model to experimental data. The pair of ^{15}O-concentration–time curves were recorded at the inflow and outflow to an isolated blood-perfused rabbit heart preparation. The perfusate was 20% hematocrit with erythrocytes suspended in standard Krebs-Henseleit physiological saline. A bolus of 1.4 ml of blood containing ^{15}O-oxygen was injected into the inflow stream via a switching arrangement and the radioactivity in the inflow and outflow recorded via flow-through plastic scintillation detectors for the *β* + -emissions were counted over 1 s intervals. For the analysis, the input to the blood-tissue exchange model was the observed upstream dilution curve. The flow used for the modeling was the measured flow. The endothelial cells were assumed to have a negligible contribution to the waveform and the consumption of oxygen and so those parameters were set to zero. Consequently, the estimates of PS_{g} represent the sum of permeation through interendothelial clefts plus permeation across the endothelial cells. The volumes of distribution were estimated from the water spaces in cach region, for example, from the steady-state tracer estimates of Gonzalez and Bassingthwaighte [21], modified to provide for the binding of oxygen by hemoglobin in red blood cells and by myoglobin in the myocytes (parenchymal cells). Values for PS_{rbc} were set to high values in accord with the observations of Forster et al. [22]. The ratio of small vessel to large vessel hematocrit, Hct_{ratio}, was set to 0.98 for this solution. The axial diffusion coefficients in all regions were set to an arbitrarily high value of 10^{−4} cm^{2} s^{−1} to degrade axial gradients to some extent as would diffusional exchanges between capillary-tissue units with offset beginnings and endings. The adjustable or free parameters were therefore only PS_{g} and PS_{pc} for water and oxygen, and the consumption or transformation parameter, *G*_{pc}, for oxygen within myocytes.

Outflow concentration-time curves for experimental data on ^{15}O-tracer in the effluent from an isolated perfused heart followed bolus injection and model. Upper panel: Input concentration-time curve for ^{15}O-oxygen in blood. Lower panel: Outflow dilution **...**

Fig. 6 shows the pronounced effect of the differing velocities between plasma and red blood cells. In this figure, we see for a Hct_{ratio} = 0.7, the velocity in the plasma is *u*_{p} = 0.01573 cm s^{−1} and for the red blood cells the velocity is *u*_{rbc}= 0.02415 cm s^{−1}. As the Hct_{ratio} increases, these two velocities become more nearly equal until at Hct_{ratio} = 1.0, they are the same. The values of the other parameters are the same as those used in Fig. 5 except all the PS* _{i}* = 0. Erthrocyte precession exacerbates early tracer loss from erythrocytes into tracer-free plasma (left curves). This does not occur when both travel at the same velocity (the rightmost curve).

Outflow concentration–time curves for oxygen within red blood cells for various hematocrit ratios. For these solutions, the PS's are all zcro, so the curves represent intravascular transport only. Whcn erythrocytes travel faster than plasma, more **...**

Fig. 7 shows the simulated effect of various hemoglobin binding values on Xenon transfer. The various values of the partition coefficient *θ*_{rbc} were chosen to be similar to that found experimentally. See for example, Carlin and Chien [23] or Goresky et al. [24]. Higher *θ*_{rbc}, equivalent to higher partitioning into erythrocytes, gives curves closer to erythrocyte curves, low *θ*_{rbc} gives waterlike curves. The change of *θ* changes the blood-tissue partition coefficient. Because the PS's used for these solutions are so high, the tracer is virtually flow-limited in its exchange. This means that they will be almost superimposed on each other by time scaling each relative to its own mean transit time [25,26], which is a simple but revealing test for flow-limited blood-tissue exchange.

The numerical methods are designed for high accuracy. Related models to which this model can be reduced, e.g., those of Bassingthwaighte et al. [27], with three regions and Bassingthwaighte et al. [16], with four regions, are apparently accurate also, since they can be fitted very closely by this model and have been tested against analytical models. The analytical models, either three-region ([28], water model) or four-region [16], are not at all accurate at long times: the latter, for example, requires two convolution integrations and one double convolution integration, all containing slowly convergent Bessel functions. Thus, we regard this new model as a standard for accuracy.

A main reason for the numerical strategy is that the scheme lends itself to further developments in treating nonlinear systems. For solutes at finite chemical levels, where the linearization due to the use of tracers at effectively zero chemical concentrations is no longer valid, then the PS's and *G*'s and even the volumes of distribution become concentration-dependent. In this more demanding situation, the basic numerical method can be anticipated to work well.

We have developed a second order accurate computational model of oxygen transport and reaction from red blood cells to mitochondria. Our model treats the flow of red blood cells and the flow of plasma separately. The model is useful for analyzing residue data from positron emission tomography or outflow data from multiple indicator dilution experiments.

We would like to thank Joseph Chan and Gary Raymond for their helpful discussions concerned with the development of this model. Research has been supported by NIH grants RR 1243 (Simulation Resource Facility: http://nsr.bioeng.washington.edu/) and HL38736 and by training grant T32 HL07403. Work supported by NIH grants HL38736 and RR1243.

1. Krogh A. The number and distribution of capillaries in muscles with calculation of the oxygen pressure head necessary for supplying the tissue. J Physiol (Lond) 1919;52:409–415. [PubMed]

2. Gutierrez G. The rate of oxygen release and its effect on capillary O_{2} tension: a mathematical analysis. Respir Physiol. 1986;63:79–96. [PubMed]

3. Federspiel WJ, Sarelius IH. An examination of the contribution of red cell spacing to the uniformity of oxygen flux at the capillary wall. Microvasc Res. 1984;27:273–285. [PubMed]

4. Federspiel WJ, Popel AS. A theoretical analysis of the effect of the particulate nature of blood on oxygen release in capillaries. Microvasc Res. 1986;32:164–189. [PubMed]

5. Ellsworth ML, Popel AS, Pittman RN. Assessment and impact of heterogeneties of convective oxygen transport parameters in capillaries of striated muscle: experimental and theoretical. Microvasc Res. 1988;35:341–362. [PubMed]

6. Groebe K, Thews G. Theoretical analysis of oxygen supply to contracted skeletal muscle. Adv Exp Med Biol. 1986;200:495–514. [PubMed]

7. Groebe K, Thews G. Effects of red cell spacing and red cell movement upon oxygen release under conditions of maximally working skeletal muscle. Adv Exp Med Biol. 1989;248:175–185. [PubMed]

8. Hellums JD. The resistance to oxygen transport in the capillaries relative to that in the surrounding tissue. Microvasc Res. 1977;13:131–136. [PubMed]

9. Popel AS. Theory of oxygen transport to tissue. Crit Rev Biomed Eng. 1989;17:257–321. [PubMed]

10. Rose CP, Goresky CA. Limitations of tracer oxygen uptake in the canine coronary circulation. Circ Res. 1985;56:57–71. [PubMed]

11. Groebe K, Thews G. Role of geometry and anisotropic diffusion for modeling *p*O_{2} profiles in working red muscle. Respir Physiol. 1990;79:255–278. [PubMed]

12. Deussen A, Bassingthwaighte JB. Modeling ^{15}O-oxygen tracer data for estimating oxygen consumption. Am J Physiol Heart Circ Physiol. 1996;270:H1115–H1130. [PMC free article] [PubMed]

13. Mintun MA, Raichle ME, Martin WRW, Herscovitch PP. Brain oxygen utilization measured with O-15 radiotracers and positron emission tomography. J Nucl Med. 1984;25:177–187. [PubMed]

14. Huang SC, Phelps ME. Error in parameter estimates with variations in flow: measurement of oxygen utilization with positron-emission tomography. Circulation. 1985;72(Suppl IV):77–80. [PubMed]

15. Larson KB, Markham J, Raichle ME. Tracer-kinetic models for measuring cerebral blood flow using externally detected radiotracers. J Cereb Blood Flow Metab. 1987;7:443–463. [PubMed]

16. Bassingthwaighte JB, Wang CY, Chan IS. Blood-tissue exchange via transport and transformation by capillary endothelial cells. Circ Res. 1989;65:997–1020. [PMC free article] [PubMed]

17. LeVeque RJ, Oliger J. Numerical methods based on additive splittings for hyperbolic partial differential equations. Math Comput. 1983;40:469–497.

18. LeVeque RJ. Intermediate boundary conditions for time-split methods applied to hyperbolic partial differential equations. Math Comput. 1986;47:37–54.

19. Colella P. Glimm's method for gas dynamics. SIAM J Sci Stat Comput. 1982;3(1)

20. Kahaner D, Moler C, Nash S. Numerical Methods and Software. Prentice-Hall; New Jersey: 1989.

21. Gonzalez F, Bassingthwaighte JB. Heterogeneities in regional volumes of distribution and flows in the rabbit heart. Am J Physiol Heart Circ Physiol. 1990;258:H1012–H1024. [PMC free article] [PubMed]

22. Forster RE, Goodwin CW, Itada N. A new approach to the measurement of mean tissue *p*O_{2}. In: Grote J, et al., editors. Oxygen Transport to Tissue II. Plenum Press; New York: 1976. pp. 41–46.

23. Carlin R, Chien S. Partition of xenon and iodoantipyrine among erythrocytes, plasma, and myocardium. Circ Res. 1977;40:497–504. [PubMed]

24. Goresky CA, Schwab AJ, Rose CP. Xenon handling in the liver: red cell capacity effect. Circ Res. 1988;63:767–778. [PubMed]

25. Goresky CA. A linear method for determining liver sinusoidal and extravascular volumes. Am J Physiol. 1963;204:626–640. [PubMed]

26. Bassingthwaighte JB. Physiology and theory of tracer washout techniques for the estimation of myocardial blood flow: flow estimation from tracer washout. Prog Cardiovasc Dis. 1977;20:165–189. [PMC free article] [PubMed]

27. Bassingthwaighte JB, Chan IS, Wang CY. Computationally efficient algorithms for convection-permeation-diffusion models for blood-tissue exchange. Ann Biomed Eng. 1992;20:687–725. [PubMed]

28. Rose CP, Goresky CA, Bach GG. The capillary and sarcolemmal barriers in the heart: an exploration of labeled water permeability. Circ Res. 1977;41:515–533. [PubMed]

29. Anderson E, Bai Z, Bischof C, Demmel J, Dongarra J, Du Croz J, Greenbaum A, Hammarling S, McKenney A, Ostrouchov S, Sorensen D. LAPACK Users Guide. second. SIAM; Philadelphia, PA: 1995.

30. Strikwerda JC. Finite Difference Schemes and Partial Differential Equations. Wadsworth and Brooks/Cole; Pacific Grove, CA: 1989.

31. Yanenko NN. The Method of Fractional Steps; the Solution of Problems of Mathematical Physics in Several Variables. Springer; Berlin: 1971.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |