PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Math Biosci. Author manuscript; available in PMC 2017 September 22.
Published in final edited form as:
PMCID: PMC5609729
NIHMSID: NIHMS857492

Analysis of Capillary-Tissue Diffusion in Multicapillary Systems

Abstract

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.

1. INTRODUCTION

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.

2. MATHEMATICAL STATEMENT OF THE PROBLEM

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

Fig. 1
Three-dimensional capillary network. Arrows indicate blood flow directions.

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

D(2cx2+2cy2)+Ds2cs2-g=0,
(1)

where c is the local concentration of the substance, D is the diffusion coefficient in a plane perpendicular to the capillaries, Ds is the diffusion coefficient in the direction along the capillaries, and g is the local consumption rate (g > 0) or production rate (g < 0) of the substance.

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

Qjdds[c¯j+χ(c¯j)]=-Jj,
(2)

where the index j identifies the capillary, Qj = πrjūj is the volumetric blood flow rate in the capillary j, rj is the capillary radius, ūj is the average velocity of blood in the capillary, cj is the intracapillary concentration of the substance dissolved in the plasma, χ(cj) is the concentration of the substance bound to hemoglobin, and Jj is the volume flux of the substance through the capillary wall per unit length of the capillary.

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

  1. Diffusion in the s-direction in the tissue can be neglected. This can be justified if either Ds = 0 or r* [double less-than sign] L*, where r* is the characteristic intercapillary distance and L* is the characteristic capillary length (cf. [15]).
  2. The consumption/production rate is constant, g = const.
  3. The capillary radii are small in comparison with the intercapillary distances, rj [double less-than sign] r*.
  4. 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

|(x-xj)2+(y-yj)2|rj,sjssj+Lj,
(3)

where xj,yj are the coordinates of the center of the capillary, sj is the coordinate of the “lower” end of the capillary, and Lj is the capillary length. If Qj is the capillary blood flow (Qj > 0 if the flow is in the positive s-direction, Qj <0 if it is in the negative s-direction), and cj0 is the concentration of the substance at the capillary inlet (i.e., at the arterial end of the capillary), the set of geometrical and physical parameters characterizing the capillary will be written as

(xj,yj,sj,Lj,Qj,cj0).
(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 ≤ 2a, 0 ≤ y ≤ 2b, −∞ < s < ∞) described by the set (4) there is an infinite number of capillaries outside of Ω, characterized by the sets

(±xj+4am,±yj+4bn,sj,Lj,Qj,cj0),
(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,

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).
(6)

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

cn=0.
(7)

On the capillary-tissue interface Γj,

Jj=-2πrjDΓjcndΓj,
(8)
c=νc¯j.
(9)

Here n is the unit vector normal to Γj pointing in the external direction, and v = S/Sb where S and Sb are the solubility coefficients of the substance in the tissue and blood, respectively.

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

c¯j(sj)=cj0ifQj>0,
(10)

c¯j(sj+Lj)=cj0ifQj<0,
(11)

where cj0 is the concentration of the substance at the capillary inlet. In the following sections we will seek an approximate solution of Eqs. (1) and (2) with the boundary conditions (7)–(9) and (10) or (11).

3. SOLUTION OF THE DIFFUSION EQUATION

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 Ji per unit length. Therefore, Eq. (1) can be rewritten in the form

2cx2+2cy2-G+jqjδ(x-xj,y-yj)=0,
(12)

where

G=gD,qj=JjD,
(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

νc¯j=1ΓjΓjcdΓj,
(14)

where Γ̃j denotes the portion of the capillary-tissue interface confined within ω, and |Γ̃j| denotes the length of the corresponding arch.

A. SOLUTION OF THE BOUNDARY-VALUE PROBLEM

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

xa,0y2b.
(15)

Substitution of

c=Gx22+u(x,y)
(16)

into (12) leads to the Poisson equation

2ux2+2uy2=-jqjδ(x-xj,y-yj)
(17)

with the Neumann-type boundary conditions

un=0aty=0,y=2b,
(18)

un=Gaatx=±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)=-Γu(z)nN(z,z)dΓ+ωF(z)N(z,z)dxdy+c,
(20)

