## Abstract

A mathematical method is presented for analysing oxygen transport in skeletal muscle that includes both the interaction of very large number of capillaries and the effect of axial diffusion in the tissue. The analysis takes into account differences in the blood flow velocity among the capillaries and in the oxygen concentration and haematocrit of the blood entering the capillaries at the arterial end. First, the oxygen concentration is determined within a functional unit consisting of a single capillary surrounded by a region of tissue, in which a flux is prescribed on the outer boundary of the region. The method of matched asymptotic expansions is used to include the effect of axial diffusion in the tissue. The flux that is prescribed on the boundary of the functional unit is a result of the interaction among all the capillaries comprising the vascular bed. These fluxes are found by matching the oxygen concentration along the boundaries of adjacent functional units. This leads to a system of ordinary differential equations for the capillary concentrations coupled with a system of algebraic equations for the fluxes. Results are presented for a two-dimensional array of capillaries when each capillary has a different blood flow velocity. Axial diffusion causes the oxygen concentration to drop more rapidly at the arterial end and remain below the level determined when axial diffusion is neglected. This difference decreases toward the venous end, because the excess oxygen that left the capillaries at the arterial end and diffused axially is available to the tissue there. In some cases a single capillary model, which neglects interaction, incorrectly indicates that portions of the tissue become anoxic, whereas including interaction shows that all of the tissue is in fact supplied with oxygen.

## 1. Introduction

Oxygen is transported by blood dissolved in plasma and bound to haemoglobin, and diffuses from the capillaries into the tissue where it is consumed. In skeletal muscle the capillaries run parallel to the muscle fibres, and the blood flow velocity within the capillaries and the oxygen concentration and haematocrit of the blood entering the capillaries at the arterial end vary throughout the muscle. The tissue is therefore not uniformly perfused and a large-scale diffusion of oxygen occurs both in the plane perpendicular to the capillaries, and in the direction parallel to the capillaries, as oxygen moves from regions that are richly supplied to regions that are poorly supplied. Because of the parallel arrangement of the capillaries, the mathematical analysis of oxygen transport in skeletal muscle is considerably simplified by using a conceptual model due to Krogh (1919), consisting of a single capillary supplying oxygen to a surrounding concentric cylinder of tissue. Krogh simplified the analysis of his model by assuming that diffusion parallel to the capillary, axial diffusion, can be neglected, and considered only diffusion in the radial direction. This assumption was motivated by the fact that the capillary spacing is small compared with the capillary length, and so the tissue cylinder is long and thin. The result is a mathematical problem that is trivial to analyse but provides some interesting information on oxygen delivery within skeletal muscle. Krogh's initial work has been extended in many ways to make it more realistic and to examine a variety of phenomena (see, for example the reviews by Fletcher 1978; Popel 1989; Hoofd 1992). Using numerical analysis, Whiteley *et al*. (2002) have shown that axial diffusion is significant within the tissue. Including axial diffusion leads to a much more complex mathematical problem, which has been treated analytically by Salathe *et al*. (1980) using matched asymptotic expansions.

While these single capillary models provide useful information on oxygen transport to tissue, they are inherently limited because they do not take into account the large-scale interaction among all the capillaries comprising the vascular bed. The analysis can be extended to include multiple capillaries by writing an equation for each capillary and for the surrounding tissue. This leads to a formidable mathematical problem when more than a few capillaries are included. Numerical methods have been used with a moderate number of capillaries (Popel 1978). Other approximate methods have been developed to include multiple capillaries (Hoofd 1995*a*,*b*; Titcombe & Ward 2000), but do not include axial diffusion. As an alternative, a homogenization approach was developed (Salathe 1982) in which a single nonlinear partial differential equation described convection by the capillaries, and diffusion and consumption within the tissue, and this has been used to study a number of problems of physiological interest (Salathe & Xu 1988; Salathe & Chen 1993).

Other than the homogenization approach, there is at present no effective method to analyse the interaction of a large number of capillaries when axial diffusion is included. A relatively simple approach to this problem is presented here which, unlike the homogenization method, preserves the discrete structure of the capillaries within the tissue. This is accomplished by combining an earlier analysis of axial diffusion in a single capillary model which uses the method of matched asymptotic expansions (Salathe *et al*. 1980), with a new technique that has been developed to study the interaction among many capillaries in a vascular bed, but which did not include axial diffusion (Salathe 2003). The method consists of associating each capillary with a region of tissue, in this case a box of square cross-section with the capillary located along the centre line. However, unlike the Krogh model, a flux is prescribed along the outer edge of the box because the oxygen supplied by the capillary at its centre line is consumed not only within the box and because the tissue in the box is supplied by capillaries other than the one at its centre line. These fluxes, which account for the interaction among the capillaries, are determined by matching the concentrations along the boundaries of adjacent boxes.

