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

**|**HHS Author Manuscripts**|**PMC5609729

Formats

Article sections

- Abstract
- 1. INTRODUCTION
- 2. MATHEMATICAL STATEMENT OF THE PROBLEM
- 3. SOLUTION OF THE DIFFUSION EQUATION
- 4. EXAMPLES
- 5. CONCLUDING REMARKS
- References

Authors

Related links

Math Biosci. Author manuscript; available in PMC 2017 September 22.

Published in final edited form as:

PMCID: PMC5609729

NIHMSID: NIHMS857492

ALEKSANDER S. POPEL, Department of Chemical Engineering, University of Arizona, Tucson, Arizona 85721;

The publisher's final edited version of this article is available at Math Biosci

The diffusive transport of a substance between a parallel capillary network and the surrounding tissue is investigated. The consumption/production rate of the substance in the tissue is assumed to be constant (zero-order chemical kinetics). The solution of the diffusion problem which describes the distribution of the substance in the tissue and along the capillary network is found in an analytical form. A rather general assumption regarding the symmetry of capillary network makes it possible to formulate a Neumann-type boundary-value problem in a rectangular domain. The solution of the diffusion problem in the rectangle allows the capillary-tissue fluxes to be expressed linearly in terms of the concentrations in the capillaries, and hence leads to ordinary differential equations for those concentrations. Several examples are considered with different network geometry and concurrent or countercurrent flow conditions. The solution makes it possible to investigate the effect of capillary interaction on mass transfer in various microcirculatory units.

The diffusive transport of various substances between capillary blood vessels and surrounding tissue structures has been studied extensively by numerous investigators both theoretically and experimentally. The theoretical work in this area has been recently reviewed by Middleman [1] and Leonard and Jørgensen [2].

The first simple mathematical model of the molecular transport between blood capillaries and surrounding tissue cylinders, based on the assumption that the diffusion flux vanishes at the external boundary of the tissue cylinder, was proposed in Krogh’s paper [3]. This model, due to its simplicity, is widely employed in physiological studies concerning the estimates of the supply conditions in tissues.

Krogh’s model has been modified by several authors. Thews [4], for example, took into account the axial diffusion in the tissue cylinder, and Reneau et al. [5] solved the transport equations for the tissue cylinder numerically, taking into account the intracapillary diffusion and the axial diffusion in the tissue region.

Using certain simplifying assumptions, Hyman [6] and Stewart and Morazzi [7] obtained approximate solutions within the framework of the tissue cylinder model.

In order to describe capillary interactions, Wilson [8] proposed a model of the tissue cylinder with a boundary condition on the external surface of the cylinder which differed from Krogh’s. Wilson set the concentration on this surface equal to a constant though unspecified value, and the integral of the diffusion flux through the surface equal to zero. This model leads to a concentration distribution in the tissue being significantly different from that predicted by Krogh’s model; however, it is difficult to interpret Wilson’s boundary condition in terms of the capillary network geometry.

Diemer [9], Bailey [10], and Reneau and Knisely [11] demonstrated different approaches to the problem of the mass transfer for countercurrent capillary-flow conditions. It was shown that countercurrent flow provides more homogeneous supply to the tissue than concurrent flow does.

Metzger [12, 13], Grunewald [14], and Grunewald and Sowa [15] considered the diffusion in tissues perfused by two- or three-dimensional capillary networks. The diffusion equation was solved numerically by the finite-difference technique. Results of the computer simulation of the oxygen transport in skeletal muscle appeared in a recent paper [15].

In the present work we will derive an analytical solution to the problem of the diffusive transport in tissues for a three-dimensional capillary bed. Certain simplifying assumptions make it possible to present the solution in a closed form and to investigate a number of examples that are of importance for physiological applications.

We consider the problem of molecular transport in a tissue supplied by an infinite array of parallel capillaries [Fig. 1(a)]. A rectangular coordinate system (*x*,*y*,*s*) is introduced such that the *s*-axis is parallel to the capillaries. The capillary distribution in the *x–y* plane is illustrated in Fig. 1(b).

The steady-state diffusion in the tissue is governed by the equation

$$D\phantom{\rule{0.16667em}{0ex}}\left(\frac{{\partial}^{2}c}{\partial {x}^{2}}+\frac{{\partial}^{2}c}{\partial {y}^{2}}\right)+{D}_{s}\frac{{\partial}^{2}c}{\partial {s}^{2}}-g=0,$$

(1)

where *c* is the local concentration of the substance, *D* is the diffusion coefficient in a plane perpendicular to the capillaries, *D _{s}* is the diffusion coefficient in the direction along the capillaries, and

Assuming that the axial diffusion of the substance in the blood can be neglected and that the concentration is uniform in the cross-section of the capillary (cf. Aroesty and Gross [16]), the convective transport in the capillaries can be described by the equations

$${Q}_{j}\frac{d}{ds}\left[{\overline{c}}_{j}+\chi \left({\overline{c}}_{j}\right)\right]=-{J}_{{j}^{\prime}},$$

(2)

where the index *j* identifies the capillary, *Q _{j}* =

In order to find an analytical solution to the problem, the following assumptions are made:

- Diffusion in the
*s*-direction in the tissue can be neglected. This can be justified if either*D*= 0 or_{s}*r***L***,*where*r** is the characteristic intercapillary distance and*L** is the characteristic capillary length (cf. [15]). - The consumption/production rate is constant,
*g*= const. - The capillary radii are small in comparison with the intercapillary distances,
*r*_{j}*r**. - A certain symmetry of the capillary bed will be assumed; in order to describe it in mathematical terms we introduce an appropriate notation.

The location of the capillary blood vessel can be given by the relationships

$$\left|{(x-{x}_{j})}^{2}+{(y-{y}_{j})}^{2}\right|\le {r}_{j},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{s}_{j}\le s\le {s}_{j}+{L}_{j},$$

(3)

where *x _{j}*,

$$({x}_{j},{y}_{j},{s}_{j},{L}_{j},{Q}_{j},{c}_{j0}).$$

(4)

With this notation in mind we can describe the symmetry of the capillary bed required for the solution in the following form: a rectangular coordinate system (*x*,*y*,*s*) with the *s*-direction parallel to the capillaries, and constants *a* and *b* can be chosen so that for any capillary in the rectangular cylinder Ω (0 ≤ *x* ≤ 2*a*, 0 ≤ *y* ≤ 2*b,* −∞ < *s* < ∞) described by the set (4) there is an infinite number of capillaries outside of Ω, characterized by the sets

$$(\pm {x}_{j}+4am,\pm {y}_{j}+4bn,{s}_{j},{L}_{j},{Q}_{j},{c}_{j0}),$$

(5)

where *m* and *n* are integers, and all possible combinations of plus and minus signs in (5) should be considered.

The symmetry condition (5) describes an infinite double periodic capillary array. It enables us to reduce the two-dimensional problem in the infinite plane (*x*, *y*), governed by Eq. (1) with the *s*-derivatives neglected, to the corresponding problem in the fundamental rectangle *ω* (cross-section of the rectangular cylinder Ω).

Taking into consideration the symmetry conditions (5), one can infer that the concentration possesses similar symmetric properties, namely,

