Search tips
Search criteria 


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

Analysis of Capillary-Tissue Diffusion in Multicapillary Systems


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

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

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


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


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


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


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


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,


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


On the capillary-tissue interface Γj,


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



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


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




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


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


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


Substitution of


into (12) leads to the Poisson equation


with the Neumann-type boundary conditions



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]


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




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


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


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


and the parameter k is a root of the transcendental equation


In the following equations the standard notation used will be


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,


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


correspond to each other.

Introducing the notation


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


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


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


has been used.

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




we can rewrite the solution (20) in the form


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


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


where [zeta] is the complex conjugate of ζ.

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


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




Therefore the concentration c is given by


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


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


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.


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










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


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


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:
    Substitution of (60) into (45) gives to the lowest order with respect to the small parameter rj/a
  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
  3. zj = ±a, ±a + 2bi. In this case, the quarter of the circle Γj is mapped approximately into semi-circle Γj in the ζ-plane:
    The average concentration can be calculated using the relationships


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.


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


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

The conservation of mass of the substance (21) yields


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



where the area of the rectangle ω is denoted by


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







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


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


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:


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.

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


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


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



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


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


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



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



we rewrite the equations (2) in the form


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





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



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

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



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


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

The concentration at the venous end of the capillary equals


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


and using (94) and (95) we obtain


The concentration at the venous ends of the capillaries is


[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


for the first and second capillary, respectively.

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


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:


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


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


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


[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:

  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.


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.


2a, 2b
dimensions of the fundamental rectangle
A, B
constants in (91), (92) and (94), (95)
concentration in the tissue
c0, cj0
inlet concentration
cA, cB
defined by (79) and (81)
outlet concentration
intracapillary concentration
parameter, introduced in Eq. (20)
diffusion coefficient
diffusion coefficient in the s-direction
function in Eq. (20)
defined in Eq. (38)
consumption/production rate
defined in Eq. (13)
capillary-tissue diffusion flux
argument of the complete elliptic integral
complete elliptic integral
K′, K′
defined by Eq. (27)
L, Lj
capillary length
characteristic capillary length
number of capillaries in the fundamental rectangle
Neumann function, (Eq. 20)
defined in Eq. (39)
o, O
orders of magnitude
defined in Eq. (13)
defined in Eq. (22)
capillary volumetric blood flow rate
r, rj
radius of the capillary
characteristic intercapillary distance
Krogh cylinder radius
defined in Eqs. (60), (64), and (65)
sn, cn, dn
Jacobian elliptic functions
S, Sb
solubility coefficients
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)
complex variable


defined by (89)
boundary of the rectangle ω
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, [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


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.