The first step is to solve the equations governing oxygen concentration, with axial diffusion included, within a single box in which a flux is prescribed along each of the four sides, then from this solution to determine the concentration along each side. The method is simplified by assuming that the flux along each of the four sides of the box are functions only of axial location along the box. They are found by matching, at each axial location, the average concentration at the boundary of adjacent boxes. Although the analysis is rather involved, and the resulting solutions are complex, the average concentrations are determined as linear functions of the four fluxes prescribed on the boundaries, with coefficients that depend only on the ratio of the capillary radius to the capillary spacing. Now that these coefficients are known, there is no need to repeat the analysis that led to them. A single ordinary differential equation governs the concentration in each capillary. These equations contain the fluxes along the sides of the box associated with the capillary, and are coupled by a system of algebraic equations for the fluxes in terms of the capillary concentrations.

The exact approach would require letting the fluxes be functions of position along each of the four boundaries at any axial location, and finding them by matching concentrations point-wise along the boundaries of adjacent units. This results in a very complex analysis, even without axial diffusion, and involves expansions in terms of non-orthogonal eigenfunctions and the inversion of infinite matrices. Several such problems have been solved neglecting axial diffusion and the exact solution compared with the solutions obtained with the approximate method (Salathe 2003). The error inherent in the approximate method was less than 1% in some cases, a few per cent in others, and no more than 5% in any of the cases studied. This is a reasonable level of accuracy considering the relatively simple approach the approximate method yields to analyse a rather difficult problem.

## 2. Analysis