$$\begin{array}{l}c(x,y,s)=c(x,-y,s)=c(-x,y,s),=c(-x,-y,s),\\ c(x,y,s)=c(x+4am,y+4bn,s).\end{array}$$

(6)

The symmetry (6) permits us to formulate the no-flux boundary condition on the boundary Γ of the fundamental rectangle in the form

$$\frac{\partial c}{\partial n}=0.$$

(7)

On the capillary-tissue interface Γ* _{j}*,

$${J}_{j}=-2\phantom{\rule{0.16667em}{0ex}}\pi \phantom{\rule{0.16667em}{0ex}}{r}_{j}D{\int}_{{\mathrm{\Gamma}}_{j}}\frac{\partial c}{\partial n}d{\mathrm{\Gamma}}_{j},$$

(8)

$$c=\nu {\overline{c}}_{j}.$$

(9)

Here *n* is the unit vector normal to Γ* _{j}* pointing in the external direction, and

The equations (2) will be solved using the boundary conditions

$${\overline{c}}_{j}({s}_{j})={c}_{j0}\text{if}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{Q}_{j}>0,$$

(10)

$${\overline{c}}_{j}({s}_{j}+{L}_{j})={c}_{j0}\text{if}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{Q}_{j}<0,$$

(11)

where *c _{j}*

Since the capillary radii are small in comparison with the intercapillary distances, we can further simplify the problem by replacing the capillaries by line sources (or sinks) with the intensity *J _{i}* per unit length. Therefore, Eq. (1) can be rewritten in the form

$$\frac{{\partial}^{2}c}{\partial {x}^{2}}+\frac{{\partial}^{2}c}{\partial {y}^{2}}-G+\sum _{j}{q}_{j}\delta \phantom{\rule{0.16667em}{0ex}}(x-{x}_{j},y-{y}_{j})=0,$$

(12)

where

$$G=\frac{g}{D},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{q}_{j}=\frac{{J}_{j}}{D},$$

(13)

and *δ*(*x*,*y*) is the two-dimensional Dirac delta function.

Generally speaking, replacing the capillaries by the line sources or sinks makes it impossible to satisfy the boundary condition (9), since the latter requires that *c* = const on Γ* _{j}*. We will approximately satisfy (9) by equating

$$\nu {\overline{c}}_{j}=\frac{1}{\mid {\stackrel{\sim}{\mathrm{\Gamma}}}_{j}\mid}{\int}_{{\stackrel{\sim}{\mathrm{\Gamma}}}_{j}}c\phantom{\rule{0.16667em}{0ex}}d{\stackrel{\sim}{\mathrm{\Gamma}}}_{j},$$

(14)

where Γ̃* _{j}* denotes the portion of the capillary-tissue interface confined within

Let us introduce a new coordinate system so that the fundamental rectangle ω is determined by the inequalities

$$\mid x\mid \phantom{\rule{0.16667em}{0ex}}\le a,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}0\le y\le 2b.$$

(15)

Substitution of

$$\frac{{\partial}^{2}u}{\partial {x}^{2}}+\frac{{\partial}^{2}u}{\partial {y}^{2}}=-\sum _{j}{q}_{j}\delta \phantom{\rule{0.16667em}{0ex}}(x-{x}_{j},y-{y}_{j})$$

(17)

with the Neumann-type boundary conditions

$$\frac{\partial u}{\partial n}=0\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}y=0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}y=2b,$$

(18)

$$\frac{\partial u}{\partial n}=Ga\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}x=\pm a,$$

(19)

where *n* is the internal unit normal to the boundary Γ.

Let *z* = *x* + *iy* be a complex variable, and *N*(*z′*,*z*) be the Neumann function of the problem (17)–(19). The solution *u*(*z*) can be written in the form [17]

$$u(z)=-{\int}_{\mathrm{\Gamma}}\frac{\partial u({z}^{\prime})}{\partial n}N\phantom{\rule{0.16667em}{0ex}}({z}^{\prime},z)d\mathrm{\Gamma}+{\int}_{\omega}F({z}^{\prime})N\phantom{\rule{0.16667em}{0ex}}({z}^{\prime},z)d{x}^{\prime}d{y}^{\prime}+{c}^{\ast},$$

(20)

where *F*(*z*) = Σ*q _{j}δ*(

The necessary condition for the resolution of the Neumann problem

$$\sum _{j}{q}_{j}^{\ast}=4\mathit{abG},$$

(21)

where

$${q}_{j}^{\ast}=\{\begin{array}{lll}{q}_{j}\hfill & \text{if}\hfill & {z}_{j}\in \omega ,\phantom{\rule{0.16667em}{0ex}}{z}_{j}\notin \mathrm{\Gamma},\hfill \\ \frac{{q}_{j}}{2}\hfill & \text{if}\hfill & {z}_{j}\in \mathrm{\Gamma},\phantom{\rule{0.16667em}{0ex}}{z}_{j}\ne \pm a,\pm a+2bi,\hfill \\ \frac{{q}_{j}}{4}\hfill & \text{if}\hfill & {z}_{j}=\pm a,\pm a+2bi.\hfill \end{array}$$

(22)

The relationship (21) has a simple physical interpretation: it expresses the conservation of mass of the substance in the rectangle (15).

In order to find the Neumann function *N*(*z′*,*z*), we build up a transformation which conformally maps the rectangle (15) onto the upper half plane Im *ζ* ≥ 0, where

$$\zeta =\xi +i\eta .$$

(23)

This transformation is effected by the Jacobian elliptic function [18]

$$\zeta =\text{sn}\left(\frac{zK}{a},k\right),$$

(24)

where *K*(*k*) is the complete elliptic integral of the first kind,

$$K(k)={\int}_{0}^{1}\frac{dt}{\sqrt{(1-{t}^{2})(1-{k}^{2}{t}^{2})}},$$

(25)

and the parameter *k* is a root of the transcendental equation

$$\frac{K(\sqrt{1-{k}^{2}})}{K(k)}=\frac{2b}{a},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}0<k<1.$$

(26)

In the following equations the standard notation used will be

$${k}^{\prime}=\sqrt{1-{k}^{2}},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{K}^{\prime}(k)=K({k}^{\prime}).$$

(27)

When the functions *K* or *K′* are written without specifying the argument, then the argument *k* is implied. The same agreement is also valid for the Jacobian elliptic functions, namely, sn *u* = sn(*u*, *k*).

The inverse transformation mapping the upper half plane Im *ζ* ≥ 0 onto the rectangle (15) in the *z*-plane is effected by the incomplete elliptic integral of the first kind,

$$z=\frac{a}{K}{\int}_{0}^{\zeta}\frac{dt}{\sqrt{(1-{t}^{2})(1-{k}^{2}{t}^{2})}}.$$

(28)

The analytic function (24) transforms the boundary Γ of the rectangle *ω* into the straight line Im *ζ* = 0 so that the points

$$\begin{array}{l}z=0\leftrightarrow \zeta =0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}z=\pm a\leftrightarrow \zeta =\pm 1,\\ z=\pm a+2bi\leftrightarrow \zeta =\pm \frac{1}{k},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}z=2bi\leftrightarrow \zeta =\infty \end{array}$$

(29)

correspond to each other.

Introducing the notation