where F(z) = Σqjδ(zzj), c* is an arbitrary constant, and the second integral in (20) is taken over the area of the rectangle ω.

The necessary condition for the resolution of the Neumann problem

jqj=4abG,
(21)

where

qj={qjifzjω,zjΓ,qj2ifzjΓ,zj±a,±a+2bi,qj4ifzj=±a,±a+2bi.
(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

ζ=ξ+iη.
(23)

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

ζ=sn(zKa,k),
(24)

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

K(k)=01dt(1-t2)(1-k2t2),
(25)

and the parameter k is a root of the transcendental equation

K(1-k2)K(k)=2ba,0<k<1.
(26)

In the following equations the standard notation used will be

k=1-k2,K(k)=K(k).
(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=aK0ζdt(1-t2)(1-k2t2).
(28)

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

z=0ζ=0,z=±aζ=±1,z=±a+2biζ=±1k,z=2biζ=
(29)

correspond to each other.

Introducing the notation

Z=X+iY=zKa,
(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

ξ=sn(X,k)dn(Y,k)cn2(Y,k)+k2sn2(X,k)sn2(Y,k),η=cn(X,k)dn(X,k)sn(Y,k)cn(Y,k)cn2(Y,k)+k2sn2(X,k)sn2(Y,k).
(31)

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

2uξ2+2uη2=-jqjδ(ζ-ζj),
(32)

where ζj = sn(zjK/a,k). The property of the delta function

δ(ζ-ζj)=|dzdζ|2δ(z-zj)
(33)

has been used.

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

(un)Γ=(un)ΓdzdζΓ,
(34)

dΓ=dzdζΓdΓ,
(35)

dxdy=dzdζ2dξdη,
(36)

we can rewrite the solution (20) in the form

u(ζ)=-Γ(un)ΓN(ζ,ζ)dΓ+ωF(ζ)N(ζ,ζ)dξdη,
(37)

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

F(ξ)=-jqjδ(ζ-ζj),
(38)
N(z,z)=N(snzKa,snzKa).
(39)

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

N(ζ,ζ)=-12πlnζ-ζζ-ζ¯,
(40)

where [zeta] is the complex conjugate of ζ.

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

|dzdζ|Γ={aK1(ξ2-1)(1-k2ξ2)as1<ξ<1k,aK1(1-ξ2)(1-k2ξ2)asξ<1orξ>1k.
(41)

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

u(ζ)=u1(ζ)+u2(ζ)+c,
(42)

where

u1=Ga22πK01ln{[(1-k2t2k-ξ)2+η2][(1-k2t2k+ξ)2+η2]}(1-t2)(1-k2t2)dt,
(43)
u2=-jqj2πln[(ξ-ξj)2+(η-ηj)2(ξ-ξj)2+(η+ηj)2].
(44)

Therefore the concentration c is given by

c=Gx22+u1(ξ,η)+u2(ξ,η)+c,
(45)

where u1 and u2 are expressed by (43) and (44), and ξ(x,y) and η(x,y) are determined by (31). The parameters qj and c* are functions of the axial coordinate s. They will be determined later.

It can be shown that

c=limcasz2bi.
(46)

In order to verify (46) we notice that as z→2bi, ζ → ∞ and the functions u1 and u2 can be expanded as

u1=Ga2πKKln(ξ2+η2)+O(1ξ2+η2),
(47)
u2=jqj2πln(ξ2+η2)+O(1ξ2+η2).
(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→2bi, 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

Fig. 2
Points in the fundamental rectangle ω for which integral (43) is calculated are marked by x.

Ga2πKKυ=Gx22+u1,
(49)

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

z=0,ζ=0υ=ln1k,
(50)

z=±a,ζ=±1:υ=lnkk,
(51)

z=±a+2bi,ζ=±1k:υ=lnkk2,
(52)

z=bi,ζ=ik:υ=12ln2(1+k)k2k,
(53)

z=±a+bi,ζ=±1k:υ=12ln2(1-k)k2k,
(54)

z=±a2,ζ=±11+k:υ=12ln2(1-k)kkk4,
(55)

z=±a2+bi,ζ=±1+k+i1-k2k:υ=12ln2(1-k)k2k,
(56)

z=±a2+2bi,ζ=±11-k:υ=12ln2(1+k)kkk4.
(57)

B. CALCULATION OF THE LINE INTEGRALS (14)

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 = zj + rje in the z-plane is transformed by the conformal mapping (24) into a curve

ζ=sn(ζj+rjKaeiθ,k),
(58)

we can use the condition rjK/a[double less-than sign]1, and expand the function (58) in a Taylor series as

ζ=ζj+rjKa|cnzjKadnzjKa|ei(θ-θj)+12(rjKa)2|snzjKa(dn2zjKa+k2cn2zjKa)|e2i(θ-θj)+,
(59)

where θj and θj are constants.

We consider three different cases:

  1. zj [set membership]ω, zj [negated set membership] Γ. The relationship (59) shows that the image of the circle Γj is to the order o(rj/a) a circle Γj:
    ζ-ζj=Rj,Rj=rjKa|cnzjKadnzjKa|.
    (60)
    Substitution of (60) into (45) gives to the lowest order with respect to the small parameter rj/a
    12πrjΓjcdΓj=Gxj22+u1(ζj)-qj2πln(2Rjηj)-ljql2πlnρjlρ¯jl+c,
    (61)
    where
    ρjl=(ξj-ξl)2+(ηj-ηl)2,ρ¯jl=(ξj-ξl)2+(ηj+ηl)2.
    (62)
  2. zj [set membership] Γ, zj ≠ ±a, ±a + 2bi. In this case, the semi-circle Γj is mapped approximately into a semi-circle Γj in the ζ-plane with the radius Rj determined by (60). The integral in the boundary condition (14) is calculated as
    1πrjΓjcdΓj=Gxj22+u1(ζj)-qj2πlnRj-ljqlπlnρjl+c.
    (63)
  3. zj = ±a, ±a + 2bi. In this case, the quarter of the circle Γj is mapped approximately into semi-circle Γj in the ζ-plane:
    z±a=rjζ±1=Rj,Rj=(rjKka2)2,
    (64)
    z±a-2bi=rj|ζ±1k|=Rj,Rj=(rjKka2k)2.
    (65)
    The average concentration can be calculated using the relationships
    2πrjΓjcdΓj=Ga2πKKlnkk-qj2πlnrjKka2-ljqlπlnρjl+catz=±a,
    (66)
    2πrjΓjcdΓj=Ga2πKKlnkk2-qj2πlnrjKka2k-ljqlπlnρjl+catz=±a+2bi.
    (67)

C. SUMMARY

The relationships (61), (63), (66), and (67) together with the boundary condition (9) and Eq. (21) permit us to express the capillary-tissue fluxes q1,q2,…, qM and the parameter c* in terms of the capillary concentrations c1,c2,…, cM, where M is the number of capillaries in the domain ω. Substitution of these expressions into the set of equations (2) yields a set of M ordinary differential equations with respect to c1,c2,…, cM. The boundary conditions for this system are given by (10) and (11).

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

4. EXAMPLES

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.

Fig. 3
Double periodic capillary distribution in the xy plane and the fundamental rectangle ω.

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

c(x,y)=Gx22+u1(ξ,η)-q18πln[(ξ+1)2+η2]-q28πln[(ξ-1k)2+η2]+c,
(68)

where u1(ξ,η) is given by (43) and ξ(x,y),η(x,y) are determined by (31).

The conservation of mass of the substance (21) yields

q1+q2=16abG.
(69)

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

νc¯1=Gσ2πlnkk-q12πlnr1Kka2-q24πln1+kk+c,
(70)

νc¯2=Gσ2πlnkk2-q14πln1+kk-q22πlnr2Kka2k+c,
(71)

where the area of the rectangle ω is denoted by

σ=4ab.
(72)

In order to simplify the following calculations we assume that the capillary radii are equal, i.e., r1 = r2 = r.

Resolving (70) and (71) with respect to q1 and q2, and using (69) we find

q1=2Gσ-πνlnμ(c¯1-c¯2),
(73)

q2=2Gσ+πνlnμ(c¯1-c¯2),
(74)

c=GσπlnrKk(1+k)a2+ν2lnμ(c¯1lnrKka2(1+k)+c¯2lnrKkka2(1+k)).
(75)

where

μ=rKkk1/4a2(1+k).
(76)

A. COMPARISON WITH THE KROGH CYLINDER MODEL

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

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

In addition, Krogh used the boundary condition

cρ=0atρ=rt,
(77)

where ρ is the distance from the capillary axis and rt is the radius of the tissue cylinder.

Assuming that the ratio r/rt is small, we can write Krogh’s result [3] in the form

c=νc¯-Grt22(lnrtr-0.5)+o(rrt)atρ=rt,
(78)

where c 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 q2 = 0 and a = b. The expression (74) with q2 = 0 yields the concentration at the point A in Fig. 4:

cA=νc¯+2GσπlnrKkk1/4a2(1+k),
(79)

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

Fig. 4
Krogh tissue cylinder and the square capillary array.
c=νc¯+2GσπlnrKk3/4k1/4a2.
(80)

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

cB=νc¯+2GσπlnrKk2a.
(81)

We introduce the Krogh cylinder radius rt so that the cross-sectional area πrt2 equals the area supplied by each capillary in Fig. 4; this definition gives

rt=4aπ.
(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

cA=νc¯-2Gσπ(ln2a2r-0.61),
(83)

cB=νc¯-2Gσπ(ln2ar-0.44),
(84)

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

c(rt)=νc¯-2Gσπ(ln4arπ-0.5).
(85)

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

cA<c(rt)<cB,
(86)

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

If the Krogh cylinder radius is chosen as rt=2a2, then c(rt)<cA, where c(rt) is calculated using (78); if rt=2a, then c(rt)>cB. Therefore we have the useful estimates

c(rt)<cA<c(rt)<cB<c(rt).
(87)

B. SOLUTION OF THE CONVECTIVE TRANSPORT EQUATIONS

For simplicity in the present paper we restrict ourselves to the case of linear functional relationships between χ and cj in (2). The particular case χ=0 corresponds to the condition in which the substance is not bound to hemoglobin. Using the notation

χ(c¯j)=ϕc¯j,
(88)

αj=DQj(1+ϕ),j=1,2,
(89)

we rewrite the equations (2) in the form

dc¯jds=αjqj,j=1,2.
(90)

Substitution of (73) and (74) into (90) yields a system of two linear ordinary differential equations with respect to c1 and c2.

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

c¯1=Aeλs-4Gσα1α2α1+α2s+2Gσπνα1(α1-α2)(α1+α2)2lnμ+B,
(91)

c¯2=-Aα2α1eλs-4Gσα1α2α1+α2s-2Gσπνα2(α1-α2)(α1+α2)2lnμ+B,
(92)

where

λ=πν(α1+α2)lnμ,
(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

c¯1=-2Gσπνα2lnμs2-2α(Gσ-πνAlnμ)s+B+A,
(94)

c¯2=-2Gσπνα2lnμs2+2α(Gσ+πνAlnμ)s+B-A,
(95)

where α = α1 = −α2, and A, B are arbitrary constants.

C. CONCURRENT FL0W

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

Fig. 5
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 c0,, the solution (91) (92) assumes the form

c¯1=c0-4Gσα1α2α1+α2s+2Gσπνα1(α1-α2)(α1+α2)2lnμ(1-eλs),
(96)

c¯2=c0-4Gσα1α2α1+α2s-2Gσπνα2(α1-α2)(α1+α2)2lnμ(1-eλs),
(97)

for 0 ≤ sL. 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

c¯1=c¯2=c0-2Gσαs
(98)

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

The concentration at the venous end of the capillary equals

c¯1(L)=c¯2(L)=c0-2Gσα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].

D. FULLY COUNTERCURRENT FLOW

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

c¯1(-L2)=c0,c¯2(L2)=c0,
(100)

and using (94) and (95) we obtain

c¯1=c0-2Gσα(s+L2)-2Gσπνα2lnμ(s2-L24),
(101)
c¯2=c0+2Gσα(s-L2)-2Gσπνα2lnμ(s2-L24).
(102)

The concentration at the venous ends of the capillaries is

c¯1(L2)=c¯2(-L2)=c0-2Gσα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

s1m=-lnμ2πα,s2m=lnμ2πα
(104)

for the first and second capillary, respectively.

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

s1m>0,s2m<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.

E. AN EXAMPLE OF MORE COMPLEX GEOMETRY

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:

1stcapillary2ndcapillaryc¯(0+)=c0-GσαL2(5-eλL/2),c¯2(0)=c0-GσαL2(3-eλL/2)(106)c¯1(L2)=c0-GσαL2(3-eλL/2),c¯2(L2-)=c0-GσαL2(5-eλL/2),(107)c¯2(L2+)=c0-GσαL2(3+eλL/2),(108)c¯1(L)=c0,c¯2(L)=c0-GσαL,(109)c¯1(32L)=c0-GσαL,c¯2(32L)=c0,(110)c¯1(2L-)=c0-GσαL2(3+eλL/2),c¯2(2L)=c0-GσαL2(3-eλL/2).(111)

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

c¯(a+)=limsas>ac¯(s),c¯(a-)=limsas<ac¯(s).

Since the capillary structure is periodic in the s-direction with period 2L, the functions c1(s) and c2(s) for 0 ≤ s ≤ 2L 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

c0-GσαL2(3+eλL/2),c0-GσαL2(5-eλ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υ=c0-2GσαL
(113)

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

5. CONCLUDING REMARKS

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:

  1. Knowing the ratio of the sides of the fundamental rectangle makes it possible to find the parameter k from Eq. (26).
  2. Equations (61) (63), (66), (67), together with (21) and the boundary condition (9), yield M linear relationships between the capillary concentrations c1, c2, …, cM, capillary-tissue diffusion fluxes q1, q2, …, qM, and the parameter 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 q1, q2, …, qM and c* as linear functions of c1, c2, …, cM.
  3. Substitution of these expressions for the mass fluxes qj = Jj/D into the convective transport equations (2) leads to a linear [if χ(cj) 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 c1(s), c2(s), …, cM(s) also determines the functions q1(s), q2(s), …, qM(s) and c*(s).
  4. The general solution (45) gives the three-dimensional distribution c(x, y, s). It should be noted that the integral (43) is given by one of the simple expressions (50)–(57) at 15 points of the fundamental rectangle marked in Fig. 2.

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, O2, CO2, 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.

Acknowledgments

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.

6. NOMENCLATURE

2a, 2b
dimensions of the fundamental rectangle
A, B
constants in (91), (92) and (94), (95)
c
concentration in the tissue
c0, cj0
inlet concentration
cA, cB
defined by (79) and (81)
cυ
outlet concentration
cj
intracapillary concentration
c*
parameter, introduced in Eq. (20)
D
diffusion coefficient
Ds
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)
Jj
capillary-tissue diffusion flux
k
argument of the complete elliptic integral
K
complete elliptic integral
K′, K′
defined by Eq. (27)
L, Lj
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
qj
defined in Eq. (13)
qj
defined in Eq. (22)
Qj
capillary volumetric blood flow rate
r, rj
radius of the capillary
r*
characteristic intercapillary distance
rt,rt,rt
Krogh cylinder radius
Rj
defined in Eqs. (60), (64), and (65)
sn, cn, dn
Jacobian elliptic functions
S, Sb
solubility coefficients
u
defined in Eq. (16)
u1, u2
defined in Eq. (42)
υ
defined in Eq. (49)
x, y, s
rectangular coordinates
X, Y, Z
defined in Eq. (30)
z
complex variable

GREEK SYMBOLS

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

References

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