A capillary of length *L* and radius *R*_{c} is located along the centre line of a box of tissue that has sides of width 2*R*_{t} and length *L*. The *z*-axis is along the centre line of this box, the *x*- and *y*-axes are parallel to the sides of the box cross-section, and the origin is at the arterial end. The distances *x* and *y* are normalized by *R*_{t} and *z* is normalized by *L*. Blood enters the capillary at the arterial end with velocity *U* and oxygen concentration *C*_{a}. As oxygen is convected through the capillary, it diffuses into the surrounding tissue where it is consumed at a constant rate *M*_{D}. Oxygen also enters the tissue across the outer boundaries of the box. At any cross-section perpendicular to the capillary, the average flux per unit area across each of the four sides of the box—right (*x*=1), left (*x*=−1), top (*y*=1) and bottom (*y*=−1)—is denoted by *F*_{R}, *F*_{L}, *F*_{T} and *F*_{B}, respectively, where *F*_{B} and *F*_{L} are fluxes into the box and *F*_{R} and *F*_{T} are fluxes out of the box. These fluxes are functions of axial location *z* only, and are non-dimensional quantities normalized by , where is a reference concentration and *D*_{t} is the oxygen diffusivity of the tissue. The oxygen concentration in the tissue, *c*_{t}, is a function of *x*, *y* and *z*, and is normalized by . Within the capillary, the red blood cells are separated by a bolus of fluid undergoing a recirculating motion, which results in convective mixing and a more uniform oxygen concentration at any cross-section. Therefore, the oxygen concentration in the capillary, *C*_{b}, is assumed to be a function only of axial location *z* and is normalized by where *α*=α_{b}/*α*_{t} is the ratio of the solubility of oxygen in tissue to its solubility in blood.

The oxygen concentration within the tissue, *c*_{t}(*x*, *y*, *z*), is governed by the equation(2.1)where and *ϵ*=*R*_{t}/*L*. The oxygen flux specified at the outer boundary of the box is taken to be the average value along each of the four sides at any axial location *z*. There is no oxygen flux out of the ends of the box, *z*=0, 1. Therefore, in dimensionless form, the boundary conditions along the sides of the box are(2.2)and at the ends the boundary conditions are(2.3)

At the capillary wall, the oxygen concentration in the tissue and in the capillary are related by(2.4)where *R*=*R*_{c}/*R*_{t}. The rate of decrease in oxygen concentration along the capillary is related to the rate oxygen diffuses into the tissue, according to the equation(2.5)where and *θ*=arctan(*y*/*x*) are the polar coordinates. The oxygen capacity of blood, *N*_{D}, is proportional to the haematocrit and , *D*=2*D*_{t}*L*/*αUR*_{c}*R*_{t}. The oxyhaemoglobin dissociation relationship is given by(2.6)for suitable constants *K*^{*} and *p*. This is one of the most commonly used forms of the oxyhaemoglobin dissociation relationship, although it loses accuracy at low concentration. A discussion of alternative forms has been given by Meldon (1987). The capillary concentration at *z*=0 satisfies the initial condition(2.7)

The rate at which oxygen within the capillary diffuses along the length of the capillary is negligible compared with the rate at which it is convected, and so axial diffusion is neglected in equation (2.5).

This system of equations will be solved for *R*≪1, *ϵ*≪1 using the method of matched asymptotic expansions. Writing(2.8)the zero order terms satisfy the equation(2.9)with boundary conditions(2.10)and the matching condition at the capillary(2.11)

Because axial diffusion, and therefore the term in equation (2.1), has been neglected, it is not possible for the zero order equations to satisfy the no-flux boundary conditions at the ends *z*=0 and *z*=1, equation (2.3). Equation (2.8) therefore represents a singular perturbation expansion, and is valid only in the region bounded away from the ends *z*=0 and *z*=1. The solution to equations (2.9)–(2.11) has been obtained previously (Salathe 2003) and is(2.12)where(2.13)(2.14)and , , , .

This solution satisfies the boundary condition at the capillary wall, equation (2.4), only to *O*(*R*), and so is valid for *R*≪1. Therefore, the strength of the logarithmic source in equation (2.14) is , which equals the net flux out through the edges of a square with sides of length 2 and consumption 4*M*, corresponding to the total area of the square, and not this area diminished by the area *πR*^{2} of the capillary at the centre.

The first order quantities satisfy the equation(2.15)with boundary conditions(2.16)and the matching condition at the capillary(2.17)

Again, the zero flux boundary conditions at *z*=0 and *z*=1 are not satisfied, because the expansion, equation (2.8), is not valid in the neighbourhood of *z*=0 and *z*=1. The solution to these first order equations, for *R*≪1, is (see Appendix A for details)(2.18)where(2.19)and , , , , , .

From these solutions, the average concentration along the four boundaries of the region at any axial location *z* can be found. For the zero order solution these are(2.20)and for the first order solutions they are(2.21)where(2.22)

This functional unit and the capillary associated with it are now regarded as one of many interacting units forming a vascular bed. A cross-section perpendicular to the capillaries is assumed to be such that it can be divided into rows and columns to form squares with each square containing one capillary; and so that the capillaries, which are at the centre of these squares, are aligned in both the *x* and *y* directions. The functional units are also arranged so that the arterial ends all start in the same plane, perpendicular to the capillaries. A heterogeneous distribution of capillaries in the *x*−*y* plane can be modelled by extending the analysis of the functional unit to the case where the capillary does not lie along its centre line. There are *N* rows and *M* columns of these units, and *C* ^{n,m} denotes the concentration *C*_{b} in the capillary occupying the *n*th row and *m*th column. The flux in the *x*‐direction out of the unit in the *n*th row and *m*th column is denoted by *F*^{n,m}, and the flux in the *y*‐direction into this unit is denoted by *G*^{n,m}. Therefore, for the unit in the *n*th row and *m*th column, *F*_{R}=*F*^{n,m}, *F*_{L}=*F*^{n,m−1}, *F*_{T}=*G*^{n−1,m}, *F*_{B}=*G*^{n,m}. Matching the average concentration along the edge *x*=1 of the unit in the *n*th row and *m*th column to the average concentration along the edge *x*=−1 of the unit in the *n*th row and (*m*+1)th column gives, to zero order in *ϵ*,(2.23)Matching the average concentration along the edge *y*=−1 of the unit in the *n*th row and *m*th column to the average concentration along the edge *y*=1 of the unit in the (*n*+1)th row and *m*th column gives, to zero order in *ϵ*,(2.24)

These equations determine the zero order fluxes in terms of the zero order capillary oxygen concentrations. Corresponding matching to first order in *ϵ* gives equations for the first order fluxes in terms of the first order capillary oxygen concentrations and the zero order fluxes and concentrations. Matching in the *x*-direction gives(2.25)while matching in the *y*-direction gives(2.26)

To determine the capillary concentrations, equation (2.5) is written in the form(2.27)where ℜ is the region −1≤*x*≤1, −1≤*y*≤1, *x*^{2}+*y*^{2}≥*R*^{2}. This result follows from integrating equation (2.1) over the volume of tissue occupying the region between *z* and *z*+d*z*, using the divergence theorem to obtain the left-hand side in terms of the fluxes on the boundaries of this region, and then taking the limit d*z*→0. It follows that the zero order concentration in the capillary occupying the *n*th row and *m*th column is governed by the equation(2.28)and the first order concentration in this capillary satisfies(2.29)where *D*^{n,m}=2*D*_{t}*L*/α*U*^{n,m}*R*_{t}*R*_{c}, *N*^{n,m} is the oxygen capacity and *U*^{n,m} is the flow velocity in the capillary occupying the *n*th row and *m*th column. Expanding the left-hand side and substituting the solution obtained for the zero order tissue concentration, equation (2.12), into the integral over ℜ reduces equation (2.29) to(2.30)where(2.31)