$$Z=X+iY=\frac{zK}{a},$$

(30)

and using the known properties of the Jacobian elliptic functions sn(*z*, *k*), cn(*z*, *k*), and dn(*z*, *k*) [18], one can write (24) in terms of the real variables

$$\begin{array}{l}\xi =\frac{\text{sn}(X,k)\text{dn}(Y,{k}^{\prime})}{{\text{cn}}^{2}(Y,{k}^{\prime})+{k}^{2}{\text{sn}}^{2}(X,k){\text{sn}}^{2}(Y,{k}^{\prime})},\\ \eta =\frac{\text{cn}(X,k)\text{dn}(X,k)\text{sn}(Y,{k}^{\prime})\text{cn}(Y,{k}^{\prime})}{{\text{cn}}^{2}(Y,{k}^{\prime})+{k}^{2}{\text{sn}}^{2}(X,k){\text{sn}}^{2}(Y,{k}^{\prime})}.\end{array}$$

(31)

After transformation to the new variables (*ξ*, *η*), (Eq. (17) takes the form

$$\frac{{\partial}^{2}u}{\partial {\xi}^{2}}+\frac{{\partial}^{2}u}{\partial {\eta}^{2}}=-\sum _{j}{q}_{j}\delta \phantom{\rule{0.16667em}{0ex}}(\zeta -{\zeta}_{j}),$$

(32)

where *ζ _{j}* = sn(

$$\delta \phantom{\rule{0.16667em}{0ex}}(\zeta -{\zeta}_{j})={\left|\frac{dz}{d\zeta}\right|}^{2}\delta \phantom{\rule{0.16667em}{0ex}}(z-{z}_{j})$$

(33)

has been used.

Denoting the image of the boundary Γ as Γ* and taking into account the relationships

$${\left(\frac{\partial u}{\partial n}\right)}_{{\mathrm{\Gamma}}^{\ast}}={\left(\frac{\partial u}{\partial n}\right)}_{\mathrm{\Gamma}}{\mid \frac{dz}{d\zeta}\mid}_{{\mathrm{\Gamma}}^{\ast}},$$

(34)

$$d\mathrm{\Gamma}={\mid \frac{dz}{d\zeta}\mid}_{{\mathrm{\Gamma}}^{\ast}}d{\mathrm{\Gamma}}^{\ast},$$

(35)

$$dx\phantom{\rule{0.16667em}{0ex}}dy={\mid \frac{dz}{d\zeta}\mid}^{2}d\xi d\eta ,$$

(36)

we can rewrite the solution (20) in the form

$$u(\zeta )=-{\int}_{{\mathrm{\Gamma}}^{\ast}}{\left(\frac{\partial u}{\partial n}\right)}_{{\mathrm{\Gamma}}^{\ast}}{N}^{\ast}({\zeta}^{\prime},\zeta )d{\mathrm{\Gamma}}^{\ast}+{\int}_{{\omega}^{\ast}}{F}^{\ast}({\zeta}^{\prime}){N}^{\ast}({\zeta}^{\prime},\zeta )d\xi d\eta ,$$

(37)

where *ω** denotes the domain Im *ζ* ≥ 0, and

$${F}^{\ast}(\xi )=-\sum _{j}{q}_{j}\delta \phantom{\rule{0.16667em}{0ex}}(\zeta -{\zeta}_{j}),$$

(38)

$$N\phantom{\rule{0.16667em}{0ex}}({z}^{\prime},z)={N}^{\ast}\left(\text{sn}\frac{{z}^{\prime}K}{a},\text{sn}\frac{zK}{a}\right).$$

(39)

The Neumann function for the upper half plane Im ζ ≥ 0 can be constructed using the method of images [17]:

$${N}^{\ast}({\zeta}^{\prime},\zeta )=-\frac{1}{2\pi}ln\mid {\zeta}^{\prime}-\zeta \mid \mid {\zeta}^{\prime}-\overline{\zeta}\mid ,$$

(40)

where is the complex conjugate of *ζ*.

Calculation of the distortion coefficient of the conformal mapping (28) yields

$${\left|\frac{dz}{d\zeta}\right|}_{{\mathrm{\Gamma}}^{\ast}}=\{\begin{array}{lll}\frac{a}{K}\frac{1}{\sqrt{({\xi}^{2}-1)(1-{k}^{2}{\xi}^{2})}}\hfill & \text{as}\hfill & 1<\phantom{\rule{0.16667em}{0ex}}\mid \xi \mid \phantom{\rule{0.16667em}{0ex}}<\frac{1}{k},\hfill \\ \frac{a}{K}\frac{1}{\sqrt{(1-{\xi}^{2})(1-{k}^{2}{\xi}^{2})}}\hfill & \text{as}\hfill & \mid \xi \mid \phantom{\rule{0.16667em}{0ex}}<1\phantom{\rule{0.16667em}{0ex}}\text{or}\phantom{\rule{0.16667em}{0ex}}\mid \xi \mid \phantom{\rule{0.16667em}{0ex}}>\frac{1}{k}.\hfill \end{array}$$

(41)

In view of (37)–(41) and (34)–(36), after appropriate transformations, the solution can be written in the form

$$u(\zeta )={u}_{1}(\zeta )+{u}_{2}(\zeta )+{c}^{\ast},$$

(42)

where

$${u}_{1}=\frac{G{a}^{2}}{2\pi K}{\int}_{0}^{1}\frac{ln\left\{\left[{\left(\frac{\sqrt{1-{{k}^{\prime}}^{2}{t}^{2}}}{k}-\xi \right)}^{2}+{\eta}^{2}\right]\phantom{\rule{0.16667em}{0ex}}\left[{\left(\frac{\sqrt{1-{{k}^{\prime}}^{2}{t}^{2}}}{k}+\xi \right)}^{2}+{\eta}^{2}\right]\right\}}{\sqrt{(1-{t}^{2})(1-{{k}^{\prime}}^{2}{t}^{2})}}dt,$$

(43)

$${u}_{2}=-\sum _{j}\frac{{q}_{j}^{\ast}}{2\pi}ln\left[\sqrt{{(\xi -{\xi}_{j})}^{2}+{(\eta -{\eta}_{j})}^{2}}\sqrt{{(\xi -{\xi}_{j})}^{2}+{(\eta +{\eta}_{j})}^{2}}\right].$$

(44)

Therefore the concentration *c* is given by

$$c=\frac{G{x}^{2}}{2}+{u}_{1}(\xi ,\eta )+{u}_{2}(\xi ,\eta )+{c}^{\ast},$$

(45)

where *u*_{1} and *u*_{2} are expressed by (43) and (44), and *ξ*(*x*,*y*) and *η*(*x*,*y*) are determined by (31). The parameters *q _{j}* and

It can be shown that

$${c}^{\ast}=limc\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{as}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}z\to 2bi.$$

(46)

In order to verify (46) we notice that as *z*→2*bi*, *ζ* → ∞ and the functions *u*_{1} and *u*_{2} can be expanded as

$${u}_{1}=\frac{G{a}^{2}}{\pi}\frac{{K}^{\prime}}{K}ln({\xi}^{2}+{\eta}^{2})+O\phantom{\rule{0.16667em}{0ex}}\left(\frac{1}{{\xi}^{2}+{\eta}^{2}}\right),$$

(47)

$${u}_{2}=\sum _{j}\frac{{q}_{j}^{\ast}}{2\pi}ln({\xi}^{2}+{\eta}^{2})+O\phantom{\rule{0.16667em}{0ex}}\left(\frac{1}{{\xi}^{2}+{\eta}^{2}}\right).$$

(48)

The condition of conservation of mass (21) causes the coefficient at ln(*ξ*^{2} + *η*^{2}) in (45) to vanish. Hence, the sum of the first three terms in (45) goes to zero as *z*→2*bi*, which proves the statement (46).

The integral (43) can be expressed in a simpler form at certain points in the fundamental rectangle marked in Fig. 2. Introducing the function *υ*(*x*,*y*) as

$$\frac{G{a}^{2}}{\pi}\frac{{K}^{\prime}}{K}\upsilon =\frac{G{x}^{2}}{2}+{u}_{1},$$

(49)

and using the tables of elliptic integrals [19], we obtain after some simple though lengthy calculations:

$$z=0,\zeta =0\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =ln\frac{1}{k},$$

(50)

$$z=\pm a,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm 1:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =ln\frac{{k}^{\prime}}{k},$$

(51)

$$z=\pm a+2bi,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm \frac{1}{k}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =ln\frac{{k}^{\prime}}{{k}^{2}},$$

(52)

$$z=bi,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\frac{i}{\sqrt{k}}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =\frac{1}{2}ln\frac{2(1+k)}{{k}^{2}\sqrt{k}},$$

(53)

$$z=\pm a+bi,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm \frac{1}{\sqrt{k}}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =\frac{1}{2}ln\frac{2(1-k)}{{k}^{2}\sqrt{k}},$$

(54)

$$z=\pm \frac{a}{2},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm \frac{1}{\sqrt{1+{k}^{\prime}}}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =\frac{1}{2}ln\frac{2(1-{k}^{\prime}){k}^{\prime}\sqrt{{k}^{\prime}}}{{k}^{4}},$$

(55)

$$z=\pm \frac{a}{2}+bi,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm \frac{\sqrt{1+k}+i\sqrt{1-k}}{\sqrt{2k}}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =\frac{1}{2}ln\frac{2(1-k)}{{k}^{2}\sqrt{k}},$$

(56)

$$z=\pm \frac{a}{2}+2bi,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\zeta =\pm \frac{1}{\sqrt{1-{k}^{\prime}}}:\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\upsilon =\frac{1}{2}ln\frac{2(1+{k}^{\prime}){k}^{\prime}\sqrt{{k}^{\prime}}}{{k}^{4}}.$$

(57)

In order to satisfy the boundary condition (14) we have to integrate the concentration *c* expressed by (45) over the capillary-tissue interface Γ* _{j}*. However, this is very difficult to do because of the complexity of relationships (31).

In order to avoid this computation we will again use the fact that the capillaries are narrow in comparison with the intercapillary distances and, therefore, with the dimensions of the rectangle *ω*. Since a circle *z* = *z _{j}* +

$$\zeta =\text{sn}\left({\zeta}_{j}+{r}_{j}\frac{K}{a}{e}^{i\theta},k\right),$$

(58)

we can use the condition *r _{j}K*/

$$\zeta ={\zeta}_{j}+\frac{{r}_{j}K}{a}\left|\text{cn}\frac{{z}_{j}K}{a}\text{dn}\frac{{z}_{j}K}{a}\right|{e}^{i(\theta -{\theta}_{j}^{\prime})}+\frac{1}{2}{\left(\frac{{r}_{j}K}{a}\right)}^{2}\left|\text{sn}\frac{{z}_{j}K}{a}\left({\text{dn}}^{2}\frac{{z}_{j}K}{a}+{k}^{2}{\text{cn}}^{2}\frac{{z}_{j}K}{a}\right)\right|{e}^{2i(\theta -{\theta}_{j}^{\prime \prime})}+\cdots ,$$

(59)

where ${\theta}_{j}^{\prime}$ and ${\theta}_{j}^{\prime \prime}$ are constants.

We consider three different cases:

*z*_{j}*ω*,*z*Γ. The relationship (59) shows that the image of the circle Γ_{j}is to the order_{j}*o*(*r*/_{j}*a*) a circle ${\mathrm{\Gamma}}_{j}^{\ast}$:$$\mid \zeta -{\zeta}_{j}\mid ={R}_{j},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{R}_{j}=\frac{{r}_{j}K}{a}\left|\text{cn}\frac{{z}_{j}K}{a}\text{dn}\frac{{z}_{j}K}{a}\right|.$$(60)Substitution of (60) into (45) gives to the lowest order with respect to the small parameter*r*/_{j}*a*$$\frac{1}{2\pi {r}_{j}}{\int}_{{\mathrm{\Gamma}}_{j}}c\phantom{\rule{0.16667em}{0ex}}d{\mathrm{\Gamma}}_{j}=\frac{G{x}_{j}^{2}}{2}+{u}_{1}({\zeta}_{j})-\frac{{q}_{j}}{2\pi}ln(2{R}_{j}{\eta}_{j})-\underset{l\ne j}{{\sum}^{\prime}}\frac{{q}_{l}^{\ast}}{2\pi}ln{\rho}_{jl}{\overline{\rho}}_{jl}+{c}^{\ast},$$(61)where$$\begin{array}{l}{\rho}_{jl}=\sqrt{{({\xi}_{j}-{\xi}_{l})}^{2}+{({\eta}_{j}-{\eta}_{l})}^{2}},\\ {\overline{\rho}}_{jl}=\sqrt{{({\xi}_{j}-{\xi}_{l})}^{2}+{({\eta}_{j}+{\eta}_{l})}^{2}}.\end{array}$$(62)*z*Γ,_{j}*z*≠ ±_{j}*a*, ±*a*+ 2*bi*. In this case, the semi-circle Γis mapped approximately into a semi-circle ${\mathrm{\Gamma}}_{j}^{\ast}$ in the_{j}*ζ*-plane with the radius*R*determined by (60). The integral in the boundary condition (14) is calculated as_{j}$$\frac{1}{\pi {r}_{j}}{\int}_{{\mathrm{\Gamma}}_{j}}c\phantom{\rule{0.16667em}{0ex}}d{\mathrm{\Gamma}}_{j}=\frac{G{x}_{j}^{2}}{2}+{u}_{1}({\zeta}_{j})-\frac{{q}_{j}}{2\pi}ln{R}_{j}-\underset{l\ne j}{{\sum}^{\prime}}\frac{{q}_{l}^{\ast}}{\pi}ln{\rho}_{jl}+{c}^{\ast}.$$(63)*z*= ±_{j}*a*, ±*a*+ 2*bi*. In this case, the quarter of the circle Γis mapped approximately into semi-circle ${\mathrm{\Gamma}}_{j}^{\ast}$ in the_{j}*ζ*-plane:$$\mid z\pm a\mid ={r}_{j}\to \mid \zeta \pm 1\mid ={R}_{j},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{R}_{j}={\left(\frac{{r}_{j}K{k}^{\prime}}{a\sqrt{2}}\right)}^{2},$$(64)$$\mid z\pm a-2bi\mid ={r}_{j}\to \left|\zeta \pm \frac{1}{k}\right|={R}_{j},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{R}_{j}={\left(\frac{{r}_{j}K{k}^{\prime}}{a\sqrt{2k}}\right)}^{2}.$$(65)The average concentration can be calculated using the relationships$$\frac{2}{\pi {r}_{j}}{\int}_{{\mathrm{\Gamma}}_{j}}c\phantom{\rule{0.16667em}{0ex}}d{\mathrm{\Gamma}}_{j}=\frac{G{a}^{2}}{\pi}\frac{{K}^{\prime}}{K}ln\frac{{k}^{\prime}}{k}-\frac{{q}_{j}}{2\pi}ln\frac{{r}_{j}K{k}^{\prime}}{a\sqrt{2}}-\underset{l\ne j}{{\sum}^{\prime}}\frac{{q}_{l}^{\ast}}{\pi}ln{\rho}_{jl}+{c}^{\ast}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}z=\pm a,$$(66)$$\frac{2}{\pi {r}_{j}}{\int}_{{\mathrm{\Gamma}}_{j}}c\phantom{\rule{0.16667em}{0ex}}d{\mathrm{\Gamma}}_{j}=\frac{G{a}^{2}}{\pi}\frac{{K}^{\prime}}{K}ln\frac{{k}^{\prime}}{{k}^{2}}-\frac{{q}_{j}}{2\pi}ln\frac{{r}_{j}K{k}^{\prime}}{a\sqrt{2k}}-\underset{l\ne j}{{\sum}^{\prime}}\frac{{q}_{l}^{\ast}}{\pi}ln{\rho}_{jl}+{c}^{\ast}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}z=\pm a+2bi.$$(67)