Equations (2.23), (2.24) and (2.28) are a system of ordinary differential equations that determine the zero order capillary oxygen concentrations in terms of the zero order fluxes, coupled with a system of algebraic equations for the zero order fluxes in terms of the zero order capillary concentrations. Similarly, with the zero order solutions determined, equations (2.25), (2.26) and (2.30) are an analogous system of ordinary differential equations coupled with algebraic equations for the first order fluxes and the first order capillary concentrations. The initial condition on the capillary concentration, equation (2.7), cannot be imposed on *C*_{0} and *C*_{1} because the expansion, equation (2.8), is not valid in the region near *z*=0. In addition, the no-flux boundary conditions on *c*_{t} at the ends of the tissue region, *z*=0 and *z*=1, equation (2.3), are not satisfied by *c*_{0} and *c*_{1}. A boundary layer solution in these regions must be constructed and matched to the expansion already obtained, which is valid in the region bounded away from *z*=0 and *z*=1. The boundary layer expansion in the region near *z*=0, obtained in Appendix B, determines the initial conditions to be imposed on *C*_{0} and *C*_{1}. This involves the method of matched asymptotic expansions, which has been used in a similar manner for a cylindrical functional unit, or Krogh cylinder, when there is no flux on the outer boundary of the cylinder (Salathe *et al*. 1980). The initial conditions are shown in the appendix to be(2.32)(2.33)where(2.34)and is the initial concentration in the capillary occupying the *n*th row and *m*th column. The constants *f*^{n,m}, *f*^{n,m−1}, *g*^{n−1,m} and *g*^{n,m} are determined from the system of algebraic equations (B 28) and (B 29) in Appendix B.

Equations (2.23) and (2.24) have the form **A x**

_{0}=

*b*_{0}, where

*x*_{0}is a (2

*MN*−

*M*−

*N*) dimensional vector of the unknown zero order fluxes. (There are

*N*×(

*M*−1) fluxes

*F*

^{n,m}and (

*N*−1)×

*M*fluxes

*G*

^{n,m}, since

*F*

^{n,0},

*F*

^{n,M},

*G*

^{0,m}and

*G*

^{N,m}, occurring along the outer boundary of the array, must be specified. In the examples, they are taken to be zero, corresponding to the outer boundary of a closed system.) The matrix

**A**has dimension (2

*MN*−

*M*−

*N*)×(2

*MN*−

*M*−

*N*) and contains only the known constants

*a*

_{0},

*b*

_{0}and

*d*

_{0}, and the vector

*b*_{0}is made up of the

*N*×

*M*capillary concentrations . Specifying all the capillary concentrations at

*z*=0 determines the vector

*b*_{0}, and the fluxes are found from

**A**

^{−1}

*b*_{0}. These fluxes provide the right-hand side of equation (2.28), and the initial concentrations determine the factor , appearing in that equation, so that is determined at

*z*=0. The concentrations at

*z*+Δ

*z*are . These concentrations are used to compute a new vector

*b*_{0}at

*z*+Δ

*z*and then new fluxes from

**A**

^{−1}

*b*_{0}and, from equation (2.28), the derivative , so that the process can be repeated, moving forward in steps of size Δ

*z*. At each stage of this process of moving forward, the first order fluxes and concentrations are found in a similar manner. The system of equations for the first order fluxes, equations (2.25) and (2.26), has the form

**A**

*x*_{1}=

*b*_{1}+

*b*^{*}, where

*x*_{1}and

*b*_{1}are the same as the vectors

*x*_{0}and

*b*_{0}, but with the first order fluxes and concentrations replacing the corresponding zero order quantities, and

*b*^{*}is determined from the zero order solutions. At

*z*=0 the vector

*b*_{1}is formed using the concentrations given by equation (2.33). The first order fluxes are then found from

**A**

^{−1}(

*b*_{1}+

*b*^{*}), which, when substituted into equation (2.30) determines and therefore the first order concentrations at

*z*+Δ

*z*. The first order fluxes and concentrations can therefore be found in the same manner as the zero order quantities by moving forward in steps of size Δ

*z*. However, for the first order solution, the vector

*b*^{*}must be determined at each step. This vector is constructed from the second derivative of the zero order fluxes and concentrations. In addition, the equation for contains the second derivative of the zero order fluxes and concentrations and the quantity . This first derivative of is found at each step as just described, but the second derivatives must also be found before moving to the next step is possible. If the vector

**′**

*b*_{0}is formed using the derivatives in place of the concentrations in

*b*_{0}, then yields the derivatives of the fluxes . For the first step, at

*z*=0, this yields , , which along with the already found, provides the right-hand sides in equations (B 28) and (B 29) for

*f*

^{n,m},

*f*

^{n,m−1},

*g*

^{n−1,m}and

*g*

^{n,m}appearing in equation (2.33). Differentiating equation (2.28) with respect to

*z*gives(2.35)which provides the second derivative of in terms of the quantities already determined. Then, forming

*b*_{0}″ using the second derivatives in place of the concentrations in

*b*_{0}, the second derivative of the fluxes is given by

**A**

^{−1}

*b*_{0}″.

The rate at which oxygen enters at the arterial end minus the rate at which it is consumed within the tissue equals the rate at which oxygen exits at the venous end. In the absence of interaction, this balance holds for each individual functional unit. However, when all the functional units interact, the balance must be taken over the entire vascular bed. The rate at which oxygen enters at the arterial end of all the capillaries is(2.36)and the rate at which it is consumed within all the functional units is(2.37)

Without axial diffusion, the rate at which oxygen is convected out at the venous end through all the capillaries is(2.38)With axial diffusion, the rate at which oxygen is convected out at the venous end is not given by equation (2.38) with replaced by , since the outer expansion is not valid at *z*=1. In order to find the capillary concentration at *z*=1, the boundary layer expansion at the venous end must be obtained. This is obtained in Appendix C. The capillary concentration at *z*=1, denoted , is shown there to be given by(2.39)where(2.40)and(2.41)The constants , , , are determined from the system of algebraic equations (C 25) and (C 26) in Appendix C. Therefore, with axial diffusion, the rate at which oxygen is convected out through all the capillaries at the venous end is(2.42)

Balancing the net rate of oxygen flux into the region with the rate of consumption using the above formulae provides a check on the accuracy of the solutions.

## 3. Results

Results will be illustrated for a 5×5 array of capillaries with different flow velocities but with the same initial concentration at *z*=0 and the same haematocrit, so that and . Since the oxygen solubility in plasma and in tissue are nearly equal (Weerappuli & Popel 1989; Whiteley *et al*. 2002), the ratio *α* will be assumed to be equal to one. The maximum blood flow velocity is *U*^{1,1}, the velocity in the capillary occupying the first row and the first column. Values for *C*_{a}, *U*^{1,1} and the other parameters used are shown in table 1. The flow velocities in all the capillaries relative to *U*^{1,1} are(3.1)

The capillary occupying the second row and second column has the lowest flow velocity, with a value one fifth that of the fastest flowing capillary.

The oxygen concentration as a function of distance along the capillary is shown for the fastest and slowest flowing capillaries in figure 1. The four curves that originate at *z*=0 with the value 1.0 are the solutions without axial diffusion. They are the zero order solutions corresponding to *ϵ*=0, and they satisfy the initial condition equation (2.32), . The top curve of this group corresponds to the fastest flowing capillary and the bottom curve corresponds to the slowest flowing capillary when no interaction occurs among the capillaries. The solutions without interaction are obtained from equations (2.1)–(2.7) with the fluxes *F*_{R}, *F*_{L}, *F*_{T} and *F*_{B} equal to zero. Not enough oxygen enters the slowest flowing capillary at the arterial end to supply the oxygen demand of the tissue in its functional unit, and the concentration falls to zero about half way along its length. The middle two curves of the four show the concentration in these two capillaries when there is interaction. With interaction, the concentration falls more rapidly in the fast flowing capillary, as it supplies oxygen to the more poorly perfused regions, while the concentration falls more slowly in the slow flowing capillary. The concentration in the slow flowing capillary does not fall to zero because some of the tissue surrounding it is supplied by neighbouring capillaries. These results illustrate that interaction produces a much more uniform distribution of oxygen among the capillaries. All the tissue is supplied with oxygen, although a simplified model that does not include interaction incorrectly predicts that some of the tissue becomes anoxic. The four curves that do not start at *z*=0 correspond to the solution with axial diffusion, and represent . At *z*=0 they satisfy the initial condition given in equations (2.32) and (2.33). These solutions are not valid in the neighbourhood of *z*=0 and *z*=1; they represent the oxygen concentrations only in the region bounded away from the ends. Again, the top curve of this group corresponds to the fastest flowing capillary and the bottom curve corresponds to the slowest flowing capillary when no interaction occurs among the capillaries, and the middle two curves show the effect of interaction.

The tissue concentrations are not illustrated here, but the concentration in the tissue does not deviate significantly from the values in the capillaries. The detailed variation in tissue concentration can be found using equations (2.12) and (2.18). However, the location where the capillary concentration falls to zero is a good indication of where anoxia first develops in the tissue. For larger substrates, where there is resistance to movement across the capillary wall, the tissue concentration may be substantially lower than the concentration in the capillary, and the capillary concentration is not a good indication of the tissue concentration. The tissue can become depleted of substrate long before the capillary concentration falls to zero, and this must be taken into account.