The relationships (61), (63), (66), and (67) together with the boundary condition (9) and Eq. (21) permit us to express the capillary-tissue fluxes *q*_{1},*q*_{2},…, *q _{M}* and the parameter

It should be noted that in some cases integration of (2) can be performed analytically, in particular when the function *χ*(* _{j}*.) is linear.

Let us consider a capillary distribution in the planes *s*=const with the symmetry shown in Fig. 3(a). The corresponding fundamental rectangle is shown in Fig. 3(b). This configuration can describe different types of structures in the *s*-direction, which will be specified later.

According to the formulas derived in the preceding section, the concentration distribution in the tissue can be written in the form

$$c(x,y)=\frac{G{x}^{2}}{2}+{u}_{1}(\xi ,\eta )-\frac{{q}_{1}}{8\pi}ln\left[{(\xi +1)}^{2}+{\eta}^{2}\right]-\frac{{q}_{2}}{8\pi}ln\left[{\left(\xi -\frac{1}{k}\right)}^{2}+{\eta}^{2}\right]+{c}^{\ast},$$

(68)

where *u*_{1}(*ξ*,*η*) is given by (43) and *ξ*(*x*,*y*),*η*(*x*,*y*) are determined by (31).

The conservation of mass of the substance (21) yields

$${q}_{1}+{q}_{2}=16\mathit{abG}.$$