With axial diffusion, the concentration initially falls more rapidly at the arterial end because oxygen leaving the capillary there not only supplies the tissue at that axial location, but also diffuses toward the venous end. Throughout the length of the capillary the concentration remains below the levels found by neglecting axial diffusion. However, the concentration decreases more slowly as the venous end is approached because the extra oxygen that left the capillary at the arterial end and diffused axially is available to the tissue at the venous end. As a result, the concentration at the venous end attains the value given in equation (2.39). These concentrations are shown for each of the capillaries in figure 2 without interaction, and in figure 3 with interaction. The values in figure 2 are the same with or without axial diffusion, since without interaction the concentration at *z*=1 can be determined for each capillary from the influx of oxygen at the arterial end and the amount of oxygen consumed within each individual functional unit. It therefore follows from equation (2.39) that when there is no interaction. This can be verified by integrating equation (2.30) from zero to one (with the left-hand side written in the form shown in equation (2.29)), using equation (2.33), setting the fluxes equal to zero, and recalling that , . Figure 3 includes both interaction and axial diffusion. Comparing figures 2 and 3 shows that neglecting interaction incorrectly predicts that several of the capillaries become depleted of oxygen, whereas in fact oxygen is supplied to all of the tissue throughout the entire region.

With interaction, the expected value of the capillary concentration at *z*=1 cannot be determined independently from the balance of input, output and consumption within each functional unit since the oxygen entering a unit at the arterial end does not necessarily stay within that unit. However, the accuracy of the perturbation analysis can be examined by determining the extent to which the total flux of oxygen out at the venous end from all functional units equals the total flux in at the arterial end, diminished by the amount of oxygen consumed in all of the tissue. The accuracy can therefore be defined, using the quantities in equations (2.36) and (2.37) and in either equation (2.38) when axial diffusion is neglected or in equation (2.42) when axial diffusion is included, as(3.2)For the example illustrated, this has magnitude 0.0002 without axial diffusion and 0.0006 with axial diffusion.

The concentrations with axial diffusion closely approach the corresponding curves without axial diffusion in the example illustrated in figure 1, and the venous boundary layer produces only a slightly different concentration at *z*=1 than that given by . A more stringent test of the accuracy of the perturbation analysis is provided by using an array having the same data as the example considered previously, but with capillary lengths of 0.015 cm, or one quarter that given in table 1. As *z*→1 the solution with axial diffusion, , is not close to the corresponding solution without axial diffusion, , and accurately determining the concentrations near *z*=1 requires an analysis of the venous boundary layer. The 5×5 matrices shown below give the concentration for each capillary at *z*=1 when interaction occurs. The first matrix gives the values corresponding to no axial diffusion. The second matrix gives the concentrations , and, as noted, are not the correct concentrations at *z*=1. If equation (2.39) is used to determine the venous concentrations when axial diffusion is included, the correct values shown in the third matrix are obtained. Although the concentrations in the individual capillaries shown in the first and third matrices are different, they yield identical values for total output, while the values in the second matrix yield a total output that is about 5% less than the correct value. Using the values in the first and third matrices to determine the accuracy according to equation (3.2) gives an error of about a tenth of a per cent, while using the values in the second matrix gives an error of about 4%,

## 4. Discussion

An effective method for determining oxygen distribution within an array of interacting capillaries that includes axial diffusion has been presented. The governing equations and their solutions are complicated; however, these solutions, equations (2.12) and (2.18), are needed only to determine the constants shown in equation (2.22). Now that these constants are known, the oxygen concentrations can be found by solving an ordinary differential equation for the concentration in each capillary, equations (2.28) and (2.30), together with a system of algebraic equations, equations (2.23)–(2.26), that couples these differential equations and account for the interaction. Equations (2.28) and (2.30) were solved using a simple first order Euler method of marching forward in steps along the length of the capillary. Other numerical approaches, such as Runge–Kutta and higher order Euler methods, have been used and give comparable results. When there is no interaction, the concentration at the end of the capillary, *z*=1, is the same with or without axial diffusion, and can be determined directly, since for each functional unit the rate of oxygen flowing out at the venous end must equal the rate it enters at the arterial end minus the rate at which it is consumed within the tissue. Therefore, the accuracy of the Euler method used here can be evaluated by the extent to which the concentration attained at the end of the capillary, *z*=1, agrees with the known value there. A step size Δ*z*=0.005 resulted in a very short computational time and yielded results that were accurate to very many decimal places and could not be distinguished from other more accurate integration techniques.

The concentrations all start at the value 1.0, but with axial diffusion a boundary layer analysis is required to obtain the concentration profiles near *z*=0. These boundary layer solutions are not valid in the region bounded away from the arterial end and are not shown in figure 1, as they have not been completely obtained. The boundary layer analysis was carried out only to the extent necessary to obtain the initial condition for the perturbation expansion valid in the region bounded away from the ends. This initial condition gives a concentration profile that begins at *z*=0 but is not valid until some distance from *z*=0. This is emphasized in figure 1 by not starting the curves representing the solutions with axial diffusion at *z*=0. The venous boundary layer matches this solution and provides the profiles from the location where the outer solution ceases to be valid to the end *z*=1. The curves in figure 1 corresponding to axial diffusion have not been extended to *z*=1 because the boundary layer solution at the venous end has not been obtained. The venous boundary layer analysis was carried out only to the extent necessary to obtain the correct concentration at *z*=1. If the solutions for the arterial and venous boundary layer are obtained, they can be combined with the solution in the central region using the techniques of matched asymptotic expansions to construct a solution that is uniformly valid in the region 0≤*z*≤1. For the case of no interaction, such uniformly valid asymptotic expansions were obtained by Salathe *et al*. (1980).

Figures 2 and 3 illustrate the extent to which interaction affects the oxygen concentration in the capillaries and the surrounding tissue. Oxygen diffusion from richly perfused regions into poorly perfused regions can prevent the development of anoxia. This effect is further enhanced by the presence of myoglobin in skeletal muscle, which facilitates the diffusion of oxygen at low concentration. Myoglobin facilitated diffusion in a multicapillary system has been studied previously using the homogenization approach (Salathe & Chen 1993). Myoglobin was shown to be effective in retarding the development of anoxia following occlusion, by enhancing large-scale interaction among the capillaries. Hoofd (1995*a*,*b*) studied the heterogeneous distribution of the capillaries within the cross-section, and compared the resulting concentration profiles to those resulting from a homogeneous distribution of capillaries. Heterogeneous distribution of capillaries can be included in the present model by placing the capillaries off the centre line of the box, and the extent to which such heterogeneous distribution affects large-scale interaction within the capillary bed could be examined. Axial diffusion results in a more rapid initial decrease in capillary concentration, but without interaction the concentration returns to the same value at the venous end with or without axial diffusion. However, the altered concentrations owing to axial diffusion affect the extent of interaction among the capillaries, and therefore the concentrations at the venous end. Although the solutions here are not valid near the arterial end, they are able to describe the concentrations beyond this region and to indicate the extent to which large-scale oxygen diffusion may prevent the formation of anoxia.

## Appendix

Equations (2.15)–(2.17) for *c*_{1} are solved by writing(A1)where *c*_{P} provides the right-hand side of equation (2.15) and *c*_{H} provides the boundary conditions, equation (2.16), for that equation. Therefore, *c*_{H} satisfies the equation(A2)and the boundary conditions(A3)with(A4)

By analogy with the zero order solution, it follows that(A5)

The portion *c*_{p} satisfies(A6)

(A7)

(A8)

Writing *c*_{p} as(A9)then *c*_{R}, *c*_{L}, *c*_{T} and *c*_{B} each satisfy Laplace's equation and the boundary condition(A10)Zero boundary conditions are imposed for *c*_{R}, *c*_{L}, *c*_{T} and *c*_{B} along the three edges for which a flux is not specified in equation (A 10). The function *h*(ξ) was defined previously (equation (2.19)) and satisfies the condition(A11)as does the function cos *mπξ*, also occurring in the boundary conditions, equation (A 10). Therefore, *c*_{R}, *c*_{L}, *c*_{T} and *c*_{B} each satisfy a well posed boundary value problem and can be solved by standard methods using Fourier series. They are made unique by requiring each to vanish at *x*=*y*=0, therefore satisfying the boundary condition at the capillary wall, equation (A 8), only approximately. The result yields the solution given by equation (2.18).

## Appendix