(69)

Using the formulas (66), (67) and the boundary conditions (9) we obtain

$$\nu {\overline{c}}_{1}=\frac{G\sigma}{2\pi}ln\frac{{k}^{\prime}}{k}-\frac{{q}_{1}}{2\pi}ln\frac{{r}_{1}K{k}^{\prime}}{a\sqrt{2}}-\frac{{q}_{2}}{4\pi}ln\frac{1+k}{k}+{c}^{\ast},$$

(70)

$$\nu {\overline{c}}_{2}=\frac{G\sigma}{2\pi}ln\frac{{k}^{\prime}}{{k}^{2}}-\frac{{q}_{1}}{4\pi}ln\frac{1+k}{k}-\frac{{q}_{2}}{2\pi}ln\frac{{r}_{2}K{k}^{\prime}}{a\sqrt{2k}}+{c}^{\ast},$$

(71)

where the area of the rectangle *ω* is denoted by

$$\sigma =4ab.$$

(72)

In order to simplify the following calculations we assume that the capillary radii are equal, i.e., *r*_{1}
*= r*_{2}
*= r.*

Resolving (70) and (71) with respect to *q*_{1} and *q*_{2}, and using (69) we find

$${q}_{1}=2G\sigma -\frac{\pi \nu}{ln\mu}({\overline{c}}_{1}-{\overline{c}}_{2}),$$

(73)

$${q}_{2}=2G\sigma +\frac{\pi \nu}{ln\mu}({\overline{c}}_{1}-{\overline{c}}_{2}),$$

(74)

$${c}^{\ast}=\frac{G\sigma}{\pi}ln\frac{rK\sqrt{{k}^{\prime}(1+k)}}{a\sqrt{2}}+\frac{\nu}{2ln\mu}\left({\overline{c}}_{1}ln\frac{rK{k}^{\prime}}{a\sqrt{2(1+k)}}+{\overline{c}}_{2}ln\frac{rK{k}^{\prime}\sqrt{k}}{a\sqrt{2(1+k)}}\right).$$

(75)

where

$$\mu =\frac{rK{k}^{\prime}{k}^{1/4}}{a\sqrt{2(1+k)}}.$$

(76)

The following assumptions led Krogh to the solution of the problem [3]:

- Flow in the capillaries is concurrent.
- Axial diffusion in the blood can be neglected.
- The blood is well mixed in the capillary cross-section.
- Axial diffusion in the tissue can be neglected.
- Metabolic consumption (or production) rate is constant.

In addition, Krogh used the boundary condition

$$\frac{\partial c}{\partial \rho}=0\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\rho ={r}_{t},$$

(77)

where *ρ* is the distance from the capillary axis and *r _{t}* is the radius of the tissue cylinder.

Assuming that the ratio *r/r _{t}* is small, we can write Krogh’s result [3] in the form

$$c=\nu \overline{c}-\frac{G{r}_{t}^{2}}{2}\left(ln\frac{{r}_{t}}{r}-0.5\right)+o\phantom{\rule{0.16667em}{0ex}}\left(\frac{r}{{r}_{t}}\right)\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\text{at}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\rho ={r}_{t},$$

(78)

where is the capillary concentration,

Apparently Krogh’s assumptions (2)–(5) are consistent with those we made in Sec. 2.

The case of concurrent flow in capillaries located in the corners of squares as shown in Fig. 4 is described by the relationships (68), (73)–(75) with *q*_{2} = 0 and *a* = *b*. The expression (74) with *q*_{2} = 0 yields the concentration at the point *A* in Fig. 4:

$${c}_{A}=\nu \overline{c}+\frac{2G\sigma}{\pi}ln\frac{rK{k}^{\prime}{k}^{1/4}}{a\sqrt{2(1+k)}},$$

(79)

where the subscript 1 at the capillary concentration is omitted. The relationship (75) yields, for *c**,

$${c}^{\ast}=\nu \overline{c}+\frac{2G\sigma}{\pi}ln\frac{rK{{k}^{\prime}}^{3/4}{k}^{1/4}}{a\sqrt{2}}.$$

(80)

The solution (68) together with (80) gives the concentration at the point *B* in Fig. 4:

$${c}_{B}=\nu \overline{c}+\frac{2G\sigma}{\pi}ln\frac{rK{k}^{\prime}}{2a}.$$

(81)

We introduce the Krogh cylinder radius *r _{t}* so that the cross-sectional area
$\pi {r}_{t}^{2}$ equals the area supplied by each capillary in Fig. 4; this definition gives

$${r}_{t}=\frac{4a}{\sqrt{\pi}}.$$

(82)

Using tables of the complete elliptic integrals [20] to determine the values of *k* and *K,* one can write (79) and (81) in the form

$${c}_{A}=\nu \overline{c}-\frac{2G\sigma}{\pi}\left(ln\frac{2a\sqrt{2}}{r}-0.61\right),$$

(83)

$${c}_{B}=\nu \overline{c}-\frac{2G\sigma}{\pi}\left(ln\frac{2a}{r}-0.44\right),$$

(84)

whereas (78) to the same order can be rewritten as

$$c({r}_{t})=\nu \overline{c}-\frac{2G\sigma}{\pi}\left(ln\frac{4a}{r\sqrt{\pi}}-0.5\right).$$

(85)

It can be shown that for the case *g* > 0 (the substance is consumed in the tissue),

$${c}_{A}<c({r}_{t})<{c}_{B},$$

(86)

which could be expected from the geometrical picture shown in Fig. 4.

If the Krogh cylinder radius is chosen as ${r}_{t}^{\ast}=2a\sqrt{2}$, then $c({r}_{t}^{\ast})<{c}_{A}$, where $c({r}_{t}^{\ast})$ is calculated using (78); if ${r}_{t}^{\ast \ast}=2a$, then $c({r}_{t}^{\ast})>{c}_{B}$. Therefore we have the useful estimates

$$c({r}_{t}^{\ast})<{c}_{A}<c({r}_{t})<{c}_{B}<c({r}_{t}^{\ast \ast}).$$

(87)

For simplicity in the present paper we restrict ourselves to the case of linear functional relationships between *χ* and * _{j}* in (2). The particular case

$$\chi ({\overline{c}}_{j})=\varphi {\overline{c}}_{j},$$

(88)

$${\alpha}_{j}=\frac{D}{{Q}_{j}(1+\varphi )},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}j=1,2,$$

(89)

we rewrite the equations (2) in the form

$$\frac{d{\overline{c}}_{j}}{ds}={\alpha}_{j}{q}_{j},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}j=1,2.$$

(90)

Substitution of (73) and (74) into (90) yields a system of two linear ordinary differential equations with respect to _{1} and _{2}.

If *α*_{1} + *α*_{2}≠0, the general solution of (90) can be expressed in the form

$${\overline{c}}_{1}=A{e}^{\lambda s}-4G\sigma \frac{{\alpha}_{1}{\alpha}_{2}}{{\alpha}_{1}+{\alpha}_{2}}s+\frac{2G\sigma}{\pi \nu}\frac{{\alpha}_{1}({\alpha}_{1}-{\alpha}_{2})}{{({\alpha}_{1}+{\alpha}_{2})}^{2}}ln\mu +B,$$

(91)

$${\overline{c}}_{2}=-A\frac{{\alpha}_{2}}{{\alpha}_{1}}{e}^{\lambda s}-4G\sigma \frac{{\alpha}_{1}{\alpha}_{2}}{{\alpha}_{1}+{\alpha}_{2}}s-\frac{2G\sigma}{\pi \nu}\frac{{\alpha}_{2}({\alpha}_{1}-{\alpha}_{2})}{{({\alpha}_{1}+{\alpha}_{2})}^{2}}ln\mu +B,$$

(92)

where

$$\lambda =\frac{\pi \nu ({\alpha}_{1}+{\alpha}_{2})}{ln\mu},$$

(93)

and *A, B* are arbitrary constants. In the case *α*_{1} + *α*_{2}=0, which according to (89) describes countercurrent flow in the capillaries with equal absolute values of the blood flow, the general solution of (90) is

$${\overline{c}}_{1}=-\frac{2G\sigma \pi \nu {\alpha}^{2}}{ln\mu}{s}^{2}-2\alpha \left(G\sigma -\frac{\pi \nu A}{ln\mu}\right)s+B+A,$$

(94)

$${\overline{c}}_{2}=-\frac{2G\sigma \pi \nu {\alpha}^{2}}{ln\mu}{s}^{2}+2\alpha \left(G\sigma +\frac{\pi \nu A}{ln\mu}\right)s+B-A,$$

(95)

where *α* = *α*_{1} = −*α*_{2}, and *A, B* are arbitrary constants.

We consider a microcirculatory unit with concurrent capillary blood flow shown in Fig. 5(a) (the three-dimensional capillary structures similar to those shown in Fig. 5 were considered earlier by Grunewald; cf. [14).

Rectangular cylinders Ω representing three types of the microcirculatory units with concurrent and countercurrent flow. The horizontal arrows denote arterial inflow *a* and venous outflow *υ*.

If the arterial ends of the capillaries are located in the plane *s* = 0 and the concentration of the substance at the capillary inlets is *c*_{0},, the solution (91) (92) assumes the form

$${\overline{c}}_{1}={c}_{0}-4G\sigma \frac{{\alpha}_{1}{\alpha}_{2}}{{\alpha}_{1}+{\alpha}_{2}}s+\frac{2G\sigma}{\pi \nu}\frac{{\alpha}_{1}({\alpha}_{1}-{\alpha}_{2})}{{({\alpha}_{1}+{\alpha}_{2})}^{2}}ln\mu (1-{e}^{\lambda s}),$$

(96)

$${\overline{c}}_{2}={c}_{0}-4G\sigma \frac{{\alpha}_{1}{\alpha}_{2}}{{\alpha}_{1}+{\alpha}_{2}}s-\frac{2G\sigma}{\pi \nu}\frac{{\alpha}_{2}({\alpha}_{1}-{\alpha}_{2})}{{({\alpha}_{1}+{\alpha}_{2})}^{2}}ln\mu (1-{e}^{\lambda s}),$$

(97)

for 0 ≤ *s* ≤ *L.* It is easy to show that in this case the concentration varies monotonically along the capillaries.

If the blood flow rates in both capillaries are identical, so that *α*_{1}=*α*_{2}= *α*, then (96), (97) reduce to

$${\overline{c}}_{1}={\overline{c}}_{2}={c}_{0}-2G\sigma \alpha s$$

(98)

in accordance with the Krogh-cylinder-model solution [1].

The concentration at the venous end of the capillary equals

$${\overline{c}}_{1}(L)={\overline{c}}_{2}(L)={c}_{0}-2G\sigma \alpha L.$$

(99)

The relationship (98) can easily be obtained by consideration of the mass balance between a capillary and the corresponding volume of tissue supplied by this capillary. In contrast, the solution (96) (97) reflects more complicated phenomena, namely, the interaction between the adjacent capillaries. The formulas (96) and (97) conclude the solution of the three-dimensional mass-transport problem for the geometry specified in Fig. 5(a); substitution of (96) and (97) into (68) together with (73)–(75) makes it possible to calculate the spatial distribution *c*(*x*,*y*,*s*).

The derivation of the numerical values of the concentration at an arbitrary point in the fundamental rectangle involves computations of the Jacobian elliptic functions and complete elliptic integral. Rapidly converging series expansions can be used for this purpose [20]; detailed tables of these functions are also available [21].

The countercurrent flow conditions are realized in the microcirculatory unit shown in Fig. 5(b). We assume that the blood flow rates in both capillaries are identical; *α*_{1} + *α*_{2}=0. Placing the origin of the coordinate: system at the midpoint between the arterial and venous ends of the capillaries, applying the boundary conditions

$${\overline{c}}_{1}\left(-\frac{L}{2}\right)={c}_{0},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\overline{c}}_{2}\left(\frac{L}{2}\right)={c}_{0},$$