The boundary layer at the arterial end is examined by introducing the boundary layer variable *Z*=*z*/*ϵ*, defining , , , , , , and rewriting the governing equations in terms of these quantities:(B1)

(B2)

(B3)

(B4)

These equations must be augmented with a suitable matching condition for the boundary layer solution as *Z*→∞. To determine this matching condition, the earlier solution is rewritten in terms of the boundary layer variable *Z* and expanded as a series in *ϵ*. From equation (2.27), the solution for the zero order capillary concentration is(B5)

From this it follows that the capillary concentration when expressed as a function of *Z* has an expansion as a series in *ϵ* of the form(B6)where(B7)

Expanding the tissue solution and the fluxes as a series in *ϵ* to order *ϵ* gives(B8)

These expansions of the outer solution show that the boundary layer expansion is of the form(B9)

Substituting these expansions into the boundary layer equations and boundary conditions, equations (B 1)–(B 4), yields for *Ψ*,(B10)(B11)(B12)and the matching conditions that as *Z*→∞,(B13)(B14)and for *P*(*Z* ) the equation(B15)with the initial condition(B16)and the matching condition(B17)

From equations (B 15) and (B 16) it follows that(B18)which, together with the matching condition, equation (B 17), and the definition of *ϕ* and *ψ*, yields(B19)

Integrating equation (B 10) over the region of tissue between *Z* and *Z*+d*Z*, and using the divergence theorem, shows that the sum of all the fluxes out of this region equals zero. Taking the limit d*Z*→0 yields a relationship that reduces equation (B 19) to(B20)which, in view of the matching condition, equation (B 13), gives(B21)

The fluxes *H*_{R}, *H*_{L}, *H*_{T} and *H*_{B} are found by solving equations (B 10)–(B 13) for *Ψ* and matching the average concentrations along the boundaries of adjacent units at each value of *Z*. However, the integrals of the functions appearing in equation (B 21), and therefore the initial concentration, *C*_{1}(0), can be determined without solving for *Ψ*(*x*,*y*,*Z*). Writing(B22)and letting *w*(*x*,*y*,*s*) be the Fourier cosine transform of *W*,(B23)*w* satisfies the equation(B24)and the boundary conditions(B25)

The function *h*_{R}(*s*) is the Fourier cosine transform of ,(B26)with corresponding definitions for *h*_{L}(*s*), *h*_{T}(*s*) and *h*_{B}(*s*). Solving for *w*(*x*,*y*,*Z*) and matching average values at the boundaries of adjacent units yields a system of algebraic equations from which *h*_{R}(*s*), *h*_{L}(*s*), *h*_{T}(*s*) and *h*_{B}(*s*) can be determined. The inverse transforms yield the functions appearing in the integrals in equation (B 21) for *C*_{1}(0). However, these integrals are , , and , and so it is necessary to solve equations (B 24) and (B 25) only for *s*=0. When *s*=0, they are a system of equations for that are identical to equations (2.15)–(2.17) for *c*_{1}(*x*,*y*,*z*), with *C*_{1}(*z*)=0, with *C*″_{0}(*z*), *A*″_{0}(*z*), *B*″_{0}(*z*), *Q*″_{0}(*z*) and *G*″_{0}(*z*) replaced by *C*′_{0}(0), *A*′_{0}(0), *B*′_{0}(0), *Q*′_{0}(0) and *G*′_{0}(0), respectively (recall *φ*=*C*′_{0}(0)), and with *F*_{R1}(*z*), *F*_{L1}(*z*), *F*_{T1}(*z*) and *F*_{B1}(*z*) replaced with , , and . The average values of along each of the four boundaries of the functional unit are therefore those given for *c*_{1} in equation (2.21) with the substitutions described above. For the unit in the *n*th row and *m*th column, let , , and . Then the initial condition on the first order concentration for the capillary in this unit is(B27)where the *f*^{n,m}, *g*^{n,m} are found from(B28)and(B29)which is a system of algebraic equations for *f*^{n,m}, *g*^{n,m} analogous to equations (2.25) and (2.26) for *F*^{n,m}, *G*^{n,m}.

## Appendix

The boundary layer at the venous end, *z*=1, is treated in a similar manner to the boundary layer at the arterial end, using the boundary layer variable . Defining , , , , , , and rewriting the governing equations in terms of these quantities, yields(C1)

(C2)

(C3)

(C4)

No value is specified for at analogous to the condition in equation (B 4), since *C*_{V}, the capillary concentration at *z*=1, is not known, but is to be determined. These equations must be augmented with a suitable matching condition for the boundary layer solution as . To determine this matching condition, the earlier solution is rewritten in terms of the boundary layer variable and expanded as a series in *ϵ*. From equation (2.27), the solution for the zero order capillary concentration is(C5)

From this it follows that the capillary concentration when expressed as a function of has an expansion as a series in *ϵ* of the form(C6)where(C7)

Expanding the tissue solution and the fluxes to order *ϵ* gives(C8)

These expansions of the outer solution show that the venous boundary layer expansion is of the form(C9)

(C10)

(C11)

Substituting these expansions into the boundary layer equations and boundary conditions yields, for *Φ*,(C12)

(C13)

(C14)

with the matching conditions that as (C15)(C16)and for *Γ*(*Z*) the equation(C17)with the matching condition(C18)

From equation (C 17) it follows that(C19)which, together with the matching condition, equation (C 18), and the definitions of *ϕ*_{V} and *ψ*_{V}, yields(C20)

Integrating equation (C 1) over the region of tissue between and , and using the divergence theorem, shows that the sum of all the fluxes out of this region equals zero. Taking the limit yields a relationship that reduces equation (C 20) to(C21)which, in view of the matching condition, equation (C 15), gives(C22)

The integrals appearing in equation (C 22) can be determined in the same manner as for the arterial boundary layer. Let denote the Fourier cosine transform of , with corresponding definitions for , and . For the unit in the *n*th row and *m*th column, let , , and . Then the concentration at the venous end for the capillary in this unit is, from equation (C 9),(C23)where is given by(C24)

Then and are determined from the algebraic equations(C25)and(C26)

## Footnotes

- Received May 20, 2003.
- Accepted July 19, 2004.

- © 2005 The Royal Society