(100)

$${\overline{c}}_{1}={c}_{0}-2G\sigma \alpha \left(s+\frac{L}{2}\right)-\frac{2G\sigma \pi \nu {\alpha}^{2}}{ln\mu}\left({s}^{2}-\frac{{L}^{2}}{4}\right),$$

(101)

$${\overline{c}}_{2}={c}_{0}+2G\sigma \alpha \left(s-\frac{L}{2}\right)-\frac{2G\sigma \pi \nu {\alpha}^{2}}{ln\mu}\left({s}^{2}-\frac{{L}^{2}}{4}\right).$$

(102)

The concentration at the venous ends of the capillaries is

$${\overline{c}}_{1}\left(\frac{L}{2}\right)={\overline{c}}_{2}\left(-\frac{L}{2}\right)={c}_{0}-2G\sigma \alpha L$$

(103)

[cf. (98)].

In contrast to the case of concurrent flow (96), (97), the concentrations (101), (102) may vary non-monotonically along the capillaries. If the extremum exists, then it is located at the points

$${s}_{1m}=-\frac{ln\mu}{2\pi \alpha},\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{s}_{2m}=\frac{ln\mu}{2\pi \alpha}$$

(104)

for the first and second capillary, respectively.

Since (76) yields *μ* < 1, it follows from (104) that

$${s}_{1m}>0,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{s}_{2m}<0.$$

(105)

Therefore in the case of fully countercurrent flow the capillary concentration can achieve an extremal value at the point between the midpoint of the capillary and the venous end, but not between the arterial end and the midpoint.

In the microcirculatory unit shown in Fig. 5(c), the flow in the capillaries is concurrent or countercurrent depending on the spatial region.

The solutions (91), (92) and (94), (95) enable us to describe a situation with different capillary-inlet concentrations and capillary blood flow rates. For brevity, we consider only the case of equal capillary-inlet concentrations and equal blood flow rates. We write out the values of the capillary concentration at certain points of the unit, separated by the distance *L*/2 in the *s*-direction:

$$\begin{array}{lll}1\text{st}\phantom{\rule{0.16667em}{0ex}}\text{capillary}\hfill & 2\text{nd}\phantom{\rule{0.16667em}{0ex}}\text{capillary}\hfill & \phantom{\rule{0.16667em}{0ex}}\hfill \\ \overline{c}(0+)={c}_{0}-\frac{G\sigma \alpha L}{2}(5-{e}^{\lambda L/2}),\hfill & {\overline{c}}_{2}(0)={c}_{0}-\frac{G\sigma \alpha L}{2}(3-{e}^{\lambda L/2})\hfill & (106)\hfill \\ {\overline{c}}_{1}(\frac{L}{2})={c}_{0}-\frac{G\sigma \alpha L}{2}(3-{e}^{\lambda L/2}),\hfill & {\overline{c}}_{2}(\frac{L}{2}-)={c}_{0}-\frac{G\sigma \alpha L}{2}(5-{e}^{\lambda L/2}),\hfill & (107)\hfill \\ \phantom{\rule{0.16667em}{0ex}}\hfill & {\overline{c}}_{2}(\frac{L}{2}+)={c}_{0}-\frac{G\sigma \alpha L}{2}(3+{e}^{\lambda L/2}),\hfill & (108)\hfill \\ {\overline{c}}_{1}(L)={c}_{0},\hfill & {\overline{c}}_{2}(L)={c}_{0}-G\sigma \alpha L,\hfill & (109)\hfill \\ {\overline{c}}_{1}({\scriptstyle \frac{3}{2}}L)={c}_{0}-G\sigma \alpha L,\hfill & {\overline{c}}_{2}({\scriptstyle \frac{3}{2}}L)={c}_{0},\hfill & (110)\hfill \\ {\overline{c}}_{1}(2L-)={c}_{0}-\frac{G\sigma \alpha L}{2}(3+{e}^{\lambda L/2}),\hfill & {\overline{c}}_{2}(2L)={c}_{0}-\frac{G\sigma \alpha L}{2}(3-{e}^{\lambda L/2}).\hfill & (111)\hfill \end{array}$$

Here the notation (*a* +) or (*a* −) means

$$\overline{c}(a+)=\underset{\begin{array}{l}s\to a\\ s>a\end{array}}{lim}\overline{c}(s),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\overline{c}(a-)=\underset{\begin{array}{l}s\to a\\ s<a\end{array}}{lim}\overline{c}(s).$$

Since the capillary structure is periodic in the *s*-direction with period 2*L*, the functions _{1}(*s*) and _{2}(*s*) for 0 ≤ *s* ≤ 2*L* completely solve the problem.

The expressions (106.1), (107.2), (108), and (111.1) show that the blood going into venules is a mixture of blood from two kinds of capillaries with the concentrations

$${c}_{0}-\frac{G\sigma \alpha L}{2}(3+{e}^{\lambda L/2}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{c}_{0}-\frac{G\sigma \alpha L}{2}(5-{e}^{\lambda L/2}).$$

(112)

Since the blood flow rate in all capillaries is the same, the average concentration in the venular blood can be found as the arithmetic average of (112):

$${c}_{\upsilon}={c}_{0}-2G\sigma \alpha L$$

(113)

[cf. (98) and (103)], which is in accordance with the Fick law.

A model of steady-state mass transport in tissue has been developed, which can be considered as a generalization of the Krogh cylinder model.

The solution has been obtained with the assumption that the capillary blood vessels are parallel to each other and that their distribution in a plane *s* = const is double periodic, so that the diffusion problem can be stated in a rectangle instead of the infinite plane. It was also assumed that the capillary radii are small in comparison with the intercapillary distances, and the latter are in turn small in comparison with the capillary lengths; the latter assumption permits us to neglect the axial diffusion in the tissue.

Practical computations of the concentration distribution can be done according to the following format:

- Knowing the ratio of the sides of the fundamental rectangle makes it possible to find the parameter
*k*from Eq. (26). - Equations (61) (63), (66), (67), together with (21) and the boundary condition (9), yield
*M*linear relationships between the capillary concentrations_{1},_{2}, …,, capillary-tissue diffusion fluxes_{M}*q*_{1},*q*_{2}, …,*q*, and the parameter_{M}*c**, where*M*is the number of capillaries with different values of concentration in the fundamental rectangle. These linear relationships and the additional relationship (21) can be used to express*q*_{1},*q*_{2}, …,*q*and_{M}*c** as linear functions of_{1},_{2}, …,._{M} - Substitution of these expressions for the mass fluxes
*q*=_{j}*J*/_{j}*D*into the convective transport equations (2) leads to a linear [if χ() is a linear function] or non-linear set of ordinary differential equations which can be solved with the boundary conditions of the type (10) or (11). The solution_{j}_{1}(*s*),_{2}(*s*), …,(_{M}*s*) also determines the functions*q*_{1}(*s*),*q*_{2}(*s*), …,*q*(_{M}*s*) and*c**(*s*).

The solution presented is only valid if *c* > 0, since this condition enables us to assume a constant consumption/production rate *g*. If *c* = 0 in certain domains, the problem becomes non-linear, and the derived solution is no longer valid.

The present model makes it possible to study the mass transport of various substances (glucose, O_{2}, CO_{2}, etc.) in the tissue. Further, the model takes into account such important physiological effects as interaction between capillaries with concurrent and countercurrent flows, interaction between “fast” and “slow” capillaries, and mass transport between tissue regions supplied by different arteriolar vessels.

The author is indebted to Dr. J, F. Gross, Dr. H. D. Papenfuss, and Dr. H. S. Chen for helpful discussions and critical comments, and to Mrs. Rhoda G. Miller for typing the manuscript.

This study was supported by NIH Grants HL 17421 and NO I-CB-63981.

- 2
*a,*2*b* - dimensions of the fundamental rectangle
*A, B*- constants in (91), (92) and (94), (95)
*c*- concentration in the tissue
*c*_{0},*c*_{j}_{0}- inlet concentration
*c*,_{A}*c*_{B}- defined by (79) and (81)
*c*_{υ}- outlet concentration
_{j}- intracapillary concentration
*c**- parameter, introduced in Eq. (20)
*D*- diffusion coefficient
*D*_{s}- diffusion coefficient in the
*s*-direction *F*- function in Eq. (20)
*F**- defined in Eq. (38)
*g*- consumption/production rate
*G*- defined in Eq. (13)
*J*_{j}- capillary-tissue diffusion flux
*k*- argument of the complete elliptic integral
*K*- complete elliptic integral
*K′, K′*- defined by Eq. (27)
*L, L*_{j}- capillary length
*L**- characteristic capillary length
*M*- number of capillaries in the fundamental rectangle
*N*- Neumann function, (Eq. 20)
*N**- defined in Eq. (39)
- o,
*O* - orders of magnitude
*q*_{j}- defined in Eq. (13)
- ${q}_{j}^{\ast}$
- defined in Eq. (22)
*Q*_{j}- capillary volumetric blood flow rate
*r, r*_{j}- radius of the capillary
*r**- characteristic intercapillary distance
*r*${r}_{t}^{\ast},{r}_{t}^{\ast \ast}$_{t},- Krogh cylinder radius
*R*_{j}- defined in Eqs. (60), (64), and (65)
- sn, cn, dn
- Jacobian elliptic functions
*S, S*_{b}- solubility coefficients
*u*- defined in Eq. (16)
*u*_{1},*u*_{2}- defined in Eq. (42)
*υ*- defined in Eq. (49)
*x, y, s*- rectangular coordinates
*X, Y, Z*- defined in Eq. (30)
*z*- complex variable

*α*_{j}- defined by (89)
- Γ
- boundary of the rectangle
*ω* - Γ
_{j} - capillary-tissue interface
- Γ*
- image of Γ
*δ*- Dirac delta function
*ζ*- complex variable
- complex conjugate of
*ζ* *λ*- defined by (93)
*μ*- defined by (76)
*ν*- ratio of the solubility coefficients
*ξ, η*- defined by (23)
*ρ*- distance from the capillary axis
*ρ*_{jl},_{jl}- defined in Eqs. (62)
*σ*- area of
*ω* *ϕ*- defined in Eq. (88)
*χ*- concentration of the substance bound to hemoglobin
*ω*- fundamental rectangle
*ω**- image of
*ω* *Ω*- rectangular cylinder

1. Middleman S. Transport Phenomena in the Cardiovascular System. Chapter 3 Wiley-Interscience; New York: 1972.

2. Leonard EF, Jorgensen SB. The analysis of convection and diffusion in capillary beds. Annual Rev Biophys Bioeng. 1974;3:293–339. [PubMed]

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

4. Thews G. Über die mathematische Behandlung physiologischer Diffusionsprozess in Zylinderförmigen Objekten. Acta Biotheor (Leiden) 1953;10:105–136.

5. Reneau DD, Jr, Bruley DF, Knisely MH. A mathematical simulation of oxygen release, diffusion, and consumption in the capillaries and tissue of the human brain. In: Hershey D, editor. Chemical Engineering in Medicine and Biology. Plenum; New York: 1967. pp. 135–241.

6. Hyman WA. A simplified model of the oxygen supply function of capillary blood flow. Advances Exp Med Biol. 1973;37b:835–841. [PubMed]

7. Stewart RR, Morrazzi CA. Oxygen transport in the human brain—analytical solution. Advances Exp Med Biol. 1973;37b:843–848. [PubMed]

8. Wilson TA. Capillary spacing for minimum entropy production. J Theoret Biol. 1966;11:436–445. [PubMed]

9. Diemer K. Eine verbesserte Modellvorstellung zur O_{2}-Versorgung des Gehirns. Naturwissenschaften. 1963;50:617–618.

10. Bailey HR. Oxygen exchange between capillary and tissue: some equations describing countercurrent and nonlinear transport. In: Reeve EB, Guyton AC, editors. Physical bases of circulatory transport. Regulation and Exchange. Saunders; Philadelphia: 1967. pp. 353–366.

11. Reneau DD, Knisely MH. A mathematical simulation of oxygen transport in the human brain under conditions of countercurrent capillary blood flow. Chem Engineering Progr. 1971;67:18–27.

12. Metzger H. Distribution of oxygen partial pressure in a two-dimensional tissue supplied by capillary meshes and concurrent and countercurrent systems. Math Biosci. 1969;5:143–154.

13. Metzger H. Geometric considerations in modelling oxygen transport processes in tissue. Advances Exp Med Biol. 1973;37b:761–772. [PubMed]

14. Grunewald WA. Computer calculation for tissue oxygenation and the meaningful presentation of the results. Advances Exp Med Biol. 1973;37b:783–792. [PubMed]

15. Grunewald WA, Sowa W. Capillary structures and O_{2} supply to tissue. An analysis with a digital diffusion model as applied to the skeletal muscle. Rev Physiol Biohem Pharmacol. 1977;77:149–209. [PubMed]

16. Aroesty J, Gross JF. Convection and diffusion in the microcirculation. Microvasc Res. 1970;2:247–267. [PubMed]

17. Morse PM, Feshbach H. Methods of Theoretical Physics. McGraw-Hill; New York: 1953. p. 2.

18. Bateman H, Erdélyi A. Higher Transcendental Functions. Vol. 3. McGraw-Hill; New York: 1955.

19. Bierens de Haan D. Nouvelles Tables D’Intégrales Définies. G. E. Stechert; New York: 1939.

20. Abramowitz M, Stegun IA. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Chapters 16, 17 Dover; New York: 1972.

21. Milne-Thomson LM. Jacobian Elliptic Function Tables. Dover; New York: 1956.

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. |