MATHEMATICAL ANALYSIS OF BLOOD FLOW AND OXYGEN

TRANSPORT IN MICROCIRCULATION

By

Jun Liu

A DISSERTATION

Submitted to

Michigan State University

in partial fulfillment of the requirements

for the degree of

Applied Mathematics - Doctor of Philosophy

2020

ABSTRACT

MATHEMATICAL ANALYSIS OF BLOOD FLOW AND OXYGEN TRANSPORT IN

MICROCIRCULATION

By

Jun Liu

This dissertation presents a mathematical analysis of oxygen convection, diffusion and

consumption in a uniform and multi-capillary muscle tissue region. Capillaries are tiny ves-

sels connecting arterioles with venules and forming networks through out the body. The

diffusion of substrate, such as oxygen, from the microcirculation through capillaries and

the its consumption by tissue cells are basic in human physiology. Blockage or shortage of

blood in one or more capillaries due to thrombosis or systemic hypoperfusion could lead to

pathological conditions such as stroke. However, transient process of capillary-tissue oxy-

gen transport is poorly understood. The object of this thesis is to develop a mathematical

and computationally efficient model that estimates transient oxygen concentration in three

dimensional striated muscles (e.g. cardiac muscle), which is a four dimensional unsteady

convection-diffusion-consumption moving boundary problem. In particular, we aim to sim-

ulate the consequential stages of pathological conditions such as hyperoxia or hypoxia, by

determining oxygen concentration levels and its transport in order to better understand and

predict adverse health effects. Our computational method is based on a compartmental

modeling concept. It is simple, efficient and can be applied to boundaries of arbitrary shape.

In Chapter 1 the problem is introduced. In Chapter 2 the discrete compartmental modeling

is described. In Chapter 3 the discrete compartmental method is developed for 1D and 2D

problems. In Chapter 4 the method is extended to 3D problems and is used to simulate the

development and the recovery process of anoxia. In Chapter 5 theoretical analysis of the

method is given. Chapter 6 contains discussions, conclusion and possible future work.

Copyright by
JUN LIU
2020

I dedicate this dissertation to my parents, Shaojun Chen and Dr. Tieliang Liu, my husband
Dr. Kaixu Yang, and many friends, who all have been supporting me during this hard time.

v

ACKNOWLEDGMENTS

Here, I would like to express my deepest gratitude to my advisor Dr. Chang Yi Wang for his

guidance towards my studies and researches. Dr. Wang is extremely kind and knowledgeable.

He always provide constructive insights and suggestions that help me make great progress.

Without Dr. Wang’s guidance, I wouldn’t be able to have such a profound understanding

in this gorgeous field.

I would also like to extend my sincere appreciation to my dissertation committee mem-

bers, Dr. Milan Miklavcic, Dr. Moxun Tang and Dr. Guowei Wei. Their comments and

suggestions are extremely beneficial for my research.

I am also grateful to the help I obtained from all the professors in the Department of

Mathematics, both academically and in other aspects. During my Ph.D. life, I learn a lot

from the courses offered by these excellent professors, and I’m able to apply the knowledge

I learn to my research. The professors who I have worked with as a teaching assistant also

impress me a lot by the way they deal with different kinds of difficulties.

During my Ph.D. life, I made a lot of friends and I never feel lonely because of them.

Many thanks to my friends at Michigan State University. They either have graduated and

become successful faculty members or scientists in big companies, or are still working hard

pursuing their Ph.D. degree. I sincerely wish all of them a wonderful future.

Last but not least, I would like to express my sincere thanks to my parents for their

support and concerns, as well as my husband for his accompany and encouragements.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

xi

KEY TO SYMBOLS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv

KEY TO ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi

Chapter 1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1 Biological Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Objectives and Plan of the Current Research . . . . . . . . . . . . . . . . . .

Chapter 2 Discrete Compartment Method (DCM) . . . . . . . . . . . . . .
2.1
Idea of the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Multi-capillary Oxygen Transport in Striated Muscle . . . . . . . . . . . . .
2.2.1 The General Problem . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Discrete Compartmental Systems . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Basic Terms and Definitions . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Graphic Representation . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Application of the DCM . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Discussion and Comparison with Previous Models . . . . . . . . . . . . . . .

Chapter 3 Unsteady Diffusion-Consumption-Convection . . . . . . . . . .
3.1 One Dimensional Unsteady Diffusion-Consumption . . . . . . . . . . . . . .
3.1.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1.2 The Crank-Gupta Problem . . . . . . . . . . . . . . . . . . . . . . . .
3.1.3 Three Capillaries of Uneven Strength . . . . . . . . . . . . . . . . . .
3.2 Two Dimensional Unsteady Diffusion-Consumption . . . . . . . . . . . . . .
3.2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Two capillaries of uneven strength in a circular domain . . . . . . . .
3.2.3 Nine capillaries of uneven strength in a rectangular domain . . . . . .
3.3 Two Dimensional Unsteady Diffusion-Consumption-Convection . . . . . . . .
3.3.1 Formulation: The Space and Time Step Matching Technique . . . . .
3.3.2 A Simplified Model: Unsteady

Convection-Permeation-Consumption . . . . . . . . . . . . . . . . . .
3.3.3 Unsteady Convection-Diffusion-Consumption . . . . . . . . . . . . . .

Chapter 4 The Three Dimensional Modeling . . . . . . . . . . . . . . . . . .
4.1 Three-Dimensional Unsteady Diffusion-Consumption-Convection . . . . . . .

1
3
6
14

16
17
18
18
19
20
21
25
27
31

36
36
38
47
54
59
61
70
74
77
80

82
88

96
96

vii

4.1.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
99
4.1.2 Three capillaries of uneven strength in a cuboid domain . . . . . . . . 105
4.1.3 The Krogh Tissue Cylinder
. . . . . . . . . . . . . . . . . . . . . . . 111
4.2 Simulating 3D Multi-Capillary Oxygen Transport using DCM . . . . . . . . 113
4.2.1 Anoxia Following Capillary Occlusion . . . . . . . . . . . . . . . . . . 117
4.2.1.1 Recovery Through Eliminating Capillary Blockage
. . . . . 121
4.2.1.2 Recovery Through Decreasing Consumption Rate . . . . . . 123
4.2.1.3 Recovery Through Increasing Capillary Blood Velocity . . . 125
4.2.1.4 Recovery Through Increasing Oxygen Concentration in Blood 127
4.2.2 Anoxia Following Increased Tissue Oxygen Consumption . . . . . . . 129
4.2.2.1 Recovery Through Decreasing the Consumption Rate . . . . 131
4.2.2.2 Recovery Through Increasing Capillary Blood Velocity . . . 133
4.2.2.3 Recovery Through Increasing Oxygen Concentration in Blood 136

Chapter 5 The Analytic Theory of the Discrete Compartmental Method 138
5.1 Numerical Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
5.2 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

Chapter 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
6.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
6.2 Future Works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164

viii

LIST OF TABLES

Table 3.1: 1D schematic diagram of different type of boundaries and their labels. La-
bel “BC21” = the left boundary compartment and label “BC12” = the right
boundary compartment. “CC1” = capillary compartment of the (m = 1)
capillary. “IC” = inactivated compartment.
. . . . . . . . . . . . . . . .

Table 3.2: Index table of different types of compartments. The corresponding relation-
ships of the compartment label and the coefficients is given. Label “CCm”
represents the mth capillary compartments, label “TC” represents the inte-
rior tissue compartments, label “BC21” represents the left boundary com-
partment and label “BC12” represents the right boundary compartment.
K indicates the value of the consumption term.
. . . . . . . . . . . . . . .

Table 3.3: Compartmental map for the Crank-Gupta problem.
. . . . . . . . . . . .
Table 3.4: Comparison for C × 106 at different times. On each time level, the first row
shows values taken from Hansen & Hougaard [29], the second row shows
values from Crank and Gupta [12] with ∆x = 0.05, the third row gives
values obtained from DCM with ∆x = 0.05 and the forth row the DCM
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
with ∆x = 0.0125.
Table 3.5: Comparison for C × 106 at different times. On each time level, the first row
shows values taken from Hansen & Hougaard [29], the second row shows
values from Crank and Gupta [12] with ∆x = 0.05, the third row gives
values obtained from DCM with ∆x = 0.05 and the forth row the DCM
with ∆x = 0.0125.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

47

50

51

52

Table 3.6: Comparison of results at the sealed surface.

. . . . . . . . . . . . . . . . .

53

Table 3.7: Comparison of the oxygen depletion time with previous works.

. . . . . .

53

Table 3.8: Comparison of oxygen depletion time (t*) in different mesh size using DCM. 54

Table 3.9: Schematic diagram of label changing following capillary blockage in the
1D example. Label “CCM” = capillary compartment corresponding to the
oxygen concentration of the mth capillary. “CC0” = capillary compartment
. . . . . . . . . . . . . . . . . . . .
corresponding to a blocked capillary.

Table 3.10: Discrete compartmental map of 2D unsteady diffusion-consumption prob-
lem over rectangular shaped domain. Label “TC” = tissue compartment;
label “BC##”=boundary compartment of type ##; label “IC” = inacti-
vated compartment.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

67

ix

Table 3.11: Discrete compartmental map of 2D unsteady diffusion-consumption prob-
lem over irregular shaped domain. Label “TC” = tissue compartment; label
“BC##”=boundary compartment of type ##; label “IC” = inactivated
compartment.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Table 3.12: Locations and flux strength of nine capillaries used in the example. . . . .

75

Table 3.13: Discrete compartmental map of unsteady convection-permeation-consumption

problem. Label “CC”=capillary compartment;
label “CCA”=capillary
compartment on the arterial end; label “CCV”=capillary compartment on
the venous end; label “STC” = single layer tissue compartment; label “IC”
= inactivated compartment. . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3.14: Discrete compartmental map of 2D unsteady convection-diffusion-consumption

problem. Label “CC”=capillary compartment;
label “CCA”=capillary
compartment on the arterial end; label “CCV”=capillary compartment on
the venous end; label “TC” = tissue compartment; label “CTC” = capillary
neighboring tissue compartment; label “BC##”=boundary compartment
. . . . . . . . . . . .
of type ##; label “IC” = inactivated compartment.

85

89

Table 4.1: Locations and plasma oxygen concentration on the arterial end of three

randomly distributed capillaries used in the example. . . . . . . . . . . . . 105

Table 4.2: Time to fully perfused (tf ) and to steady state (ts) corresponding to dif-

ferent valid choices of spacial step size (∆x) . . . . . . . . . . . . . . . . . 107
Table 4.3: Convergence for concentration value (×106) on a random point (0.243, 0.525, 0.725)

corresponding to an interior tissue compartment for different choices of
valid spacial step size (∆x) . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Table 4.4: Biophysical parameters used for the simulation of oxygen transport in

multi-capillary tissue bed. R(cid:48)
functional capillary density; v(cid:48)
gen concentration; D(cid:48), tissue oxygen diffusivity; κ(cid:48)
consumption; κ(cid:48)

c, capillary radius; L(cid:48), capillary length; FCD,
b, normal RBC velocity; Ca, arterial oxy-
b, basal tissue oxygen

max, maximal tissue oxygen consumption.

. . . . . . . . . 113

Table 4.5: Locations and plasma oxygen concentration on the arterial end of eight

capillaries.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

x

LIST OF FIGURES

Figure 1.1: Top 10 global causes of deaths, 2016.

. . . . . . . . . . . . . . . . . . . .

5

Figure 1.2:

(a) A single Krogh Tissue Cylinder. Estimation of P O2 within a circular
cylinder of tissue surrounding a capillary of radius Rc; r is the distance
from the capillary center; Rt denotes the cylinder radius where oxygen flux
becomes zero. (b) P O2 at capillary declines monotonically both around
and within the fibre. (c) Krogh’s view of circular tissue cylinder stacking,
where tissue supply voids are inevitable. This arrangement is also known
as concurrent flow model.
. . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.3:

(a) Oxygen concentration curves of seven unevenly distributed capillaries
of unequal strength over a circular tissue region. Total flux given equals
to total area of the circular region times the tissue oxygen consumption
rate. At the point on upper right boundary, the concentration is zero. (b)
. . . . . . . . . . . . . . . .
The corresponding capillary supply regions.

Figure 3.1: Compartmentalization, map and labels of a 1D discrete compartmental
model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

11

39

Figure 3.2: 1D compartmental model

. . . . . . . . . . . . . . . . . . . . . . . . . .

40

Figure 3.3: Calculation of capillary oxygen concentration after blocking the source.
. . . .

The three darker cells denote contiguous capillary compartments.

Figure 3.4: Oxygen diffusion from three capillaries of different strengths into a con-
suming region at t=0.001, t=0.01 and t=0.1, compared with the analytic
steady state solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.5: Concentration profile at t=0.001, t=0.02 and t=0.2 times after a blockage
occurred on middle capillary, compared with the analytic steady state
solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.6: Concentration profile at t=0.002, t=0.02 and t=0.03 times after a blockage
. . . . . . . . . . . . . . . . . . . . . . .

occurred on the right capillary.

Figure 3.7:

schematic diagram of 2D planar capillary-tissue region and blockage on
one capillary. Shaded area in part (b) represents the possible anoxic areas
due to blockage on the right capillary.
. . . . . . . . . . . . . . . . . . .

43

55

56

57

59

Figure 3.8: Example of 2D capillary shapes. . . . . . . . . . . . . . . . . . . . . . . .

62

xi

Figure 3.9: The 2D compartmental model.

. . . . . . . . . . . . . . . . . . . . . . .

63

Figure 3.10: Transient process of oxygen diffusion from two capillaries of different
strengths into a circular consuming region. Dark regions show the zero
oxygen area. The vertical and horizontal axes shows the compartment
numbers used in each direction.
. . . . . . . . . . . . . . . . . . . . . . .

Figure 3.11: Transient process of oxygen diffusion from two capillaries of different
strengths into a circular consuming region. Dark regions show the zero
oxygen area. The vertical and horizontal axes shows the compartment
numbers used in each direction.
. . . . . . . . . . . . . . . . . . . . . . .

Figure 3.12: Concentration profile of two unequal strength capillaries in a circular do-
main, comparing to the analytic solution[68]. . . . . . . . . . . . . . . . .

Figure 3.13: Flux and supply boundary of two capillaries of unequal strength in a
. . . . . . . . .

circular domain, comparing to the analytic solution [68].

Figure 3.14: Oxygen concentration level curves of nine capillaries with different strengths

at steady state. The vertical and horizontal axes shows the compartment
numbers used in each direction.
. . . . . . . . . . . . . . . . . . . . . . .

Figure 3.15: Schematic diagram for 2D capillary-tissue exchange with convection, dif-
fusion and consumption with nondimensionalized parameters. vb is the
blood velocity; D = 1 is the normalized diffusion coefficient; κ is the tis-
sue consumption coefficient; L is the length of the tissue region along the
z-axis; W is the thickness of the tissue region along the x-axis; oxygen
concentration on the arterial end of the capillary is normalized by itself
and therefore equals to 1.
. . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.16: Lagrangian sliding fluid element method. With each ∆t, oxygen in the
capillary compartment slide instantaneously one segment downstream, af-
. . . . . . . . .
ter which the diffusion and consumption are computed.

Figure 3.17: Compartment model for capillary-tissue exchange in well-perfused organ.

Figure 3.18: Oxygen convection into a consuming region using a two region model.
Red dots represent oxygen concentration flowing in capillary, blue dots
represent oxygen concentration in adjacent tissue region with constant
. . . .
consumption and solid lines are analytic solutions at steady state.

72

73

74

74

76

78

81

82

87

Figure 3.19: Compartment model for capillary-tissue exchange with diffusion.

. . . .

88

xii

Figure 3.20: Concentration curves of a single capillary with blood velocity at t=0.1,
t=0.25, t=1, t=2.5, t=5, t=7.5 and t=10. The vertical and horizontal
axes shows the compartment numbers used in each direction. . . . . . . .

93

Figure 3.21: Concentration profile on the capillary boundary and the tissue boundary,
compared with the analytic solution. Circles are numerical solution of
oxygen concentration on the tissue boundary. Triangles are numerical so-
lution of oxygen concentration on the capillary wall. Solid lines represents
the oxygen concentration profile obtained from the analytic solution. Hor-
izontal axis shows the normalized capillary length and vertical axis show
the normalized concentration value. . . . . . . . . . . . . . . . . . . . . .

94

Figure 4.1: A 3D compartment for deriving the difference equations for DCM. . . . . 100

Figure 4.2: Concentration profile of three capillaries of uneven strength in a cube

domain at t = 0.1, t = 0.2, t = 0.5 and t = 2.2 (steady state).

. . . . . . 106

Figure 4.3: Convergence of the 3D algorithm on interior tissue compartments for dif-
ferent choices of valid spacial step size (∆x), at time t = 0.1. The vertical
axis represents the concentration value(×106) of the compartment and the
horizontal axis is the inverse of spacial step size (∆x). Vertical black lines
divides the plane into sub-regions due to the sliding period corresponding
to the compartment size. Short horizontal lines are the average of values
in each rectangular region and the long horizontal line is calculated using
the uniformly minimum variance unbiased estimator from the end of the
smallest ∆x obtained.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

Figure 4.4: The Krogh tissue cylinder . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 4.5: Concentration curves of the Krogh tissue cylinder at (C = 0.6, 0.7, 0.8, 0.9)
at steady state from the numerical solution. The vertical and horizontal
axes shows the compartment numbers used in each direction. . . . . . . . 112

Figure 4.6: Concentration profile of the Krogh tissue cylinder on the capillary bound-
ary and the tissue boundary. Circles are numerical solution of oxygen
concentration on the capillary wall. Triangles are numerical solution of
oxygen concentration on the tissue boundary. Solid lines represents the
oxygen concentration profile obtained from the analytic solution. Hori-
zontal axis shows the normalized capillary length and vertical axis show
the normalized concentration value.

. . . . . . . . . . . . . . . . . . . . 112

xiii

Figure 4.7: Oxygen concentration profile of eight capillaries of uneven strength and
distribution in a cuboid domain at steady state. Biophysical parameters
collected in Table 4.4 is used to simulate the actual biological process. Un-
der capillary blood velocity 1150µm/s and basal tissue oxygen consump-
tion rate 1.67 × 10−4 (ml O2)(ml tissue)−1 s−1, oxygen concentration in
this capillary-tissue system reaches the steady state at a normalized time
of t = 0.6.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure 4.8: The development of oxygen deficit by blocking one capillary closer to the
middle of the tissue domain, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic
areas are shown in white for contrast. . . . . . . . . . . . . . . . . . . . . 118

Figure 4.9: The development of oxygen deficit by blocking one capillary near the tissue
domain boundary, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic areas are
shown in white for contrast.

. . . . . . . . . . . . . . . . . . . . . . . . . 120

Figure 4.10: Recovery process of tissue oxygen distribution through unclog the blocked
capillary, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown in
white for contrast.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Figure 4.11: Recovery process of tissue oxygen distribution following a reduction of the
tissue consumption, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are
shown in white for contrast.

. . . . . . . . . . . . . . . . . . . . . . . . . 124

Figure 4.12: Recovery process following an increment in the capillary blood velocity,
at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown in white for
contrast. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

Figure 4.13: Recovery process of tissue oxygen distribution following increasing oxygen
concentration in the arterial blood, at time t = 0.001, 0.05, 0.1 and 0.2.
Anoxic areas are shown in white for contrast.

. . . . . . . . . . . . . . . 128

Figure 4.14: The development of oxygen deficit by increasing the tissue oxygen con-
sumption rate, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic areas are shown
in white for contrast. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

Figure 4.15: Recovery process of tissue oxygen distribution following reducing the tis-
sue consumption rate to original state, at time t = 0.001, 0.05, 0.1 and
0.2. Anoxic areas are shown in white for contrast. . . . . . . . . . . . . . 132

Figure 4.16: Recovery process of tissue oxygen distribution following increasing the
capillary blood velocity to 1500µm/s, at time t = 0.001, 0.01, 0.05 and
0.1. Anoxic areas are shown in white for contrast. . . . . . . . . . . . . . 134

xiv

Figure 4.17: Recovery process of tissue oxygen distribution following doubling the cap-
illary blood velocity, at time t = 0.001, 0.01, 0.05 and 0.1. Anoxic areas
are shown in white for contrast. . . . . . . . . . . . . . . . . . . . . . . . 135

Figure 4.18: Recovery process of tissue oxygen distribution following increasing oxygen
concentration in the arterial blood, at time t = 0.001, 0.01, 0.05 and 0.1.
Anoxic areas are shown in white for contrast.

. . . . . . . . . . . . . . . 137

xv

KEY TO SYMBOLS

Normalized capillary radius

Normalized tissue radius of the Krogh-Cylinder Model

Normalized capillary length

Normalized width of the capillary-tissue region

Normalized height of the capillary-tissue region

Normalized capillary-tissue domain

Normalized domain of the mth capillary

Normalized total capillary domain

Normalized total tissue domain

The ith normalized anoxia tissue region

Normalized total anoxia tissue region

Number of capillaries in the capillary-tissue region

Normalized capillary blood velocity

Normalized arterial oxygen concentration

Normalized tissue oxygen diffusion coefficient

Normalized tissue oxygen consumption coefficient

Normalized oxygen concentration

Normalized oxygen concentration in the capillary region

xvi

Rc

Rt

L

W

H

Ω

Ωm

Ωc

Ωt

Oi

O

M

vb

Ca

D

κ

C

Cc

Ct

Σnb

p

b

A

V

N

J

∆x

∆t

∆t+

∆t−
t∗

tf

ts

nCC

q

P

Normalized oxygen concentration in the tissue region

Operator representing the summation of compartment neighborhoods

The sliding period of capillary compartments

The sliding frequency of capillary compartments

The cross section area of a compartment

The volume of a compartment

Mass of oxygen molecules

Oxygen flux between adjacent compartments

The spacial step size

The time step size

The beginning of the time interval

The end of the time interval

The oxygen depletion time of the tissue region

The tissue fully perfused time of the tissue region

The time spent to reach steady state from transient state

Number of compartments representing the cross sectional area of a capillary

Capillary oxygen flux strength

Oxygen permeability through the capillary membrane wall

P O2

Oxygen partial pressure

xvii

i, j, k, m, n

Numerical indices

Ω(cid:48)

Ω(cid:48)

m

Ω(cid:48)

c

Ω(cid:48)

t

O(cid:48)

i

O

R(cid:48)

c

L(cid:48)

W(cid:48)

H(cid:48)

v(cid:48)

b

C(cid:48)

a

D(cid:48)

κ(cid:48)

C(cid:48)

κ(cid:48)

b

κ(cid:48)

max

Capillary-tissue domain (1D: µm, 2D: µm2, 3D: µm3)

Domain of the mth capillary (1D: µm, 2D: µm2, 3D: µm3)

Total capillary domain (1D: µm, 2D: µm2, 3D: µm3)

Total tissue domain (1D: µm, 2D: µm2, 3D: µm3)

The ith anoxia tissue region (1D: µm, 2D: µm2, 3D: µm3)

Total anoxia tissue region (1D: µm, 2D: µm2, 3D: µm3)

Capillary radius (µm)

Capillary length (µm)

The width of the capillary-tissue region (µm)

The height of the capillary-tissue region (µm)

Capillary blood velocity (µm/s)

Arterial oxygen concentration ((ml O2)/(ml blood))

Tissue oxygen diffusion coefficient (µm2/s)
Tissue oxygen consumption coefficient ((ml O2)(ml tissue)−1 s−1)

Oxygen concentration ((ml O2)/(ml blood))
Basal tissue oxygen consumption coefficient ((ml O2)(ml tissue)−1 s−1)
Maximal tissue oxygen consumption coefficient ((ml O2)(ml tissue)−1 s−1)

xviii

KEY TO ABBREVIATIONS

DCM

Discrete Compartmental Method

SEM

Sliding Element Method

CC

TC

BC

IC

MI

Capillary Compartment

Tissue Compartment

Boundary Compartment

Inactive Compartment

Myocardial Infarction

AMI

Acute Myocardial Infarction

FCD

Functional Capillary Density

FVM

Finite Volume Method

PDE

Partial Differential Equations

RBC

Red Blood Cell

xix

Chapter 1

Introduction

Our body constantly needs oxygen to stay alive. Oxygen were delivered and consumed

by tissue through blood flow inside the arteries. However, the arteries have thick vessel

walls almost impermeable to oxygen. The diffusion of oxygen (and nutrients) to tissue

actually occurs in the more permeable connecting capillaries. We are interested in the

possible unsteady change in oxygen concentration inside the capillary-tissue region.

Such a study is important because low oxygen delivery from the blood flow will result in

insufficient oxygen concentration for the muscle tissue. The cause of low oxygen could be due

to disease, low oxygen environment or excertion. In the muscle, low oxygen (hypoxia) may

cause soreness or pain. In the heart muscle, it may cause a heart attack and in the brain, a

stroke. Up to now, it is still difficult to measure the distribution of oxygen in a living body.

The theory or its variations remains constrained by a simple mathematical model developed

a century ago, with many deficiencies. The present research is to develop a more advanced

theoretical model.

Traditional mathematical modeling of biological processes is as follows. The problem is

described by partial differential equations which are then solved by numerical methods such

as the finite difference methods or finite element methods. For our problem, the modeling is

difficult because

1. The geometry is random and three dimensional due to natural dividing organism

1

2. The process is unsteady. This means time is an extra variable

3. The process involve diffusion - oxygen diffuses into the tissue area

4. The process also involve convection - oxygen is transported by blood flow in the cap-

illary

5. Consumption - the tissue will constantly consume the oxygen

6. Moving boundary problem - if certain regions is devoid of oxygen supply, the concen-

tration will have a moving boundary

Thus, it is a three-dimensional transient convection-diffusion-consumption moving boundary

problem. It would be a severe challenge to solve this problem with PDE or even numerical

method constructed from PDE.

This thesis presents a new method, names discrete compartmental method (DCM). The

region is divided into many small cubes, called compartments, each represents part of a

tissue or capillary based on their location. Based on their type, they are assigned of specific

difference equations, constructed directly from the physical laws, to update their value in

the next time step. By attaching the cubes together, we are able to simulate the oxygen

concentration change of the entire domain.

In addition, using a sliding element method (SEM) together with the time matching

method, the time derivative is naturally “eliminated” and the dimension of the problem was

reduced from 4D to 3D. Since the method is constructed based on the real biological process,

it provides accurate results.

2

1.1 Biological Background

The understanding of oxygen transport to tissue is basic and crucial in human physiology.

A sufficient supply of oxygen to the cells is critical to any form of energy production [52].

In mammals, oxygen is extracted from the atmospheric air in the lungs, and transported

by erythrocytes and dissolved in plasma through the circulation to the tissue. However,

actual oxygen supply from blood to tissue only occurs in the microcirculation, where oxygen

diffuses from capillaries into the tissue and is consumed there within the mitochondria.

1. Microcirculation

The microcirculation is the circulation of the blood in the smallest blood vessels, includ-

ing terminal arterioles, metarterioles, capillaries, and venules. Within organ tissues,

arterioles carry oxygenated blood to the capillaries, and blood flows out of the capillar-

ies through venules into veins [66]. In the microcirculation, substances such as oxygen

and nutrients are transported by flow convection [48]. During which process, the sub-

strate leaves the plasma, crosses the vessel wall, and diffuses into tissue. Within the

tissue, the substrate diffuses in a three dimensional region across many different his-

tological structures and fluid spaces. Eventually, the substrate reaches the metabolic

site in the cell and is consumed there [53].

2. Oxygen Transport in Striated Muscles

Capillaries are the smallest blood vessels connecting arterioles with venues and forming

networks through the human body. The size of capillaries are approximately 5-8 µm

in diameter [60]. There are approximately seven percent of the body’s blood is in

the capillaries which continuously exchange substances with the surrounding tissues,

3

named capillary exchange [66]. In striated muscles, capillaries are long and parallel

[38]. Glucose, amino acids, oxygen and other molecules exit capillaries by diffusion to

reach the tissues and is consumed there [66]. This process depends on the difference

of concentrate between the interstitial fluid and blood, with molecules moving in the

direction of declining concentration of that substrate in solution [35].

3. Significant Problem: Hypoxia

Any form of energy production depends on a continuous supply of oxygen to the living

cells. Blockage of blood flow in capillaries due to thrombosis, embolism or trauma will

cause ischemia, leading to oxygen deficiency in tissue (hypoxia), which, if untreated,

leads to absence of oxygen (anoxia), and eventually tissue death (necrosis). According

to the recent result from the World Health Organization (WHO) [49], ischemic heart

disease such as Myocardial infraction (MI) or acute myocardial infarction (AMI), com-

monly known as a heart attack, due to insufficient blood flow at part of the heart

causing damage to the heart muscle have remained the top leading cause of death in

the world during the past decade.

4

Figure 1.1: Top 10 global causes of deaths, 2016.

Ischemic heart disease occurs when blood flow is insufficient in a part of the heart caus-

ing tissue death to the heart muscle. The most common symptom is chest pain, termed

angina pectoris, due to the sudden ischemia of the heart muscle. The discomfort may

quickly radiate into the shoulder, arm, back, neck, or jaw. Loss of consciousness or

sudden death also occurs in MIs. Under such circumstance, time is of crucial impor-

tance. Acetylsalicylic acid is an appropriate immediate treatment for a suspected MI;

however, it do not improve overall recovery. Thus, future research, especially theory

for the occurrence and properties of MI is sorely needed. This thesis attempts to an-

swer the following basic questions, which may help in the understanding, diagnosis and

subsequent treatment of MI.

1. What is the transient oxygen concentration distribution in tissue supplied by a

bunch of capillaries under both oxygen sufficient and insufficient conditions?

2. Where, when and how would hypoxia occur if one or more capillary were blocked?

5

3. Could the neighboring capillaries supply sufficient oxygen to the anoxic region?

How soon and how much?

1.2 Literature Review

Representing the microcirculation process using mathematical modeling, especially on de-

scribing the concentration and consumption of oxygen in muscles, has received a lot of

attention early in the beginning of 1910s and the precursor of this field is August Krogh[38]

[39]. From observing the cross-section of striated muscle, he theorized that the rate of oxy-

gen transport was related to the number and arrangement of capillaries in tissue and to

the oxygen permeability of the capillary walls. Krogh’s work in the early stage was only

based on experimental and statistical results. However, his collaboration with mathemati-

cian Agner Erlang in 1919 opened the gate of mathematically modeling the microcirculation

[38]. They developed the first theoretical analysis of oxygen supply within muscle tissue,

which predicts the oxygen tension (partial pressure) distribution around a single capillary.

Krogh’s model assumed each capillary could be regarded as parallel to the muscle tissues

and supplying an independent concentric tissue region surrounding that capillary. He also

determined the average radius of such hypothetical tissue cylinder by counting the number of

capillaries in a certain cross-section area. With the assistance of the mathematician Erlang,

the mathematical description of such an ideal physical geometric arrangement, shown in Fig.

1.2, predicting the oxygen tension along the tissue cylinder, was known as the Krogh-Erlang

equation

p (r) = p (Rc) − M0
4k

−(cid:16)

r2 − R2

c

(cid:17)(cid:21)

,

(cid:20)

R2

t log

r2
R2
c

6

(1.1)

Figure 1.2: (a) A single Krogh Tissue Cylinder. Estimation of P O2 within a circular cylinder
of tissue surrounding a capillary of radius Rc; r is the distance from the capillary center; Rt
denotes the cylinder radius where oxygen flux becomes zero. (b) P O2 at capillary declines
monotonically both around and within the fibre. (c) Krogh’s view of circular tissue cylin-
der stacking, where tissue supply voids are inevitable. This arrangement is also known as
concurrent flow model.

where Rt and Rc are the tissue and capillary radii with Rc < r < Rt, k is Krogh’s

oxygen (O2) diffusion coefficient in tissue, and M0 is a constant tissue demand for O2.

Other assumptions of the Krogh Tissue Cylinder model includes:

1. Tissue O2 consumption is constant and uniform.

2. Tissue partial pressure or O2 concentration (P O2) at the capillary wall equals average

capillary P O2.

3. Tissue O2 solubility and diffusivity are uniform.

4. Axial (or longitudinal) diffusion of O2 is not significant in comparison to radial diffu-

sion.

5. Micro-vascular O2 transport phenomena is steady.

6. All capillaries are parallel, unbranched, and equally spaced.

7. All capillaries receive equal convective O2 supply.

Based on Krogh’s initial steady state model, A. V. Hill [30] derived unsteady state equa-

tions in 1928. G. Thews in 1950 [64] and F. W. J. Roughton in 1952 [56], developed steady

7

and unsteady state equations for Krogh tissue cylinder. In 1957, Seymour S. Kety [36] de-

veloped the first model including intra-capillary conditions. He proposed for the first time

that concentration of oxygen in capillaries should drop gradually from the arterial to the

venous side. Jackob Blum [6], in 1960, gives another innovative contribution by adding the

permeability of capillary wall and the metabolic rate in the tissue space to the Krogh tissue

cylinder. However, Blum did not considered bound substrate such as O2 in his work. His

result is not adaptive to oxygen but can be applied to diffusible substrate such as glucose.

The increasing availability of computers by the 1960s encouraged many investigators

started working not only on generalizing and completing the description of Krogh cylin-

der but also on developing new geometric models as well as exploring adaptive numerical

techniques. Gonzalez-Fernandez[26], Filtcher [16] [17] and Reneau[55] in the period 1965 to

1975 studied the unsteady Krogh model. Reneau’s work in 1966[55] introduced the finite

difference equations to treat the nonlinear unsteady system for the first time. Salthe and

Wang [59] studied the change of oxygen concentration following a capillary occlusion using

the method of matched asymptotic expansions. More detailed presentations and discussions

of the Krogh cylinder can be found in Middleman [44] and Lightfoot [42] and the reviews of

Fletcher in 1977 [18] and Popel in 1989 [52].

Over more than half a century, the original framework build by Krogh, the Krogh tissue

cylinder, has always been the foundation of studying and modeling oxygen transport. How-

ever, the concept of a“lethal corner” (i.e. tissue area close to the cylinder boundry on the

venuous end) has been troublesome to many physiologists, since they are not able to isolate

and identify a particular capillary and its associated cylinder of tissue, and then make the

measurements in such tissue space. Although Silver [61], Whalen[70] and several others have

measured low values of oxygen tension in the tissues of small animals, yet no experimental

8

measurements could disprove or confirm the single tissue cylinder model.

For the purpose of determine if a single capillary model is valid for the predictions of the

oxygen distribution, several attempts to model large capillary networks were made in the

early 1970s. Hutten, Thews, and Vaupel [31] modeled considerably large capillary networks

with an electronic analog model and concluded that single capillary may supply larger or

smaller cross-sectional tissue area then that predicted from a single capillary model. Levit

[41], using the method of Johnson [34], considered diffusional and time dependent interaction

of adjacent capillaries. He argues that capillaries interactions are potentially important in

determining the true venous outflow.

Most of the early modeling were based on Krogh’s finding and assumed that the mi-

crocirculation was a passive network where blood is supplied by tube to a nonintereacting

volume of tissue. However, experimental works in the late 1970s challenged this traditional

view. It was observed that the microcirculatory environment is a much more complex locally

regulated intereacting system. Goresky and Goldsmith [27] suggests that diffusional bypass

may happen in cardiac muscle due to the diffusional interaction between adjacent groups

of capillaries. Moreover, Goresky’s study and those of Fable [15], Weiss [69] and others

observed that in cardiac muscle, the blood flow and metabolic function are synchronized

with the contractile cycle. Hence, both spatial and temporal relationships must be consid-

ered for any model attempting to describe substrate supply in cardiac muscle. For skeletal

muscle, Groom, Plyley and their associates [51] found that a capillary unites with branches

from an adjoining parallel network rather than rejoin its generating branch. This group’s

observation once again reveals the importance of capillary interaction to the determination

of venous outflow. As a result, most assumptions for the Krogh model fails to be true and

new adaptive models and numerical methods once again became necessary.

9

With the development of modern technology, several physiology observations have became

definite facts and widely accepted. They are mainly about the heterogeneity of microvascu-

lature, such as 1) capillaries, even those in the skeletal muscle, are not uniformly distributed,

2) the size of capillaries are not quite the same, 3) the O2 capacity of each capillary is differ-

ent, and 4) stronger O2 capillary capacity not only support adjacent area but compensates

neighboring regions. With these new understandings on microvessel’s heterogeneity, many of

the previous assumptions made on developing mathematical model of blood-tissue exchange

may not be applicable.

Considering the uneven capillary distribution and diffusional strength, Wang [68], in 1999,

devised a method to compute the oxygen density supplied by capillaries in a simply connected

circular region and a analytic solution is given. Wang defined the tissue region serviced by

a certain capillary to be the capillary supply region. The analytical result shows that, the

shapes of the supply region depend on both the location and the relative diffusion strength of

the capillaries(Fig. 1.3). Wang argues that a partial blockage of a certain capillary may not

cause anoxia near the blockage, but the neighboring capillaries may compensate by shifting

capillary supply regions.

10

(a)

(b)

Figure 1.3: (a) Oxygen concentration curves of seven unevenly distributed capillaries of
unequal strength over a circular tissue region. Total flux given equals to total area of the
circular region times the tissue oxygen consumption rate. At the point on upper right
boundary, the concentration is zero. (b) The corresponding capillary supply regions.

Following Wang’s work, Go [1] [20] obtained the unsteady solution of oxygen distribution

in circular region. However, he considered negative solutions as the hypoxia region which is

far from accurate for the real case. Sun [62] extended Wang’s method to rectangular shaped

capillary-tissue regions. Again, his results would apply only to steady state situations.

On the order hand, in 1972, Crank and Gupta [12] [11] first discussed one dimensional

moving boundary problem describing the oxygen diffusion and consumption phenomenon in

the body tissue. In Crank and Gupta’s first paper [12], an explicit finite difference method

is used with Lagrange formula interpolation near the moving boundary.

In their second

paper [11], Forward Difference Spline (F.D.S.) Method and Forward Difference Polynomial

(F.D.P.) Method are used to shift the mesh at each step to reduce the irregularities when

crossing a mesh point.

11

Two years later, Hansen and Hougaard [29] derived an integral equation for the position

of the moving boundary and an integral formula for the concentration. The concentration is

estimated asymptotically for small times and calculated numerically for later times. Berger,

Ciment and Rogers [5] considered an equivalent problem as a non-linear equation on the fixed

interval [0, 1] and solved it numerically using finite difference method where the negative

part of the solution is set equal to zero before each step. However, their result has very

low accuracy even with small mesh spacing due to modifying the original problem. Miller,

Morton and Baines [45] used finite difference in time and finite element with adaptive mesh

in space. In the begining of 1980s, Gupta and Kumar [28] considered numerical integration.

Dahmardah & Mayers [13] obtained more accurate results by expanded the solution in a

Fourier series.

Mitchell [46] applied integral method together with a novel transformation to remove the

non-uniform initial profile. Mitchell & Vynnycky [47] connected the Keller box scheme with

the boundary immobilization method and focused on finding the oxygen depletion time in

2015. More detailed discussions about the moving boundary problem can be found in Crank

[10] [9] and Zerroukat [73].

Until recently, it has been admitted that capillaries are not evenly distributed but spa-

tially heterogeneous and the O2 capacities of each capillary are not identical but vary among

different adjacent cells [14]. Moreover, it is also recently being discovered that local capil-

lary capacity in muscle is primarily determined by muscle fiber type and fiber size [8]. For

transient problems, such as abnormal perfusion with more than one capillary, the problem

has not been studied at all. The difficulty is mainly due to the fact that the boundary of

the ischemic area is no longer stationary. For example, assume there is a blockage of one

capillary among many, the ischemic area may be decreasing over time by an increase in flux

12

from the neighboring cells. In addition, unevenly and randomly distributed capillaries lead

to restrictions on model construction which result in large amount of computation time. In

addition, Whiteley et al.

[71] recently examined the significance of axial diffusion in both

the capillary and the tissue on a single capillary model. Using dimensional analysis and nu-

merical methods, they found that axial diffusion in the capillary is negligible except within

a small boundary layer at the inflow region of the order of the capillary radius. However,

the effect of axial diffusion within the tissue was found to be small but significant.

Taking advantage of the increasing computing capability, more and more numerical meth-

ods are investigated to the study of oxygen transport in the micro-circulation. Goldman and

Popel [21] calculated the steady-state blood flow and time-dependent oxygen transport in a

three-dimensional capillary network using a finite difference based method. Using a slightly

modified method, Goldman [22] simulated the capillary network oxygen transport during

ischemia by comparing the result at different steady-states. Using the modified model, Gold-

man et al.

[25] studied the role of heterogeneous capillary spacing and blood flow. Fraser

et al.

[19] compared the generated parallel capillary arrays with three-dimensional recon-

structed capillary networks in modeling oxygen transport. Sun et al. [63] calculated the 2D

steady-state solution of multi-capillaries with uneven location and strength in rectangular

capillary-tissue beds.

Compartmantal method is a classic and commonly used mathematical analysis method on

studying biological processes. Compartments in physiology are composed of sets of intercon-

nected mixing chambers. Each component of the system being considered as homogeneous

instantly mixed, with uniform concentration. A comprehensive introduction of the compart-

mental method can be found in Jacquez [32]. More recent applications and computational

extensions of it can be found in Bassingthwaighte et al. [2].

13

The sliding (fluid) element method is a not so wellknown method which address the

convection. It was first proposed in Bassingthwaighte et al. [3] for the computation of a 2-

barrier, 3-region convection-permeation-diffusion models for blood-tissue exchange problem.

The central idea of the sliding element algorithm is to set the time step equal to the length

step divided by the fluid velocity. This computationally efficient method for convection has

never been extended to a three dimensional model.

1.3 Objectives and Plan of the Current Research

After a century of study, many important questions regarding the transient process of oxygen

concentration distribution in tissue remain open today. The mathematical modelling of the

capillary supply to tissue generates a complex system of coupled mathematical expressions. A

finely detailed geometry of capillary networks requires considerations of both physiologically

interactive spatial effects and computationally attainable algorithm, since it is clear that

interacting microcirculatory regions do not maintain a completely static condition. For

striated muscle system, a mathematical model describing oxygen distribution from a resting

state to a state of maximum activity and maintain functional state under hypoxic conditions

for long periods of time is needed. A computational model which can reflect the adaptive

process of capillary recruitment, increasing and decreasing blood flow, and metabolic activity

would be a significant advance in the state of the art in modeling microcirculation.

The general object of this thesis is to develop a three dimensional mathematical model

with possible extensions and computational efficiency, in order to provide a mathematical

explanation for the transient state of tissue oxygen distribution. In particular, to provide

a new research tool for the study of oxygen transport in striated muscles, the process and

14

degree of tissue hypoxia and, in particular, the forming process and location of an anoxia in

cardiac and skeletal muscles. This new method greatly improves the current frameworks on

both mathematical and physiological applications. Major innovations include:

1. diffusion-consumption problem with moving boundaries,

2. diffusion-consumption problem with multiple unevenly distributed sources and different

strength,

3. three dimensional space diffusion-consumption-convection problem modeling.

4. transient state solutions

In general, the new model is able to solve a three space dimensional, transient, mov-

ing boundary, convection-diffusion-consumption problem with multiple unevenly distributed

sources of different strength. With this new model, certain significant physiological phenom-

ena can be further studied, such as the approximation of local anoxia area caused by sudden

blockage of one or more capillaries and the time for full recovery.

15

Chapter 2

Discrete Compartment Method

(DCM)

In this section we develop the discrete compartment model for the transient diffusion-

consumption-convection problem. First, we define the general problem of unsteady oxygen

exchange problem in multiple capillary domain and discuss the assumptions. Second, we de-

fine the discrete compartment method (DCM), give definitions to the basic terms, introduce

the graphic representation and discuss the basic compartmental difference equations. Lastly,

we will discuss the differences of the present model with other models and give a visual

explanation for demonstrating the essence of the discrete compartmental model. Important

material on numerical accuracy will be discussed in detail later in Chapter 5.

The governing difference equations of discrete compartment methods are similar but

different to those of finite difference methods, and the graphic representation of a discrete

compartmental method are similar but different to those of finite volume methods. However,

the greatest difference is that the discrete compartmental methods are derived from the

basis of physical laws such as Fick’s diffusion law, a starting point that turns out to have

many advantages. Moreover, the discrete compartmental method can easily handle time-

dependent boundary conditions. At the same time, it is computationally efficient even for

four dimensional problems. More advantages and differences to other numerical methods

16

and mathematical models will be discussed in the last section of this chapter.

2.1

Idea of the Model

Consider the whole tissue space partitioned into finite number of contiguous tissue bins.

For sufficient oxygen supply, each bin would consume one normalized unit of oxygen per

normalized unit of time, and some of the bins are designated as capillaries which can deliver

an amount of oxygen per unit of time. If oxygen supply is less than the minimum required,

the bin is labeled hypoxic and available oxygen is consumed. On the other hand, by Fick

diffusion law through a material, oxygen would diffuse from bins with a higher concentration

into adjacent bins with a lower concentration. Under steady state, the oxygen content in

each bin is assumed to be well mixed. Therefore, the whole system is a network of physical

compartments. The resulting equations are similar to the discretized diffusion-consumption

equations but are formulated from completely different principles. In the axial direction of

the capillary, the key strategy in the approach to convection is to use a sliding fluid element

algorithm, with the time step set equal to the length step divided by the fluid velocity. By

matching the time and space steps, the capillary oxygen content flow down the capillary

vessel one segment at a time and movement of single or a bunch of capillaries may stop in

the case of capillary blockage. This not so well-known method, which address the convection,

works perfectly with the compartmental method but has never been applied to a discrete

model. A concrete introduction to the construction and application of the sliding element

method on the discrete compartmental method (DCM) can be found in section 3.3.

17

2.2 Multi-capillary Oxygen Transport in Striated

Muscle

2.2.1 The General Problem

Consider a 3D rectangular parallelopiped of tissue that has width W(cid:48), height H(cid:48), and length
L(cid:48), where the length is parallel to the z-axis (direction of blood flow), the width is parallel

to the x-axis and the height is parallel to the y-axis. The origin of the Cartesian coordinates

is at the center of the rectangular cross-section on the arterial (infusion) end. Consider n
parallel capillaries of same length L(cid:48) and radius R(cid:48)
Blood enters the capillaries at the arterial end with velocity v(cid:48)

c located randomly along the center line.

b. As oxygen is convected

through the capillary, it diffuses into the surrounding tissue both radially and axially under
a constant diffusivity D(cid:48) and is consumed in the tissue at a constant rate κ(cid:48) per volume of
tissue. The oxygen concentration in either blood or tissue is C(cid:48) which is non-negative due

to the physical fact that substrate concentration cannot be negative.
√

Normalize all length by the square root of the cross section area

by that on the arterial end C(cid:48)

W(cid:48) · H(cid:48), concentration
a, time by W(cid:48)H(cid:48)/D(cid:48) and drop all primes. In terms of normalized

parameters, the problem can be stated as: there are n parallel capillaries of same length

L and radius Rc located heterogeneously along the center line of a tissue box that has
√
√
W(cid:48)H(cid:48), W = W(cid:48)/
width W , height H, and length L. Where Rc = R(cid:48)
W(cid:48)H(cid:48), H =
c/
√
H(cid:48)/
W(cid:48)H(cid:48). Blood enters the capillaries at the arterial end with
W(cid:48)H(cid:48)/D(cid:48). While oxygen is convecting through the capillary, it diffuses
velocity vb = v(cid:48)

√
W(cid:48)H(cid:48), and L = L(cid:48)/

√

b

into the surrounding tissue both radially and axially with unit constant diffusivity and is
consumed there at a constant rate κ = κ(cid:48)W(cid:48)H(cid:48)/CaD(cid:48). The oxygen concentration in either

18

blood or tissue is C = C(cid:48)/Ca which is always non-negative. We assume that oxygen has

unlimited permeability through the thin capillary membrane.

2.2.2 Assumptions

1. The system is linear. Blood flow is constant. Both transfer and consumption

coefficients are not concentration-dependent. Capillaries are long and parallel, such

as those in the skeletal muscle or heart muscle. Moreover, the diffusivity of oxygen is

isotropic(axially uniform) in tissue and the permeability of the capillary endothelium

to the oxygen molecules is unlimited.

2. Only convection occurs in the capillary, i.e. oxygen is radially well-mixed in the

capillary. Radial diffusion is infinitely fast inside the capillary so the radial concentra-

tion gradients are zero in the capillary. The radial diffusion relocation times for small

solutes in capillaries is very short compare to capillary transit times and so negligible

due to order of magnitude.

3. The outer domain boundary is reflecting. We considered a local tissue region. We

assume similar diffusion also occurs inside of a similar neighboring domain. We assume

the adjacent tissue regions have exactly the same reflected distribution of capillaries and

identical (reflected) tissue oxygen concentration. This results in zero oxygen exchange

between neighboring regions and no net oxygen flux across the boundary. It allow us

to patch many solutions together and simulate much larger tissue regions.

19

2.3 Discrete Compartmental Systems

A discrete compartmental system is a system which is represented by finite number of in-

terconnected microscopic subsystems of specific size in order, called discrete compartments,

each of which is homogeneous and well-mixed. Each compartment interact with adjacent

neighbors by exchanging material(mass), and for all such transfers between compartments

and, to and from the environment on the boundary, the mass conservation condition holds.

The system equations are difference equations, and the system can be used to approximate

continuous systems described by partial differential equations. There is no specific compart-

mental mental system which is optimal for a continuous system, as it is in the traditional

compartmental methods, but the approximation converges as the number of compartments

is increased as the system is refined.

The discrete compartmental method (DCM) sounds similar to the finite volume method

(FVM) [40] [65] as they both construct algebraic equations based on the flux change over

the small volume surrounding each node point on a mesh. However, they are essentially two

distinct methods and has each of its own advantages. FVM is a method for representing and

evaluating partial differential equations in the form of algebraic equations, but DCM is a

method for analyzing and simulating biological processes in the form of algebraic equations

derived directly from physical laws. In FVM, fluxes at the surfaces of each finite volume are

converted from volume integrals in a partial differential equation that contain a divergence

term, using the divergence theorem, but in DCM, the fluxes are derived from the Fick’s

diffusion law. This intrinsic difference makes DCM a more flexible method and can be applied

to simulate more complicated biological processes. Due to the difference on the construction

method, another remarkable difference between the two methods are for FVM each non-

20

overlapping volume can be of arbitrary size and shape, but for DCM each compartment

mush have the same size and shape.

2.3.1 Basic Terms and Definitions

A number of terms are used frequently in discrete compartmental analysis and in its appli-

cation to the oxygen exchange problem in capillary beds. Variation of these terms in usage

in different literature may cause confusion and overlap in the meaning of some of the terms.

Here we define the commonly used terms more rigorously and will adhere to these definitions

in the remainder of this thesis.

(a) Discrete Compartment and Compartmentalization

A discrete compartment is a small physiological space or physical volume which con-

tains an amount of a material that acts kinetically like a distinct, homogeneous, well-

mixed material. It is to be distinguished from a classical compartment which does not

have a specific size and connectivity [33]. For example, the amount of oxygen in the

whole plasma space of the body could be a compartment in the classical compartmen-

tal system, while the amount of oxygen in a specific small volume of the capillary or

tissue space is a continuous compartment in the discrete compartment system.

Subdividing the spatial domain into homogeneous intervals in one spatial dimension,

homogeneous areas in two spatial dimension, and homogeneous volumes in three spatial

dimension, and represent each of them using discrete compartments is called compart-

mentalization. The shape of discrete compartments can only be triangular, rectangular

or hexagon. Each compartment only interact with adjacent compartments by exchang-

ing mass under a given transfer(diffusion) coefficient, and each compartment may or

21

may not consume a certain amount of mass under a given consumption coefficient. In

each time step we update the values in each compartment by balancing the mass.

(b) Cube Compartment

A cube compartment is a discrete compartment of cubic shape, which means it is three-

dimensional, bounded buy six square faces and has 12 edges and 8 vertices. Size of a

cube compartment is determined by length of each sub-interval in one space dimension,

area of each small square in two space dimension and it is exactly the volume of the

cubic sub-space in three dimensional domain.

Cube compartments are utilized here for the following reasons: 1. it could seamlessly

cover any problem domain(tissue region) without overlapping each other; 2. it is easy

to construct and set up the difference equations of mass transfer between adjacent

compartments under the Cartesian coordinate system; 3. it is computational efficient,

especially in the three dimensional transient problems.

(c) Discrete Compartmental Space

A discrete compartmental space is a three dimensional space(preferably the smallest)

formed by discrete compartments that covers the entire problem domain. For example,

a 4 by 4 cube space formed by 16 small unit cube compartments is a continuous

compartmental space for a sphere of radius 2.

It is worth to notice that in any of

the three space dimensions that make real sense, continuous compartments are in

the form of three-dimensional compartments with volume. In general, a 1D continuous

compartmental space looks like a cylinder, a 2D continuous compartmental space looks

like a slice, and a 3D continuous compartmental space looks like a block.

22

(d) Active Compartment and Inactive Compartment

There are two general types of compartments in a continuous compartmental space:

active compartment and inactive compartment. An active compartment is a com-

partment whose center is inside the problem domain. An inactive compartment is a

compartment whose center is outside the problem domain. Notice that, the concept

“center” is referring to the midpoint of a line segment in 1D space and the center of

a square area in 2D space and the centroid is a 3D space. Inactive compartments are

“ghost” compartments, they do not participate in any computation. The reason for

accepting inactive compartments is for the convenience of programming, since in the

reality computer data are stored in matrices but the problem domain may be irregular.

(e) Adjacent Neighbor and Contiguous Neighbor

Compartment B is called an adjacent neighbor of compartment A, if B has a common

surface with A. Compartment C is called an contiguous neighbor of compartment A,

if C has a common edge or common vertex with A. Each compartment has at most 2

adjacent neighbors, 0 contiguous neighbors in 1D space, at most 4 adjacent neighbors,

4 contiguous neighbors in 2D space, and at most 6 adjacent neighbors, 20 contiguous

neighbors in 3D space.

(Although a compartment only interacts with its adjacent

neighbors, contiguous neighbors may appear in the equation to satisfy some mass

conservation condition on the boundary.)

(f) Interior Compartment and Boundary Compartment

Active compartments can be further classified into two types:

interior compartment

and boundary compartment. An interior compartment is an active compartment whose

adjacent neighbors are all active compartments. A boundary compartment is an active

23

compartment whose has one or more adjacent neighbors that are inactive compart-

ments.

(g) Discrete Compartmental Domain and Domain Boundary

The set of all active compartments is called the discrete compartmental domain; it is

a subset of the discrete compartmental space. The set of all boundary compartments

is called the discrete compartmental domain boundary; it is a subset of the continuous

compartmental domain.

(h) Tissue Compartment and Capillary Compartment

Interior compartments can be eventually divide into two types: source(capillary) com-

partments and consuming(tissue) compartments, based on the position of their center

in the problem domain. For simplicity, we avoid putting a source on the boundary or

even too close to the boundary. Therefore, every boundary compartment is a consum-

ing compartment.

(i) Type of Boundary Compartments

Boundary compartment can be eventually divide into different types based on the

attendance of there active neighbors. Specific methods for classification and list of

boundary types will be discussed for each individual problem.

(j) The Moving Boundary

During a transient process, oxygen deficient regions (transient zero oxygen concen-

tration regions) may appear and their location is unknown apriori. In such a case a

moving boundary exists. Let Oi denote the ith anoxia tissue region which appeared in

24

the domain and “s” the total number of appeared oxygen-deficient regions, then

s(cid:91)

O =

Oi,

(2.1)

is the union of oxygen-deficient regions. Thus, ∂O is the collection of transient moving

i=1

boundaries at the nth time step.

(k) The Sliding Period

A sliding period is the duration of time for elements inside the source compartments

to renew once or move one grid forward once. Its size is determined by the space step

size and the velocity of element convection of the source capillary.

(l) Notation of Operators

For simplicity, we define the following operators in one, two and three dimensional

space respectively:

Σnb

: Σnbun

i = un

i−1 + un

i+1

Σnbun

i,j = un

i−1,j + un

i+1,j + un

i,j−1 + un

i,j+1

Σnbun

i,j,k = un

i−1,j,k + un

i+1,j,k + un

i,j−1,k + un

i,j+1,k + un

i,j,k−1 + un

i,j,k+1

(2.2)

2.3.2 Graphic Representation

A map is a schematic plan of the continuous compartmental space which stores a label for

each compartment indicating there type and any other additional information. Each label

corresponds to a specific difference equation which should be applied to the corresponding

compartment during the computation. The first part of a compartment label is two capital

25

letters indicating the general type of the compartment. For example, for the problem of

oxygen exchange in capillary beds, we have the following leading letters each corresponding

to a basic compartment type:

‘ CC ’ - Capillary Compartment

‘ TC ’ - Tissue Compartment

‘ BC ’ - Boundary Compartment

‘ IC ’ - Inactive Compartment

Another part of the label are numbers distinguishing a fine sorted type or other informa-

tion that indicating a compartment to be treated specifically during the computation. For

example, an interior tissue compartment can be simply marked as ”TC” and an inactive

compartment can be simply marked as ”IC” since there is no need to further distinguish-

ing them. However, two capillary compartments with different strength may be marked as

”CC001” and ”CC002” to specify the coefficients being used during the calculation. Like-

wise, boundary compartments can be marked as ”BC032” or ”BC411” where the number

indicates different boundary types.

The Map itself is stored in an individual matrix which has the same dimension as the

continuous compartmental space it serves. In each time step we check the labels in Map

one by one, and use the paired difference equation to update the value of the corresponding

compartment. This process seem to be repetitive and redundant but actually has great

advantage on the transient multi-capillary and tissue exchange problem. Consider the case

where the blood flows inside one or more capillaries are suddenly decreased by different

amounts due to a small artery occlusion, or the case where a block of tissue’s consumption

26

rate is suddenly increased due to exercise. In both of these cases we will need to make several

adjustments to the coefficients for a specific part of the domain and keep the rest unchanged.

A model which could easily process diverse time-dependent conditions is desired. The label

of a compartment on the Map can be changed at any time step to satisfy time-dependent

conditions. For example, if a capillary is closed (e.g. capillary blockage), we just need to

rewrite its label on the Map to a proper label which corresponds to zero blood velocity. If

the same capillary become open again later, the label is reestablished. If there is a change in

the tissue consumption over part of the domain due to cell hypertrophy, a new tissue label

is created, say ”TC123”, corresponding to an appropriate difference equation.

2.4 Application of the DCM

The computation of discrete compartmental method (DCM) is based on difference equations

derived directly from physical laws. Detailed formulation will be discussed at the beginning

of each individual problems in the following chapters. Here we use a 1D transient problem

to introduce the computational method of discrete compartmental modeling.

Example: 1D problem modeling

Consider a one dimensional tissue region of length L = 1 where a capillary of radius

Rc( Rc < L) is located at the center of the x-y plane in the tissue region. The normalized

diffusion coefficient D is 1, the consumption coefficient is κ, the oxygen concentration in the

capillary is Ca and the blood velocity is vb.

To apply the DCM, first determine the space step based on the size of the capillary.

The position of capillary boundary is of crucial importance to the accuracy of the method.

Thus the spacial step size is determined according to the cross-section area of the capillary.

27

In a one dimensional problem, choice of ∆x should be either equals to the diameter of the

capillary, or there should exist a positive integer whose product with the space step equals

the capillary diameter. Let X be the set of valid choices of ∆x, then

X = { ∆x | ∃ a ∈ Z+ s.t. a∆x = 2Rc }

i.e.

∆x = 2Rc, Rc,

2Rc
3

,

··· .

(2.3)

(2.4)

Determination of ∆x for two and three dimensional problems will be discussed later at the

beginning of each individual problems.

For the time step size, if convection inside the capillary is not considered, then ∆t could

be any small number provided the method converges and stable; if the convection is taken

into account, then a sliding element method is used to deal with the transient convection-

diffusion-consumption process and ∆t is determined by the sliding period. The fundamental

requirement of the sliding element method is to match the space and time step so that the

diffusion-consumption process takes exactly the same amount of time as the materials in

the capillary compartments slide one complete grid forward. Let p be the sliding period

for materials in the capillary compartment slid exactly one grid downstream into the next

compartment, which means blood flow should travel a distance of ∆x in the velocity of vb

during each p. Hence, the sliding period is determined by

p =

∆x
vb

,

(2.5)

which further determines the time step size ∆t. If the stability condition is not satisfied, one

28

should choose to use the quotient of p with a positive integer, say b. Let T be the set of

valid choices of ∆t, then

T = { ∆t | ∃ b ∈ Z+ s.t.

b∆t = p }

i.e.

∆t = p,

p
2

,

p
3

,

··· .

(2.6)

(2.7)

By doing so, the diffusion-consumption is perfumed b times every time before one slid in the

capillary. We call b the sliding frequency.

The oxygen concentration in the next time step of each compartment is determined by

difference equations derived from Ficks’ diffusion law and balance of mass. For example,

oxygen concentration in the next time step of a tissue compartment in the 1D problem is

given by

Cn+1

i

= max

Cn

i +

(cid:26)

(cid:27)
i ] − ∆t · κ, 0

∆t
(∆x)2 [ΣnbCn

i − 2Cn

.

(2.8)

Notice how the maximum function guarantees that the oxygen concentration will never

became negative, which means over each time step, a tissue compartment will either consume
∆t · κ amount of oxygen or whatever is available. (Details about the formulation will be

given at the beginning of each individual problem.)

Instead of considering a Dirichlet boundary condition on the capillary wall, capillary

blood flow is part of the problem and oxygen concentration inside each capillary is computed

simultaneously with tissue oxygen concentration. This could be a great challenge using other

numerical methods but naturally possible in DCM. The advantages of doing so includes: 1.

it provides an out flow at the end of each time step, 2. it could apply to possible cessation

29

of oxygen supply.

Consider when one of the source capillary is blocked. In such case, oxygen will continue

diffuse out from that capillary due to the difference of the partial pressures. The oxygen

concentration inside the capillary well gradually decrease. As the oxygen pressure gradient

gets smaller and the out flow of the oxygen will slow down. When the pressure difference

becomes zero this capillary no longer supply oxygen to the surrounding tissue. Assume there

is a well oxygenated capillary nearby, then the oxygen will gradually flow to regions with low

partial pressure. With the help of nearby normal capillaries, tissue cells next to the blocked

capillary may still be alive. Such complicated and physiological process can be simulated

easily using the discrete compartmental model.

Oxygen concentration of a capillary compartment at the end of the present time step

(denoted by a subscript “+”) is given by

C

n+
i = Cn

i +

∆t
(∆x)2 [ΣnbCn

i − 2Cn
i ] .

(2.9)

Notice the difference of the above equation with Eq. (2.8) is the absence of a consumption

term, and a maximum function is not needed since its magnitude will be positive. If the

capillary occupies more than one compartment, i.e. ∆x < Rc, then an average of the

capillary’s encompassing compartments will be needed to compute the oxygen concentration.

Capillary oxygen concentration will then be updated after each exchange. At the beginning
of the next step (denoted by a subscript “−”)

n+1−
i

C

= Ca.

30

(2.10)

Necessary Components for Stability

In the next few chapters we will gradually construct the governing continuous compartmen-

tal differential equation for capillary-tissue oxygen exchange problem from one dimensional

transient problems to three dimensional transient problems. One essential requirement is

that the resulting method should be convergent and stable. The essential condition that

must be satisfied while applying the continuous compartment method is:

for an d-dimensional transient problem, we require

∆t

(∆x)2

(cid:54) 1
2d

(2.11)

in the compartmental difference equations. Theoretical proof is given in Chapter 5.

2.5 Discussion and Comparison with Previous Models

Two most commonly used methods on studying oxygen transportation problem in human

muscle tissue are the methods based on the Krogh cylinder model and the finite difference

method. They both have several drawbacks on solving the oxygen delivery problem which

can be handled easily by the discrete compartmental method.

A major shortcoming of the Krogh cylinder model and its extensions is the ignorance

of the concurrent flows from neighboring capillaries, which results in unrealistic phenomena

when applying to tissue beds with multiple capillaries. For a stack of Krogh cyliners, blockage

of one capillary will only affect oxygen concentration in that single cylinder, but in reality, as

the concentration of a local area starts decreasing, oxygen will flow in from neighboring tissue

cylinders driven by the partial pressure difference, which results in the cylinder boundary

31

shifting towards the blocked capillary and the cylinder shape will no longer be maintained.

In addition, Krogh cylinder model assumes that capillaries are evenly distributed, which is

not true in reality.

Due to the development of computer technology and advance in computational speed,

numerical methods such as finite difference method has been more and more commonly used

as a tool to solve problems that can be described by partial differential equations. However,

for complicated biological processes, both the pde system and the numerical equations are

difficult to construct. For example, our oxygen transport problem is a three-dimensional

transient convection-diffusion-consumption moving boundary problem. Movement of the

oxygen is governed by different partial differential equations in different parts of the domain.

In addition, location of the moving boundary is unknown. Hence, it is difficult to set up the

boundary condition for pde equations.

The present discrete compartmental modeling is designed for solving transient convection-

diffusion-consumption problems arising from the capillary-tissue oxygen exchange problem.

It has many intrinsic advantages comparing with existing methods. First of all, the resulting

difference equations are derived directly from Fick’s law of diffusion, which makes it easy to

construct and extend to three dimensional problems with practical significance. For example,

it could deal with unevenly distributed capillaries, the occurrence of ischemic regions and

their moving boundaries easily. Second, it is easy to manipulate the time dependent bound-

ary conditions, which are essential for the physiological problem since the consequences of

increasing and decreasing blood flow and tissue consumption rate are of great importance.

Third, it is easy to deal with sources(capillaries) of uneven strengths and locations. More-

over, the difference equation system is similar to that of the finite difference method and

thus easy to compute using modern computer technology.

32

Algorithm Analogy

Imagine a row of sensible kids, each stands one foot away from another. A teacher with

a bag of candy is standing on one end of the row. In the beginning, each kid has 0 candy in

his hand, and the teacher hold 10 candies in her hand. Afterwards, for each minute, each kid

compare the number of candies he has with each one of his two neighbors and exchanges x

candies, where x is the integer difference of number of candies in their hand times 0.4. The

teacher and the kid standing on other end will have only one neighbor. After exchanging

candies, each kid eats 1 candy if he has some or 0 candy if he has nothing. If a kid has

a fraction of one candy he will eat all he has. The teacher does not eat any candy. This

completes one round. At the beginning of each round, the teacher always holds 10 candies

in his hand.

Therefore in the first round, the first kid gets 4 candies ((10 − 0) × 0.4) from the teacher
and 0 candy ((0− 0)× 0.4) from the kid on the right, and eats 1 candy; every other kid gets
0 candy and eat 0 candy. At the end of the first round, the teacher has 6 candies (10− 4) in

his hand, the first kid on the left has 3 candies in his hand, and the rest kids have 0 candies

in their hands. At the beginning of the second round, the teacher holds 10 candies in his
hand again, the first kid on the left has 3 candies (0 + 4 − 1) in his hand, and the rest kids

have 0 candies in their hand. During the second round, the first kid on the left gets 2.8
candies ((10− 3)× 0.4) from the teacher and gives the second kid 1.6 candies ((0− 4)× 0.4);

the second kid on the left gets 1.6 candies from his left neighbor and exchange 0 candy
((0 − 0) × 0.4) with his right neighbor; the rest kids get 0 candy from their neighbor. In the

second round, each of the first two kids on the left eats a candy and each of the rest kids
eats 0 candy. At the end of the second round, the teacher has 7.2 candies (10 − 2.8) in his
hand, the first kid from the left has 3.2 candies (3 + 2.8− 1.6− 1) in his hand, and the second

33

kid from the left has 0.6 candy (1.6 − 1) in his hand, and the rest kids has 0 candies in their

hands. The first round starts at time zero and finishes at the end of the first minute. The

second round starts at the beginning of the second minute and finishes at the end of it.

Notice how balance of mass is hold in the first two rounds. At the end of the second

round, there are totally 3.8 candies (3.2+0.6) left in kids hand. During the first two rounds,

3 candies (1+2) were eaten by the kids. The sum of them is equivalent to the number of

candies that the teacher gives out in the first two rounds, which is 6.8 candies (4+2.8). In

other words, the mass of candies that is being send out by the teacher is conserved with the

mass of candies eaten by the kids plus the mass of candies that left in the kids hands.

There can be more than one teacher, teachers can stand anywhere in the middle of the

line(but they may not stand together), and different teacher can hold different number of

candies in there hand at the beginning of each round. For example, consider a teacher stands

in the middle of the line. The kid on his left has 5 candies, the kid on his right has 9 candies

at the end of the previous step and this teacher hold 12 candies in his hand every time.
Then, in the next step, this teacher will give the kid on his left 2.8 candies((12 − 5) × 0.4)
and the kid on his right 1.2 candies((12 − 9) × 0.4). This teacher has 8 candies(12-2.8-1.2)

left in his hand at the end of this round, but he will take some from the candy bag and has

12 candies again in his hand at the beginning of the next round.

In the candy game, the teacher and kids plays the role of compartments, where the teacher

is a capillary compartment or source compartment that delivers the candy(material) and

each kid is a tissue compartment that consumes the candy(material). The uniform distance

between each two of them, 1 feet, is the space step, and the consistent time duration of each

round, 1 minute, is the time step. The ratio 0.4 multiplied to the candy(mass) difference is

the transfer(diffusion) coefficient, and the number of candies that a kid eats in each round

34

is the consumption coefficient. The number of candies left in each kid’s hand at the end of

each round is the updated (concentration) value of that compartment after each time step.

Notice that, if a kid is standing next to the teacher, then he will always have some candy

to eat and he will have more candies than others have who are further to the teacher; if

a kid is standing very far away from the teacher, then he will gets less candy or perhaps

even no candy at all. Due to the consumption, there is a boundary of the diffusion. For

example, if there is only one teacher standing on the left end and holding only 10 candies

at the beginning of each round, then some kid on the right and every kid on his right may

never get a candy. If a kid do not have any candy to eat over a very long time, he may starve

to death just like tissue cells in human body dying from lack of oxygen.

35

Chapter 3

Unsteady

Diffusion-Consumption-Convection

In this chapter we apply the discrete compartmental method to one and two dimensional

unsteady diffusion-consumption problems. The Crank-Gupta problem[12] is used as a test

case to check the accuracy. The example of multiple capillaries in one dimensional regions

gives an idea of how the model deals with time-dependent boundary conditions. Then we

extend the method to two dimensional space and explore the unsteady state of the capillary

supply region, which attracted much attention in the past but only the steady state solutions

have been obtained[68]. Lastly we consider the problem in a circular domain and discuss the

treatment of boundary condition on irregular boundaries.

3.1 One Dimensional Unsteady Diffusion-Consumption

In general, the problem of one dimensional unsteady diffusion and consumption describing

the oxygen transfer from capillary to tissue after normalization can be stated as follows.
Consider M capillaries with the same radius R(cid:48)
c but with different blood oxygen concentration
C(cid:48)
m, where m = 1, 2, ..., M , located irregularly in a bounded region Ω(cid:48) = [0, W(cid:48)]. The domain

of each capillary is

m =(cid:8)x(cid:48) : x(cid:48) ∈ Ω(cid:48) and |x(cid:48) − x(cid:48)m| (cid:54) R(cid:48)

Ω(cid:48)

C

(cid:9) ,

36

where x(cid:48)

m is the center of each capillary. Therefore, the total capillary domain is

M(cid:91)

Ω(cid:48)
c =

Ω(cid:48)

m

and the absorbing tissue domain is

m=1

Ω(cid:48)
t = Ω(cid:48) \ Ω(cid:48)
c.

Driven by the passive diffusion due to the partial pressure difference, oxygen molecules
permeate out from each capillary, and diffuse into the surrounding tissue at rate D(cid:48) and
were consumed there at a rate of κ(cid:48). The oxygen concentration in the tissue area can be
described as a non-negative function of space and time C(cid:48)(x(cid:48), t(cid:48)), i.e. C(cid:48)(x(cid:48), t(cid:48)) (cid:62) 0 for all x(cid:48)
and t(cid:48). Thus, the governing equation to this one dimensional diffusion-consumption problem

is

∂C(cid:48)
∂t(cid:48) = D(cid:48) ∂2C(cid:48)

∂x(cid:48)2 − κ(cid:48),

x(cid:48) ∈ Ω(cid:48)

t, t(cid:48) ∈ [0,∞) .

with associated zero flux boundary conditions

x(cid:48)(0, t(cid:48)) = C(cid:48)
C(cid:48)

x(cid:48)(W(cid:48), t(cid:48)) = 0.

(3.1)

(3.2)

Inside each capillary and on the capillary wall the concentration is C(cid:48)

m, differing in capillary

strength

C(cid:48)(x(cid:48), t(cid:48)) = C(cid:48)
m,

x(cid:48) ∈ Ω(cid:48)

m, t(cid:48) ∈ [0,∞) , m = 1, ..., M.

(3.3)

Normalize all lengths by W(cid:48), time by W(cid:48)2/D(cid:48) and concentration by the average of that

37

on the arterial end C(cid:48)

aand drop primes. In terms of normalized parameters the equations are

∂C
∂t

=

∂2C
∂x2 − κ,

C(x, t) = Cm,

x ∈ Ωt, t ∈ [0,∞)

x ∈ Ωm, t ∈ [0,∞) , m = 1,··· , M

Cx(0, t) = Cx(1, t) = 0.

(3.4)

(3.5)

(3.6)

where κ = κ(cid:48)W(cid:48)2/C(cid:48)

aD(cid:48).

If oxygen deficient occurs, the flux and the concentration value equals zero on the moving

boundaries ∂On

∂C
∂x

(cid:12)(cid:12)(cid:12)(cid:12)x∈∂On

= C|x∈∂On = 0.

(3.7)

Oxygen deficit region and its moving boundaries are defined in section 2.3.1.

3.1.1 Formulation

Instead of approaching from partial differential equations, the discrete compartmental system

is constructed directly from the physical problem. Consider a one dimensional array of

ordered identical cube compartments with equal longitudinal length ∆x and cross-section

area A. Assume there are m unevenly distributed capillaries. The ∆x is required to be

chosen carefully so that the uniform capillary diameter 2Rc can be divided evenly by ∆x.

For example, if Rc = 0.05 after normalization, then 0.1, 0.05 and 0.025 are good choices for

the ∆x, but 0.09, 0.07 or 0.03 are not.

38

Figure 3.1: Compartmentalization, map and labels of a 1D discrete compartmental model.

Align the center of the capillary-tissue region with the center of a single compartment

with width ∆x and generate the compartment array. Then, each compartment corresponds

to a sub-interval of the capillary-tissue domain of the same length, except the two end

compartments. This discretization error can be reduced by choosing a small and proper ∆x.

The key strategy of constructing the discrete compartmental model is to classify and

label the type of each compartment, depending on the location of their centers. Theses

labels will be stored in a separate array called the map and will be used as a guidance

to assign different algorithms to each compartment during the iteration. We distinguish

between tissue, capillary of different strengths and different type of boundaries. For example,

if the center of a compartment is located inside a capillary region, then this compartment is

marked as a capillary compartment. The first compartment will be marked as a left boundary

compartment (“BC21”) and the last a right boundary compartment (“BC12”). Notice that

because ∆x is appropriately chosen at the beginning, the total volume of each capillary is

kept constant during the discretization. Figure 3.1 shows an example of the construction

of maps of a 1D capillary-tissue region with two capillaries under two different sizes of ∆x.

Different numbers are used as labels for each compartment. In the present example, “TC”

39

means tissue compartment, “CCm” means the mth capillary compartment, “BC21” means

a boundary compartment whose right neighbor is a tissue compartment and “BC12” means

a boundary compartment whose left neighbor is a tissue compartment. As shown in this

figure, the second capillary is shifted a little to the left under a relatively large ∆x, but this

displacement error was reduced greatly under a smaller ∆x.

Figure 3.2: 1D compartmental model

Now we derive the governing difference equations for updating oxygen concentration in

different types of compartments in the 1D discrete compartmental system. Consider the ith
compartment. As shown in Fig 3.2, it is in between the (i− 1)∆x and i∆x. Following Fick’s

law of diffusion, each compartment exchanges oxygen with adjacent compartments. Let Ci

be the oxygen concentration of the ith compartment, then the mass of oxygen molecules in
this compartment can then be expressed as Ci · A · ∆x. Let J1 be the amount of oxygen
flows into the ith compartment from the (i − 1)th compartment through the cross section
area A at (i − 1)∆x and J2 the amount of oxygen leaving the ith compartment towards the
(i + 1)th compartment at i∆x. For simplicity, let Ni denote the mass of oxygen molecules

in the ith compartment for now. Consider small time step ∆t such that

t = ∆t · n.

40

(3.8)

J1J2ii+1i-1ACi·A·ΔxΔxThe net movement of particles across the cross section area A during a time interval ∆t from
the (i − 1)th compartment to the ith compartment is

1 = −N n
J n

i − N n
i−1
A · ∆t

.

(3.9)

Multiplying the top and bottom of the right hand side by ∆x and rewriting, we have

1 = −∆x
J n
∆t

·

(cid:21)

− N n
i−1
A · ∆x

(cid:20) N n

i

A · ∆x

.

(3.10)

The mass concentration is defined as mass per unit volume, and hence

Cn

i =

N n
i
A · ∆x

(3.11)

is the oxygen concentration in the ith compartment. Thus our expression of the net oxygen

flux from the left simplifies to

·(cid:2)Cn

i − Cn
i−1

(cid:3) .

1 = −∆x
J n
∆t

(3.12)

Similarly, the net oxygen flux from the ith compartment to the (i + 1)th compartment can

be expressed as

2 = −∆x
J n
∆t

·(cid:2)Cn

i+1 − Cn

i

(cid:3) .

(3.13)

41

1. Difference Equation of Capillary Compartments:

We start from the case ∆x = 2Rc, which means each capillary is represented by a single

compartment. Assume the ith compartment is a capillary compartment. Since the

permeability for gas through the capillary membrane is unlimited, the oxygen molecules

diffuse from the erythrocyte inside the capillary following Fick’s diffusion law under the

normalized diffusive constant D = (∆x)2/∆t (substitute back the parameters used for
normalization and their units, we have W(cid:48)2/(W(cid:48)2/D(cid:48)) = D(cid:48)). Since the capillary itself

does not consume any oxygen, the net oxygen flux difference in the ith compartment

can be found by taking difference of the inflow 3.12 and outflow 3.13

(cid:2)Cn

i−1 − 2Cn

i + Cn

i+1

(cid:3) .

(3.14)

1 − J n
J n

2 =

∆x
∆t

Balancing the mass flow though the cross section area A with the change of oxygen

concentration in the ith capillary compartment over a short time period ∆t yields

D ·(cid:104)

n−
i − Cn

i

C

(cid:105) · A · ∆x = [J n

1 − J n

2 ] · A · ∆t,

(3.15)

where n− represents the time at the end of the nth time interval. Substitute 3.15 into

3.14 and simplify, the oxygen concentration in a capillary compartment at the end of

the nth time interval is

n−
i = Cn

i +

C

∆t

(∆x)2

(cid:2)Cn

i−1 − 2Cn

i + Cn

i+1

(cid:3) .

(3.16)

The oxygen concentration of the mth capillary compartment in the begining of the

42

next time interval is

Cn+1

i

= Cm, m = 1, ..., M.

(3.17)

Now we extend the above idea to the general case, where ∆x ≤ 2Rc. One major

problem arises from this situation is the asymmetry of neighboring tissue oxygen con-

centrations on the two sides. The strategy is to construct an algorithm from the real

physiological process. Let the capillaries and tissues complete exchanging and consum-

ing oxygen using equation 3.16 during the time interval ∆t, then average the oxygen
inside the capillary at the end of each time interval t− = ∆t · n−.

Notice that, there is no need to distinguish compartments of the same capillary between

the interior and the boundary ones and assign with different difference equations. The

discrete compartmental method can self-adjust during the physiological processes.

Figure 3.3: Calculation of capillary oxygen concentration after blocking the source. The
three darker cells denote contiguous capillary compartments.

43

The MapInitial ConditionExchangeMix t = 0+ t = Δt- t = Δt+ExchangeMix t = 2Δt- t = 2Δt+ TC     TC    CC1  CC1   CC1    TC     TC   0.5     0.5      1        1        1        0       0 0.5    0.6     0.9      1       0.8    0.2       0 0.5    0.6     0.9     0.9     0.9    0.2       0 0.6    0.7     0.8     0.9     0.7    0.3     0.1 0.6    0.7     0.8     0.8     0.8    0.3     0.1Figure 3.3 demonstrates the change in oxygen concentration after blocking the source

in the case of ∆x = 2

3 Rc.

In this example, the three capillary compartments were

labeled as ”CC88” and the surrounding compartments are labeled as ”TC” in the

map. The initial oxygen concentration in each capillary compartment is 1 and the

oxygen concentration is 0.5 in the two left neighboring tissue cells but 0 in the two

right neighboring tissue cells. Due to the uneven partial pressure on the two sides,

the oxygen concentrations inside each of the three capillary compartment are different

after one step of diffusion and the concentration in the interior capillary compartment

is unaffected.

In this situation, the average of the three capillary compartments is

calculated and reassigned to all of three for the beginning of the next diffusion and

consumption step. Notice in reality the mixing process happens concurrently with the

diffusion, therefore the process does not take any extra time during the computations.

2. Difference Equation of Tissue Compartments:

Oxygen is consumed in the tissue by chemical reaction (to CO2) at a constant rate

κ. Thus, combining Fick’s first law with the mass conservation, the change of oxygen

concentration in the ith compartment over a short time period ∆t is

D ·(cid:104)

Cn+1

i

− Cn

i

(cid:105) · A · ∆x = [J n

1 A − J n

2 A − κ · A · ∆x] · ∆t.

(3.18)

Substitute in J n

1 , J n

2 from 3.12 and 3.13, oxygen concentration in an tissue compart-

ment at the next time step can be estimated by

Cn+1

i

= Cn

i +

∆t

(∆x)2

(cid:2)Cn

i−1 − 2Cn

i + Cn

i+1

(cid:3) − ∆t · κ.

(3.19)

44

However, in the case of hypoxia there may not be sufficient oxygen for tissue consump-

tion. In other words, the local available oxygen is less than the necessary consumption

rate

Cn

i +

∆t

(∆x)2

(cid:2)Cn

i−1 − 2Cn

i + Cn

i+1

(cid:3) < ∆t · κ,

(3.20)

which will result the oxygen concentration value to be negative in the next time step.

In this case, we require that each tissue compartment can only consume the amount

of oxygen that are available, which means such tissue compartment will consume all

oxygen it has and its the concentration in the next time step is set to zero

Cn+1

i

= 0.

(3.21)

Combining the two cases, the governing equation of an interior tissue compartment

can be written as

Cn+1

i

= max

(cid:26)

Cn

i +

∆t

(∆x)2

(cid:2)Cn

i−1 − 2Cn

i + Cn

i+1

(cid:27)
(cid:3) − ∆t · κ, 0

.

(3.22)

Notice that, in the above equation ∆t/(∆x)2 must be less than 0.5 for stability and

convergence in one dimensional problems. (A proof is given in Chapter 5.)

3. Difference Equation of Boundary Compartments:

There are two type boundary compartments in the one dimensional discrete compart-

mental model, the left boundary compartment and the right compartment. Table 3.1

shows a schematic diagram of the two boundary compartments and their labels. Here,

label “BC21” represents the left boundary compartment and label “BC12” represents

45

the right boundary compartment.

IC BC21 TC

TC ... TC

CC1

TC ... TC BC12

IC

Table 3.1: 1D schematic diagram of different type of boundaries and their labels. Label
“BC21” = the left boundary compartment and label “BC12” = the right boundary com-
partment. “CC1” = capillary compartment of the (m = 1) capillary. “IC” = inactivated
compartment.

Consider the left boundary compartment, labeled as “BC21”, which has oxygen con-

centration Cn

1 . Based on equation 3.22, its oxygen concentration in the next time step

is

Cn+1

1

= max

(cid:26)

Cn

1 +

∆t
(∆x)2 [Cn

0 − 2Cn

1 + Cn

,

(3.23)

(cid:27)
2 ] − ∆t · κ, 0

where Cn

0 is an inactive compartment on the left of the left boundary compartment

and is located outside the domain boundary. Recall that the outer domain boundary

is assumed to be reflecting, which implies that the oxygen concentration of the 0th

compartment is always identical to the concentration of the 2nd compartment and so

there is no net flux cross the boundary, i.e. Cn

0 ≡ Cn

2 for all n. Therefore, the difference

equation gives the oxygen concentration of the left boundary compartment in the next

time step can be generally written as

(cid:26)

(cid:2)Cn

i+1 − Cn

i

(cid:27)
(cid:3) − ∆t · κ, 0

.

(3.24)

Cn+1

i

= max

Cn

i +

2∆t
(∆x)2

Similarly, for the right boundary compartment, labeled as “BC12”, oxygen concentra-

tion in the next time step is

Cn+1

i

= max

(cid:26)

Cn

i +

(cid:2)Cn

i−1 − Cn

i

(cid:27)
(cid:3) − ∆t · κ, 0

.

(3.25)

2∆t
(∆x)2

46

Summary

There are in total four types of compartments in the one dimensional problem. The difference

equation governing each of them can be summarized to one single equation:

(cid:26)

Cn+1

i

= max

Cn

i +

(cid:2)a · Cn

i−1 − 2Cn

i + b · Cn

i+1

∆t

(∆x)2

(cid:27)
(cid:3) − ∆t · K, 0

.

(3.26)

Different parameters should be employed while applying to different type of compartments.

Table 3.2 gives a summary of the corresponding relationships of the compartment label and

the coefficients.

Compartment Type

a b K

CCm

TC

BC21

BC12

1

1

0

2

1

1

2

0

0

κ

κ

κ

Table 3.2: Index table of different types of compartments. The corresponding relationships
of the compartment label and the coefficients is given. Label “CCm” represents the mth cap-
illary compartments, label “TC” represents the interior tissue compartments, label “BC21”
represents the left boundary compartment and label “BC12” represents the right boundary
compartment. K indicates the value of the consumption term.

3.1.2 The Crank-Gupta Problem

Introduction

The one dimensional moving boundary problem describing the oxygen diffusion and con-

sumption phenomenon in the body tissue is first discussed by Crank and Gupta [12] [11]

47

due to interest in medical research concerning the uptake of oxygen by tissue. As one of

the simplest examples of a moving boundary problem, it has been used as a test case in the

development of new methods thereafter and is known as the Crank-Gupta problem. Hansen

and Hougaard [29], obtained an integral formula for the concentration. Their results are

probability the most accurate available and we shall compare our results with both Crank

and Gupta [12] and Hansen and Hougaard [29].

Statement of the problem

In the Crank-Gupta problem, oxygen is allowed to diffuse into a medium which simulta-

neously consumes the oxygen at a constant rate. The problem has two parts: first, the

oxygen concentration at the surface of the medium(could be considered as a capillary) is

kept constant until oxygen does not diffuse any further into the medium and a steady-state

is reached. The supply of oxygen is then removed abruptly and the surface is sealed. No

further oxygen passed in or out and the medium continues to absorb the available oxygen.

Consequently, the point marking the furthest depth of diffusion in the steady-state begins

to move towards the sealed surface at x = 0.

In non-dimensional form, the second part of the problem can be expressed in discrete

compartmental system as an array of I + 1 cube compartments of side length ∆x = 1/I,
each of which is centered at position ∆x · i and has oxygen concentration Cn

i at time ∆t · n.

Initially, each compartment has oxygen concentration

C0

i =

1
2

(1 − ∆x · i)2,

i = 0, 1,··· , I.

(3.27)

At the beginning, each compartment exchange oxygen with its neighbors over each small

48

time period ∆t with constant diffusion coefficient of D = 1 and consume the oxygen simul-

taneously at a constant rate of κ = 1.

Using partial differential equations, the problem is originally expressed as

∂C
∂t

=

∂2C
∂x2 − 1,

0 (cid:54) x (cid:54) s(t),

t (cid:62) 0

(3.28)

with associated boundary conditions

∂C
∂x

C =

∂C
∂x

and the initial condition

= 0,

at the sealed surface x = 0,

(3.29)

= 0,

at the moving boundary x = s(t),

(3.30)

C =

1
2

(1 − x)2,

0 ≤ x ≤ 1,

t = 0,

(3.31)

with s(0) = 1. Here C(x, t) denotes the oxygen concentration at a point, x, and time, t, and

s(t) the position of the moving boundary.

Notice that, the major challenge of the problem is to trace the location of the moving

boundary, in the meanwhile determining the oxygen concentration through the medium as a

function of time. Since not only the oxygen concentration is zero but its gradient is also zero

on the boundary, there is no explicit relationship which contains the velocity of the moving

boundary. However, if the discrete compartmental model is used to analyze the problem,

there is no need to specify the location of the boundary.

49

Implementation of the DCM

To apply the discrete compartmental method, one need to construct an appropriate com-

partmental map and to determine the label of each compartment. Table 3.3 shows the

compartmental map for the Crank-Gupta problem with the following properties: 1. only

one boundary compartment is needed on the left side, 2. the number of tissue compartments

depends on the choice of ∆x, 3. there is no need to specify the right boundary compartment,

4. the nearest compartment with zero concentration would be the moving boundary.

IC BC21 TC TC TC TC TC TC TC TC TC ...

Table 3.3: Compartmental map for the Crank-Gupta problem.

The difference equation updating the value of each compartment corresponding to each

label is given by equation (3.26) and coefficients of which is given by table 3.2.

Results and discussion

Following the discrete compartment model, computations were carried out for different com-

partment sizes (∆x = 0.05 and ∆x = 0.0125).

Interpolation is applied to improve the

accuracy of the results. Table 3.4 and table 3.5 shows the comparison of results to that

from Crank and Gupta [12] with ∆x = 0.05 and Hansen & Hougaard [29] at different times.

Before a moving boundary appears, i.e. Eq. (3.30), our result is very close to that of Crank

and Gupta [12] using a finite difference method. At all times, the values obtained using

∆x = 0.0125 are in a very good agreement with those from Hansen & Hougaard [29]. It

should be noted that both numerical solutions using relatively larger grid (∆x = 0.5) involve

large errors on the moving boundary and are unstable due to discontinuity in the gradient

there. As the moving boundary gets very close at larger t (e.g.

t = 1.95), that error is

50

enlarged since too few grids are left and are not enough to represent the remaining subtle

changes.

t\x

0.060

0.100

0.120

0.140

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

223605

217168

198777

170979

137379

101928

68291

39451

17581

4129

223746

217330

198992

171251

137684

102227

68548

39645

17705

4186

223615

217175

198783

170980

137376

101920

68283

39443

17574

4118

223605

217268

198777

170979

137379

101928

68291

39451

17579

4122

143177

139294

128082

110787

89295

65892

43018

23059

8232

603

143287

139414

128338

110996

89502

66112

43228

23232

8342

619

143181

139298

128084

110787

89289

65878

42995

23041

8271

774

143177

139294

128082

110787

89295

65892

43018

23058

8231

608

109129

106019

97025

83117

65795

46941

28667

13204

2873

109228

106125

97149

83265

65963

47115

28827

13324

2924

109131

106020

97024

83114

65790

46936

28661

13193

2838

109129

106019

97025

83117

65795

46941

28667

13203

2875

0

0

0

77850

75351

68130

56988

43197

28416

14638

4204

77937

75442

68233

57105

43322

28536

14730

4249

0

0

77850

75350

68128

56984

43190

28407

14640

4186

77850

75351

68130

56988

43197

28415

14637

4203

Table 3.4: Comparison for C × 106 at different times. On each time level, the first row shows
values taken from Hansen & Hougaard [29], the second row shows values from Crank and
Gupta [12] with ∆x = 0.05, the third row gives values obtained from DCM with ∆x = 0.05
and the forth row the DCM with ∆x = 0.0125.

51

t\x

0.150

0.160

0.180

0.190

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

63078

60845

54403

44503

32353

19583

8251

1005

63157

60928

54496

44602

32453

19668

8298

1007

0

0

63069

60837

54394

44494

32344

19570

8217

996

63077

60845

54403

44503

32353

19582

8250

1000

48823

46840

41136

32434

21927

11304

2890

48893

46912

41212

32511

21996

11346

2890

48822

46837

41132

32428

21918

11303

2874

0

0

0

48823

46840

41136

32433

21927

11303

2888

21781

20287

16066

9942

3523

21824

20328

16096

9950

3506

21776

20282

16058

9932

3506

21781

20287

16065

9941

3519

0

0

0

0

0

9021

7817

4578

799

9039

7827

4575

750

9007

7801

4554

788

9021

7817

4576

787

0

0

0

0

0

0

2884

1914

32

0.195

2880

1909

2867

1885

0

0

2883

1912

33

0

0

0

0

0

0

0

Table 3.5: Comparison for C × 106 at different times. On each time level, the first row shows
values taken from Hansen & Hougaard [29], the second row shows values from Crank and
Gupta [12] with ∆x = 0.05, the third row gives values obtained from DCM with ∆x = 0.05
and the forth row the DCM with ∆x = 0.0125.

52

Method\Time

0.05

0.1

0.16

0.18

0.19

0.195

0.197

0.1974

Crank & Gupta [12]

247841

143287

48893

21834

9030

2880

Hansen & Hougaard [29]

247687

143177

48893

21824

9039

2884

-

-

Dahmardah & Mayers [13]

247687

143177

48823

21781

9021

2883

491

DCM(∆x=0.0125)

247687

143177

48823

21781

9021

2883

503

-

-

40

39

Table 3.6: Comparison of results at the sealed surface.

Table 3.6 provides a comparison for the result at the sealed surface as obtained from (i)

Crank & Gupta [12] with ∆x = 0.5, (ii) Hansen & Hougaard [29] using integral method,

(iii) Dahmardah & Mayers [13] using Fourier expansion and (iv) DCM with ∆x = 0.0125.

Our results are in good agreement with those of Hansen & Hougaard [29] and Dahmardah

& Mayers [13] even in large times.

Method

Oxygen depletion time

Crank & Gupta [12]

0.196349

Hansen & Hougaard [29]

0.1972-0.1977

Miller, Morton & Baines [45]

0.196-0.198

Gupta & Kumar [28]

(unequal intervals)

0.19732

0.19734

Dahmardah & Mayers [13]

0.197434

Mitchell [46]

0.197432

Mitchell & Vynnycky [47]

0.197435

DCM (∆x=0.0125)

0.197437

Table 3.7: Comparison of the oxygen depletion time with previous works.

53

Table 3.7 shows the comparison of oxygen depletion time with many previous works and

Table 3.8 shows the comparison of oxygen depletion time in different mesh size using DCM

. It is seen from the tables that our result converges quickly as the size of ∆x decreases and

agree with those accurate methods with ∆x = 0.0125 and smaller.

∆x 0.1

0.05

0.025

0.0125

0.00625

0.003125

t*

0.196

0.197

0.1975

0.197437

0.197437

0.197437

Table 3.8: Comparison of oxygen depletion time (t*) in different mesh size using DCM.

3.1.3 Three Capillaries of Uneven Strength

As a multi-capillary example, consider three capillary of diameter 0.02 and uneven strength in

a one-dimensional tissue region of unit length and unit consumption rate. As an illustrative

example, let the first capillary centered at x = 0.21, the second at x = 0.61 and the third at

x = 0.91. Let the consumption rate be κ = 1. Computations in this example were carried

out with compartment size of ∆x = 0.02 and time step size ∆t = 0.0001. The analytical

solution used for comparison is found by solving differential equation (3.4) over each interval.

(a) The oxygen concentration is originally zero everywhere and the capillary concentration

is suddenly raised and kept at 0.08, 0.03 and 0.05 respectively. There is no flux on

the tissue boundary. All values here are normalized. Figure 3.4 shows the transient

concentration profiles at t = 0.001, t = 0.01 and t = 0.2. An analytic steady state
solution is obtained from solving the partial differential equations. As t → ∞, the

54

transient numerical solution approaches the analytic steady state solution



0.5(cid:0)x2 + 0.12(cid:1)
0.5(cid:0)x2 − 1.083x + 0.3499(cid:1)
0.5(cid:0)x2 − 1.77x + 0.5294(cid:1)
0.5(cid:0)x2 − 2x + 1.094(cid:1)

C(x) =

0 ≤ x ≤ 0.2

0.22 ≤ x ≤ 0.6

0.62 ≤ x ≤ 0.9

0.92 ≤ x ≤ 1.

(3.32)

Figure 3.4: Oxygen diffusion from three capillaries of different strengths into a consuming
region at t=0.001, t=0.01 and t=0.1, compared with the analytic steady state solution.

(b) The initial oxygen concentration is at the steady state of part (a). The second capillary

is blocked at t = 0. Note that a blockage does not mean the capillary concentration

suddenly drops to zero. A blocked capillary no longer provides new oxygen but the

remaining oxygen inside would gradually diffuse out. Figure 3.5 shows the transient
concentration profiles at t = 0.001, t = 0.02 and t = 0.2. As t → ∞, the transient

55

�=������=�����=�����������������������������������������solution approaches

C(x) =



0.5(cid:0)x2 + 0.12(cid:1)
0.5(cid:0)x2 − 1.208x + 0.3774(cid:1)
0.5(cid:0)x2 − 2x + 1.094(cid:1)

0 ≤ x ≤ 0.2

0.22 ≤ x ≤ 0.9

0.92 ≤ x ≤ 1.

(3.33)

which is the analytic steady state solution under the new initial conditions.

Figure 3.5: Concentration profile at t=0.001, t=0.02 and t=0.2 times after a blockage oc-
curred on middle capillary, compared with the analytic steady state solution.

(c) The initial oxygen concentration is at the steady state of part (b). The third capillary

is blocked at t = 0. Figure 3.5 shows the transient concentration profiles at t = 0.002,

t = 0.02 and t = 0.3. As t > 0.0272, the tissue compartments in the middle of the right

two capillaries begin to suffer anoxia and two moving boundaries occur. At t = 0.03, an
anoxic zone of zero oxygen concentration exists in the interval 0.6 ≤ x ≤ 1. Note that

the gradient of the oxygen concentration at the supply boundary is zero as predicted

(Eq. (3.30)).

56

�=������=�����=�����������������������������������������Figure 3.6: Concentration profile at t=0.002, t=0.02 and t=0.03 times after a blockage
occurred on the right capillary.

Table 3.9 shows how the boundary condition is easily changed in DCM by relabeling the

capillary compartments. Each raw of labels corresponds to part (a), (b) and (c) in the above

example. At the beginning, there are three capillaries labeled as “CC1”, ‘1CC2” and “CC3”

respectively, each of which determines the amount of oxygen concentration for that capillary

compartment at the beginning of each time step. While the second capillary is blocked,

label “CC2” is replaced by “CC0”, which corresponds to a particular algorithm for blocked

capillaries (diffusion only and no further update of incoming oxygen concentration). While

the third capillary is blocked, label “CC3” is replaced by “CC0” as well. Besides relabeling,

no further action is requires for changing the boundary condition.

... TC CC1

TC ... TC

CC2

TC ... TC CC3 TC ...

... TC CC1

TC ... TC

CC0

TC ... TC CC3 TC ...

... TC CC1

TC ... TC

CC0

TC ... TC CC0 TC ...

Table 3.9: Schematic diagram of label changing following capillary blockage in the 1D ex-
ample. Label “CCM” = capillary compartment corresponding to the oxygen concentration
of the mth capillary. “CC0” = capillary compartment corresponding to a blocked capillary.

57

�=������=�����=������������������������������������������Discussion

The solution of part (a) shows that with capillaries of different strengths, tissue oxygen

concentration is changing at different speeds and the supply boundaries are closer to the

weaker sources. Part (b)’s result shows that blockage of a capillary may not cause serious

oxygen deficit overtime if neighbor sources are strong enough to compensate. Part (c) shows

the method can easily treat moving boundary problems. Such problems has no analytic

solutions to compare with.

Transient solutions of this one dimensional problem illustrate many innovative advances

of DCM. However, it also shows the limits of using low dimensional models. The zero
concentration profile over the domain 0 ≤ x ≤ 0.2 never recovered. This is because there

is only one pathway in the 1D domain and tissue regions are isolated. If there are other

pathways, would the oxygen could diffuse around any blockage and replenish anoxic tissues?

We shall explore this possibility by building 2D and 3D compartmental models.

58

3.2 Two Dimensional Unsteady Diffusion-Consumption

The discrete compartment model has showed some advantages on solving the one dimensional

unsteady diffusion-consumption problem.

In 2D, it is more difficult to predict when and

where the anoxia regions will appear. The advantages of DCM will be extended to handle

such difficulty.

Figure 3.7 shows a schematic diagram of 2D capillary-tissue region with two capillaries

of uneven strength. An important question we wish to answer is, where and when would

anoxia occur if one or more capillaries was blocked. For example, if the capillary on the right

was blocked, will the ischemia area only happen around the capillary or elsewhere or both?

Figure 3.7: schematic diagram of 2D planar capillary-tissue region and blockage on one
capillary. Shaded area in part (b) represents the possible anoxic areas due to blockage on
the right capillary.

The problem of two dimensional unsteady diffusion and consumption problem describing

59

the oxygen transfer from capillary to tissue after normalization can be generally stated as
follows. Consider M capillaries of same radius R(cid:48)
located irregularly in a bounded region Ω(cid:48) ∈ R2. Define individual capillary domain

c and different oxygen concentration C(cid:48)

M

(cid:110)
(x(cid:48), y(cid:48)) : (x(cid:48), y(cid:48)) ∈ Ω(cid:48) and (x(cid:48) − x(cid:48)

Ω(cid:48)
m =

m)2 + (y(cid:48) − y(cid:48)

m)2 (cid:54) (R(cid:48)

c)2(cid:111)

,

where (x(cid:48)

m, y(cid:48)

m) are the coordinates of the center of each capillary. Thus, the total capillary

domain is

and the absorbing tissue domain is

M(cid:91)

m=1

Ω(cid:48)

m

Ω(cid:48)
c =

Ω(cid:48)
t = Ω(cid:48) \ Ω(cid:48)
c.

Driven by the passive diffusion due to the partial pressure differences, oxygen diffuse out

from each capillary into the surrendering tissue and is consumed. The governing partial

differential equation to this unsteady diffusion-consumption problem is

∂t(cid:48) = D(cid:48)(cid:18) ∂2C(cid:48)

∂C(cid:48)

∂x(cid:48)2 +

∂2C(cid:48)
∂y(cid:48)2

(cid:19)

(cid:0)x(cid:48), y(cid:48)(cid:1) ∈ Ω(cid:48)

t, t(cid:48) ∈ [0,∞)

− κ(cid:48),

(3.34)

where D(cid:48) is the diffusion constant, C(cid:48)(x(cid:48), y(cid:48), t(cid:48)) is the non-negative oxygen concentration
(or partial pressure) and k(cid:48) is the constant consumption per volume of tissue. Inside each
capillary and on the capillary wall the concentration is C(cid:48)

M , differing in the capillary strength

C(cid:48)(x(cid:48), y(cid:48), t(cid:48)) = C(cid:48)
m,

m, t(cid:48) ∈ [0,∞) , m = 1,··· , M.

(3.35)

(cid:0)x(cid:48), y(cid:48)(cid:1) ∈ Ω(cid:48)

60

On the tissue boundary ∂Ω(cid:48) there is no flux or the normal derivative of concentration is zero

(cid:12)(cid:12)(cid:12)(cid:12)(x(cid:48),y(cid:48))∈∂Ω(cid:48) = 0.

∂C(cid:48)
∂n

Normalize all length by ˜L =(cid:112)A(Ω(cid:48)), the square root of Ω(cid:48)’s area, the concentration by the

average on the arterial end C(cid:48)

a, time by ˜L2/D(cid:48) and drop primes. In terms of normalized

(cid:19)

∂2C
∂y2

− κ,

(x, y) ∈ Ωt, t ∈ [0,∞) ,

C(x, y, t) = Cm, ,

(x, y) ∈ Ωm, t ∈ [0,∞) , m = 1,··· , M,

parameters the equations are

∂x2 +

=

(cid:18) ∂2C
(cid:12)(cid:12)(cid:12)(cid:12)(x,y)∈∂Ω

∂C
∂t

∂C
∂n

= 0,

(3.36)

(3.37)

(3.38)

(3.39)

(3.40)

If oxygen deficient occurs, the flux and the concentration value equals zero on the moving

where κ = κ(cid:48) ˜L2/C(cid:48)

aD(cid:48).

boundaries ∂O

∂C
∂n

(cid:12)(cid:12)(cid:12)(cid:12)(x,y)∈∂O

= C|(x,y)∈∂O = 0

(3.41)

Oxygen deficit region and its moving boundaries are defined in section 2.3.1.

3.2.1 Formulation

Now we construct the 2D discrete compartmental system directly from the physical problem.

The general problem after normalization can be stated as follows.

Consider a two dimensional discrete compartmental space of packed ordered identical

61

compartments with equal side length ∆x and cross-section area A. The key of solving

capillary-tissue exchange problem using the discrete compartmental model is to maintain

the cross section area of each capillary. Moreover, the shape of capillary compartments need

be approximately radially symmetric. Thus, in two dimensional modeling, the choice of ∆x

must satisfy

∆x =

(cid:115)

πR2
c
nCC

(3.42)

where nCC = 1, 4, 5, 9, 13, .... So that the total area of the capillary compartments is consis-

tent with the original circular area and the direction of oxygen flux is kept radially uniform.

Figure 3.8 shows some examples of good and bad capillary compartment arrangements.

Figure 3.8: Example of 2D capillary shapes.

Now we derive the governing difference equations for updating oxygen concentration in

different types of compartments in 2D discrete compartmental system. Figure 3.9 shows the

local arrangement of nine cube compartments. The (i, j)th cell has four adjacent neighbors

(i − 1, j), (i, j − 1), (i + 1, j), (i, j + 1),

(3.43)

62

as well as four diagonal neighbors

(i − 1, j − 1), (i − 1, j + 1), (i + 1, j − 1), (i + 1, j + 1).

(3.44)

Let Cn

i,j denote the oxygen concentration of the (i, j)th compartment at time ∆t · n. The

mass change inside this closed volume V = (∆x)3 over a short time period ∆t is

(cid:104)

(cid:105) · V.

i,j − Cn
Cn+1

i,j

(3.45)

Figure 3.9: The 2D compartmental model.

Following the work in section 3.1.1, the oxygen fluxes across the cross sections of an area

A from/to the top, left, bottom or right neighbor compartments during a time interval ∆t

63

i-1,j+1i,j+1J1J2J3J4i,ji-1,j-1i+1,ji-1,j"i+1,j-1i,j-1i+1,j+1AAΔxΔxare

·(cid:104)
·(cid:104)

i,j − Cn
Cn
i−1,j
i+1,j − Cn
Cn

i,j

(cid:105)
(cid:105)

,

,

1 = −∆x
J n
∆t
3 = −∆x
J n
∆t

2 = −∆x
J n
∆t
4 = −∆x
J n
∆t

·(cid:104)
·(cid:104)

i,j − Cn
Cn
i,j−1
i,j+1 − Cn
Cn

i,j

(cid:105)
(cid:105)

,

.

(3.46)

1. Difference Equation of Capillary Compartments:

Consider the case ∆x =(cid:112)πR2

c , which means each capillary is represented by a single

compartment. Assume the (i, j)th compartment is a unobstructed capillary compart-

ment, where oxygen diffused out over each time step following the normalized diffusion

constant D = (∆x)2/∆t and the compartment itself does not consume any oxygen.

Thus, the concentration change in the (i, j)th capillary compartment over the time

period ∆t can be described as

[(J n

1 + J n

2 − J n

3 − J n

4 ) · A] · ∆t.

(3.47)

Balance of mass flow through the cross section area A with the oxygen concentration

change in the i, jth capillary compartment over the time period ∆t yields

D ·(cid:104)

n−
i,j − Cn

i.j

C

(cid:34)(cid:32) 4(cid:88)

(cid:105) · V =

(cid:33)

(cid:35)

J n
m

· A

· ∆t,

(3.48)

where n− represents time at the end of the nth time interval. Substituting in J n
1 ,

m=1

J n
2 , J n

3 , J n

4 , from (3.46) and simplifying give the oxygen concentration in the (i, j)th

capillary compartment at the end of the nth time interval

(cid:104)

ΣnbCn

i,j − 4Cn

i,j

(cid:105)

,

(3.49)

n−
i,j = Cn

i,j +

C

∆t

(∆x)2

64

where

ΣnbCn

i,j = Cn

i−1,j + Cn

i,j−1 + Cn

i+1,j + Cn

i,j+1.

(3.50)

The oxygen concentration of the mth unobstructed capillary compartment in the be-

gining of the next time interval is

Cn+1

i,j = Cm.

(3.51)

Now extend the above idea to the general case, where ∆x ≤(cid:112)πR2

c . As shown in figure

3.8(a), while one capillary is represented by more then one compartment, each of the

compartment may have different values after one exchange process and a compartment

in the middle is isolated to the tissue region whose value will not be changed during

the exchange process. Since oxygen is assumed to be well mixed inside the capillary,

an average of the oxygen concentration over each capillary’s compartments shall be
taken and reassigned back to them at the end of each time interval t− = ∆t · n−.

For capillaries with reducing strength, a new value of Cm should be used in the next
time step tn+1 = ∆t · (n + 1). For blocked capillaries,

Cn+1

i,j = C

n−
i.j .

(3.52)

2. Difference Equation of Tissue Compartments:

Consider the interior tissue compartments which consume oxygen at a constant rate

of κ. The concentration change caused by the difference between inflows, outflows and

65

consumption in the i, jth compartment over the time period ∆t can be described as

[(J n

1 + J n

2 − J n

3 − J n

4 ) · A − κ · V ] · ∆t.

Now, balance Equation (3.53) with Equation (3.45) gives

D ·(cid:104)

i,j − Cn
Cn+1

i.j

(cid:34)(cid:32) 4(cid:88)

(cid:105) · V =

m=1

(cid:33)

J n
m

· A − κ · V

(cid:35)

· ∆t.

Substituting in J n

1 , J n

2 , J n

3 , J n

4 , from (3.46) and simplifying yield

(cid:104)

ΣnbCn

i,j − 4Cn

i,j

(cid:105) − ∆t · κ,

Cn+1
i,j = Cn

i.j +

∆t

(∆x)2

where

ΣnbCn

i,j = Cn

i−1,j + Cn

i,j−1 + Cn

i+1,j + Cn

i,j+1.

Consider local oxygen deficiency, where

(cid:104)

(cid:105)

< ∆t · κ.

Cn

i,j +

∆t

(∆x)2

ΣnbCn

i,j − 4Cn

i,j

(3.53)

(3.54)

(3.55)

(3.56)

(3.57)

In this case, the tissue compartment consumes all the remaining oxygen during that

time step, which results

Cn+1

i.j = 0.

(3.58)

Combining the oxygen sufficient case with the oxygen deficient case, the governing

66

difference equation of an interior tissue compartment is

(cid:26)

(cid:104)

Cn+1

i,j = max

Cn

i,j +

∆t

(∆x)2

ΣnbCn

i,j − 4Cn

i,j

(cid:27)
(cid:105) − ∆t · κ, 0

,

(3.59)

where ∆t/(∆x)2 shall be less than 0.25 in the two dimensional problem for stability

and convergence. See section 5 for a theoretical proof.

3. Difference Equation of Boundary Compartments:

For a rectangular shaped domain in two dimensional space, there are eight type of

boundary compartments. As shown in table 3.10, these boundary compartments can

be classified into two groups: tissue compartment with one inactive neighbor(located on

the side) and tissue compartment with two inactive neighbors(located on the corner).

IC

IC

IC

IC

IC

IC BC40 BC31 BC50

IC

IC BC21

TC BC22

IC

IC BC60 BC32 BC70

IC

IC

IC

IC

IC

IC

Table 3.10: Discrete compartmental map of 2D unsteady diffusion-consumption problem over
rectangular shaped domain. Label “TC” = tissue compartment; label “BC##”=boundary
compartment of type ##; label “IC” = inactivated compartment.

Consider the boundary compartment on the middle of the left margin, labeled as

“BC21”, which has oxygen concentration Cn

2,1 at time tn = ∆t · n. Based on equation

3.59, its oxygen concentration in the next time step is

Cn+1

2,1 = max

ΣnbCn

2,1 − 4Cn

2,1

(cid:26)

(cid:104)

Cn

2,1 +

∆t

(∆x)2

67

(cid:27)
(cid:105) − ∆t · κ, 0

,

(3.60)

where

ΣnbCn

2,1 = Cn

1,1 + Cn

2,0 + Cn

3,1 + Cn

2,2.

(3.61)

Because the (2, 0)th compartment is an inactivated compartment, its value equals to the
2,0 ≡ Cn

(2, 2)th compartment for all times so that no flux crosses the boundary, i.e. Cn

2,2

for all n. Hence, the difference equation gives the oxygen concentration for boundary

compartments labeled as “BC21” is in general

(cid:26)

(cid:104)

Cn+1

i.j = max

Cn

i,j +

∆t

(∆x)2

Cn
i−1,j + Cn

i+1,j + 2Cn

i,j+1 − 4Cn

i,j

(cid:27)
(cid:105) − ∆t · κ, 0

.

(3.62)

Governing difference equations for boundary compartments labeled as “BC22”, “BC31”

and “BC32” in table 3.10 can be derived in the same way.

Now, consider the boundary compartment on the left up corner, labeled as “BC41”,
whose oxygen concentration at time tn = ∆t · n is Cn

1,1. Based on equation 3.59, its

oxygen concentration in the next time step is

(cid:26)

(cid:104)

Cn

1,1 +

∆t

(∆x)2

ΣnbCn

1,1 − 4Cn

1,1

(cid:27)
(cid:105) − ∆t · κ, 0

,

(3.63)

(3.64)

Cn+1

1,1 = max

where

ΣnbCn

1,1 = Cn

0,1 + Cn

1,0 + Cn

2,1 + Cn

1,2.

However, both the (0, 1)th and the (1, 0)th compartments are inactivated compart-

ments, which means they are outside the boundary. Due to the reflecting boundary

condition, Cn

0,1 is assumed to be identical to Cn

2,1 and Cn

1,0 is assumed to be identical

to Cn

1,2 for all n. Therefore, the difference equation updates the oxygen concentration

68

for a compartments labeled as ‘BC40” is generally

(cid:26)

(cid:104)

(cid:16)

Cn+1

i,j = max

Cn

i,j +

∆t

(∆x)2

2

i+1,j + Cn
Cn

i,j+1

(cid:17) − 4Cn

i,j

(cid:27)
(cid:105) − ∆t · κ, 0

.

(3.65)

Governing difference equations for boundary compartments labeled as “BC50”, “BC60”

and “BC70” in table 3.10 can be derived in the same way.

IC

IC

BC31 BC51

IC

IC

TC BC52

IC

IC

IC

IC

IC

IC

TC

TC

TC BC53 IC

TC BC22 IC

TC

TC

TC

Table 3.11: Discrete compartmental map of 2D unsteady diffusion-consumption problem over
irregular shaped domain. Label “TC” = tissue compartment; label “BC##”=boundary
compartment of type ##; label “IC” = inactivated compartment.

For circular tissue domain or any other irregular shaped domains, a discretization

method is used. Table 3.11 shows a local compartment map over a piece of curved

boundary. Connecting the centers of boundary compartments “BC51”, “BC52” and

“BC53” yields a discretized boundary. Consider the boundary labeled as “BC52” with

oxygen concentration Cn
i,j at time tn. A reflecting boundary condition implies that,
towards this compartment, oxygen concentration of the (i−1, j)th (top) inactive neigh-
bor is identical to that of the (i, j−1)th (left) tissue neighbor and oxygen concentration

of the (i, j + 1)th (right) inactive neighbor is equals to that of the (i + 1, j)th (bottom)

tissue neighbor. Therefore, the difference equation updates the oxygen concentration

69

Similarly, the difference equation for compartment “BC51” is

Cn+1

i,j = max

Cn

i,j +

∆t

(∆x)2

Cn
i,j−1 + 3Cn

i+1,j − 4Cn

i,j

(cid:104)

(cid:27)
(cid:105) − ∆t · κ, 0

.

(3.67)

(cid:26)

(cid:26)

(cid:26)

of this boundary compartment is

Cn+1

i.j = max

Cn

i,j +

∆t

(∆x)2

(cid:104)

2Cn

i,j−1 + 2Cn

i+1,j − 4Cn

i,j

(cid:27)
(cid:105) − ∆t · κ, 0

.

(3.66)

and for compartment “BC53” is

Cn+1

i,j = max

Cn

i,j +

∆t

(∆x)2

(cid:104)

3Cn

i,j−1 + Cn

i+1,j − 4Cn

i,j

(cid:27)
(cid:105) − ∆t · κ, 0

.

(3.68)

3.2.2 Two capillaries of uneven strength in a circular domain

Consider, after normalization, a circular region of unit radius contains two capillaries at

arbitrary locations and different diffusion strengths. Let qm be the rate of amount of oxygen

(per length of capillary) out of the mth capillary, with radius Rc centered at r = rm and

θ = θm. Let q1 = 2/3, r1 = 0.5, θ = 0, q2 = 1/3, r2 = 0.5, θ = π. The boundary and initial

conditions are that there is no flux on the circular boundary and the oxygen concentration

is zero everywhere over the tissue region when t = 0. The consumption is κ = 1 per volume
or whatever is available. As t → ∞, the transient solution should approach the steady state

solution

− N(cid:88)

j=1

qj
2

C =

r2
4

∞(cid:88)

n=0

ln pj +

rn (An cos nθ + Bn sin nθ) ,

(3.69)

70

where

and

N(cid:88)

j=1

qjam

j cos(mαj),

qjam

j sin(mαj), m (cid:62) 1,

An =

Bn =

1
2m

1
2m

N(cid:88)

j=1

(3.70)

(3.71)

as found analytically by Wang [68] Fig. 3.12.

Here we used ∆x = 0.01 and ∆t = 0.00002. The supply boundary is found by drawing

the flux lines and searching for curves that starting from any point close to the saddle point

and end close to a node point. Figure 3.10 shows the transient concentration profiles at

t = 0.0002, t = 0.003, t = 0.006, t = 0.009, t = 0.012, t = 0.018, t = 0.024 and t = 0.048. As

the figure shows, the majority of the tissue region is nourished in a very short time, but for

region away from both capillaries, it takes much more time for oxygen to be delivered. Those

regions are in great danger if the capillary supply was decreased or the tissue consumption

was increased. Actually, it takes t = 6.64 for the entire region to be fully nourished and

for the concentration profile to reach the steady state. Figure 3.12 shows the comparison of

the concentration profile at steady state with the analytical solution found by Wang [68].

Our DCM is extremely accurate as the concentration level curve agrees perfectly with the

analytic solution even near the circular boundary.

71

t = 0.0002

t = 0.003

t = 0.006

t = 0.009

Figure 3.10: Transient process of oxygen diffusion from two capillaries of different strengths
into a circular consuming region. Dark regions show the zero oxygen area. The vertical and
horizontal axes shows the compartment numbers used in each direction.

72

t = 0.012

t = 0.018

t = 0.024

t = 0.048

Figure 3.11: Transient process of oxygen diffusion from two capillaries of different strengths
into a circular consuming region. Dark regions show the zero oxygen area. The vertical and
horizontal axes shows the compartment numbers used in each direction.

73

Figure 3.12: Concentration profile of two unequal strength capillaries in a circular domain,
comparing to the analytic solution[68].

In addition, Figure 3.13 shows the plot of streamlines and the supply boundary curve

found by searching the flux line that starts from a saddle point and end at a node point.

Our result again agrees perfectly with those found analytically by Wang [68].

Figure 3.13: Flux and supply boundary of two capillaries of unequal strength in a circular
domain, comparing to the analytic solution [68].

3.2.3 Nine capillaries of uneven strength in a rectangular domain

Consider nine capillaries of uneven strength in a squared region of unit side length. Table

3.12 shows the randomly generated location and flux strength corresponding to each capillary

[62]. The boundary and initial conditions are that there is no flux on the circular boundary

and the oxygen concentration is zero everywhere over the tissue region when t = 0. The

74

consumption is κ = 1 per volume or whatever is available. The sum of flux equals to 1 here,
so that the tissue region is guaranteed be fully perfused at the steady state. As t → ∞, the

transient solution should approach the steady state solution. Note that, DCM can easily

deal with tissue oxygen deficiency cases. The purpose of this example and the above one is

to validate DCM with existing results.

m xm

ym

flux strength of the mth capillary, Qm

1

2

3

4

5

6

7

8

9

0.501

0.498

0.721

0.235

0.242

0.754

0.091

0.046

0.633

0.932

0.867

0.354

0.425

0.301

0.664

0.316

0.723

0.833

0.044

0.089

0.133

0.133

0.089

0.023

0.133

0.178

0.178

Table 3.12: Locations and flux strength of nine capillaries used in the example.

Following the discrete compartmental method, computations were carried out with ∆x =

0.01 and ∆t = 0.0000166. Figure 3.14 shows the oxygen concentration level curves at steady

state, where the tissue domain is fully nourished as expected. Our result agrees well with the

analytic result found by Sun [62]. Actually, DCM is more accurate near the tissue boundary

comparing to the analytic solution using a matching technique with first 40 terms in each

set of coefficients, as the level curves are more perpendicular to the domain boundary.

75

Figure 3.14: Oxygen concentration level curves of nine capillaries with different strengths
at steady state. The vertical and horizontal axes shows the compartment numbers used in
each direction.

Remark: If capillary flux strength is given instead of concentration, in the difference

equation for capillary compartments (3.49), we simply replace the concentration difference

in the corresponding direction by the flux strength. For example, for a tissue compartment

neighboring the mht capillary, say its (i + 1, j) neighbor is a capillary compartment, then

the governing difference equation for that tissue cell is

(cid:104)

C

n+
i,j = Cn

i.j +

∆t

(∆x)2

Cn
i−1,j + Cn

i,j−1 + Cn

(cid:105)

.

(3.72)

i,j+1 − 3Cn

i,j + Qm

76

0.0010.0010.010.010.020.020.030.030.030.040.040.040.050.050.050.060.060.060.060.070.070.080.080.090.09������������������������Discussion

In this section, we constructed the discrete compartment model in two-dimensional tissue

region with multiple capillaries of uneven location and different strength. Two examples are

demonstrated and the result at steady state agrees perfectly with the analytical solutions.

It is shown by the steady state solution that oxygen coming out from stronger capillary

will diffuse further and help nourishing low oxygen tension area. In addition to that, the

transient solutions show that the time taken for oxygen to be delivered to region away from

the capillaries could be considerable. However, a 2D model can not be used to explain the

real biological phenomena, since the real capillary-tissue region is 3D. That means there are

more pathways for oxygen to move and the oxygen concentration profile changes from the

arterial end to the veinous end. In order to extend the DCM to 3D, we will first introduce

the sliding element method and construct the 2D discrete compartmental model with the

convection of oxygen inside the capillary taking into account.

3.3 Two Dimensional Unsteady Diffusion-Consumption-

Convection

For oxygen transport to tissue, the convection due to blood flow cannot be ignored. We

consider the convection of oxygen in capillary together with diffusion and consumption in

tissue. Starting from a two-dimensional problem, introduce the sliding fluid element method

on the discrete compartment model inside capillary compartments and combine it with the

diffusion-consumption process in the tissue compartments. The three-dimensional transient

problem is introduced in the following chapter..

77

Figure 3.15: Schematic diagram for 2D capillary-tissue exchange with convection, diffusion
and consumption with nondimensionalized parameters. vb is the blood velocity; D = 1 is
the normalized diffusion coefficient; κ is the tissue consumption coefficient; L is the length
of the tissue region along the z-axis; W is the thickness of the tissue region along the x-axis;
oxygen concentration on the arterial end of the capillary is normalized by itself and therefore
equals to 1.

Consider blood with velocity v(cid:48)

b inside a long and straight capillary with radius R(cid:48)
centered along the z(cid:48)-axes, where substrates and oxygen permeate through the capillary
membrane wall with permeability P , diffuse into the nearby tissue region of length L(cid:48) and
thickness W(cid:48) with diffusivety D(cid:48), and was consumed there at constant rate κ(cid:48). Assume

c

solutes are radially well mixed inside the capillary and the axial diffusion rate inside the
capillary is neglected in comparison to transverse diffusion. Let C(cid:48)
t(x(cid:48), z(cid:48), t(cid:48))
be the oxygen concentration inside the capillary and tissue respectively, with C(cid:48)
c being non-
negative for any x(cid:48), z(cid:48), and t(cid:48). The governing partial differential equations to this unsteady

c(z(cid:48), t(cid:48)) and C(cid:48)

convection-diffusion-consumption problem are

(cid:33)

− κ(cid:48),

(3.73)

(cid:32)

∂C(cid:48)
∂t(cid:48) = D(cid:48)

t

∂2C(cid:48)
∂x(cid:48)2 +

t

∂2C(cid:48)
∂z(cid:48)2

t

78

κD=1CAPILLARYTISSUELW=1xz1vbin the tissue region, and

∂C(cid:48)
∂t(cid:48) = −v(cid:48)

b · ∂C(cid:48)
∂z(cid:48)

c

c

(3.74)

inside the capillary, where 0 (cid:54) x(cid:48) (cid:54) W(cid:48), 0 (cid:54) z(cid:48) (cid:54) L(cid:48) and t(cid:48) (cid:62) 0. Assume there is no flux

on the boundary of tissue region. The permeability H is infinite for gas solute like oxygen.

Thus the capillary and tissue oxygen concentration are equivalent on the capillary wall at
x(cid:48) = R(cid:48)
c.

Normalize all lengths by the thickness of the region W(cid:48), concentration by that of the
arterial end Ca, time by W(cid:48)2/D(cid:48) and drop primes. In terms of normalized parameters, the

governing equations and boundary conditions are

(cid:26)(cid:18)∂2Ct

∂x2 +

∂2Ct
∂z2

(cid:19)

(cid:27)

= max

− κ, 0

x ∈ (Rc, 1], z ∈ [0, L], t ∈ [0,∞)

,

∂Ct
∂t
∂Cc
∂t

= −vb · ∂Cc
∂z

,

x ∈ [0, Rc], z ∈ [0, L], t ∈ [0,∞)

Cc(0, t) = 1,

Cc(z, t) = Ct(Rc, z, t),

(cid:12)(cid:12)(cid:12)(cid:12)x=1

∂Ct
∂x

=

∂Ct
∂z

(cid:12)(cid:12)(cid:12)(cid:12)z={0,L}

= 0,

(3.75)

(3.76)

(3.77)

(3.78)

(3.79)

where κ = κ(cid:48)W(cid:48)2/D(cid:48)Ca, vb = W(cid:48)v(cid:48)

b/D(cid:48), Rc = R(cid:48)

c/W(cid:48) and L = L(cid:48)/W(cid:48). This two region and

one barrier model is illustrated in Figure 3.15.

If oxygen deficient occurs, the flux and the concentration value equals zero on the moving

boundaries ∂On

∂C
∂n

(cid:12)(cid:12)(cid:12)(cid:12)(x,z)∈∂On

= C|(x,z)∈∂On = 0

(3.80)

Oxygen deficit region and its moving boundaries are defined in section 2.3.1.

79

3.3.1 Formulation: The Space and Time Step Matching Technique

Consider a two dimensional discrete compartmental domain of ordered identical cube com-

partments with equal side length ∆x and cross-section area A. Since oxygen is assumed to

be radially well mixed inside the capillary, the choice of ∆x is arbitrary in the 2D problems.

The strategies for diffusion and consumption in the tissue region are the same as those in

section 3.2. In this section, we shall focus on developing the discrete compartmental model

for convection.

The fundamental strategy in this approach is to match the extravascular diffusion to

the intravascular convection, so that the oxygen flux changes in both regions are adequately

accounted for independent of the choice of compartment size. Figure 3.16 demonstrates the

main idea and basic steps of this blood-tissue exchange algorithm. With the chosen ∆x,

the capillary is divided into N cube compartments of equal dimensions, indexed k = 1, N
beginning at the inflow and the capillary plasma transit time t = ∆t · n. Now, we further

“split” the time step interval Let

∆t+ = begining of the time interval

∆t− = end of the time interval

(3.81)

and hence ∆t = ∆t− − ∆t+. Before the time begins, the strength of the capillary is given

by specifying the concentration on the arterial end. At beginning of each time interval ∆t+,

the contents in each capillary compartment slides one segment downstream such that the
concentration in the kth compartment is replaced by that in the (k − 1)th, at which moment

change in the tissue region is frozen. Thereafter, during the time interval ∆t, diffusion and

consumption occur in the tissue region. At the end of the current time interval ∆t−, the

80

exchange pauses again and oxygen inside the capillary compartments slide one grid down for

the next time step. The outflow at time t+ would be the concentration of the N th capillary

compartment at time t−.

Notice that, over each time interval, the sliding process in the capillary takes infinites-

imally small amount of time. The reason is that the oxygen concentration in the capillary

compartments is equivalent to the boundary condition for the concentration changes in the

tissue region. On the other hand, the diffusion-consumption process takes the entire interval

of time over each time step, ∆t, which is determined by matching the time and space steps

so that contents inside the capillary compartments slide exactly one segment at a time.

Figure 3.16: Lagrangian sliding fluid element method. With each ∆t, oxygen in the capillary
compartment slide instantaneously one segment downstream, after which the diffusion and
consumption are computed.

81

  k - 1          k         k + 1   ART t = 0+ t = Δt- t = Δt+ t = 2Δt- t = 2Δt+   ARTVEINVEINVEINVEINVEIN   ART   ART   ARTDiffuse + ConsumeDiffuse + ConsumeSlideSlideSlide3.3.2 A Simplified Model: Unsteady

Convection-Permeation-Consumption

Consider a problem with single layer of tissue region and no diffusion in the longitudinal di-

rection in order to focus on developing the discrete compartmental algorithm for convection-

permeation-consumption problem based on the sliding element method. This simplified

model serves as a basis for more complex models.

Consider a 1-barrier, 2-region model accounting for oxygen convection along a capillary

with velocity vb. Oxygen permeates into adjacent cell regions and was consumed there at

a constant rate κ. Let Cc(z, t) and Ct(z, t) be the oxygen concentration in capillary and

tissue region respectively. In non-dimensional form, this unsteady convection-permeation-

consumption problem can be expressed using partial differential equations as

∂Cc
∂t
∂Ct
∂t

+ (Ct − Cz) ,

= −vb · ∂Cc
∂z
= (Cc − Ct) − κ.

(3.82)

(3.83)

Figure 3.17: Compartment model for capillary-tissue exchange in well-perfused organ.

Now we solve this problem from a compartmental point of view. Figure 3.17 shows the

discrete compartmentalized capillary-tissue region. Oxygen is well-perfused in the tissue

82

0,k0,k+10,k-11,kACAPILLARYTISSUEΔxκ01,k-11,k+1vbregion without longitudinal diffusion. With a selected space step ∆x, the corresponding

sliding period is determined by

p =

∆x
vb

.

(3.84)

Table 3.13 shows a sample map of discrete compartmental modeling used for this prob-

lem, where interior capillary compartments, the capillary compartment on the arterial end,

the capillary compartment on the venous end, tissue compartments and inactivated compart-

ments are distinguished. Arterial capillary compartment shall be assigned an initial value at

the beginning and can be changed at anytime during the process. Depending on the size of

p, we classify the computational process into a basic case and a general case.

83

The basic case: p (cid:54) (∆x)2/2d, where d is the space dimension of the problem. This

usually happens when ∆x or vb is relatively large. In such case, we set

∆t = p.

(3.85)

and the sliding and diffusing process will happen alternatively, switching back and forth.

As time begins, contents in the capillary slides down stream one grid a time.

In the

beginning time step ∆t+

C

n+
0,k = C

n−
0,k−1.

(3.86)

Then, oxygen permeates into the tissue region due to the partial pressure difference and

consumption κ (or whatever remains) over the time step ∆t is

(cid:110)

1,k + ∆t ·(cid:16)

n+

C

n+

0,k − C

C

(cid:17)
1,k − κ

n+

(cid:111)

, 0

.

(3.87)

(n+1)−
1,k

C

= max

Combining the above two equations and drop “±” subscripts, the governing difference equa-

tion of concentration change of tissue compartments over one complete time interval ∆t

(cid:110)

1,k + ∆t ·(cid:16)

Cn

(cid:17)
1,k − κ

(cid:111)

0,k−1 − Cn
Cn

Cn+1

1,k = max

, 0

.

(3.88)

During the same time interval, capillary compartments lost an equivalent amount of oxygen

which permeated into tissue cells. Balancing the inflow and outflow, concentration inside a

capillary compartment at the end of the previous time step is

0,k + ∆t ·(cid:16)

n+

(n+1)−
0,k

C

= C

n+

1,k − C

n+
0,k

C

(cid:17)

.

(3.89)

84

Associate the above with Equation 3.86 and drop “±” subscripts gives the governing differ-

ence equation for capillary compartments

0,k−1 + ∆t ·(cid:16)

(cid:17)

.

1,k − Cn
Cn

0,k−1

(3.90)

Cn+1
0,k = Cn

At the beginning of the next time step ∆t+, the diffusion-consumption pause for a moment

and oxygen inside the capillary compartments moves another whole block forward.

CCA CC

CC

CC

CC CCV

IC

STC STC STC STC

IC

Table 3.13: Discrete compartmental map of unsteady convection-permeation-consumption
problem. Label “CC”=capillary compartment; label “CCA”=capillary compartment on the
arterial end; label “CCV”=capillary compartment on the venous end; label “STC” = single
layer tissue compartment; label “IC” = inactivated compartment.

The general case: p > (∆x)2/2d happens when one or both of ∆x and vb are relatively

small. In such case, we let

∆t =

p
b

,

(3.91)

where b is the least positive whole number such that p/b (cid:54) (∆x)2/2d. Then, Eq. 3.86 is

called for every b steps during the computation and when b = 1 it covers the basic case.

Note that, using extra sliding step will not cause convergence problems but will result in

higher concentration values, which means the computational errors may be larger.

85

Example 3.3.1. As an illustrative example, consider a capillary of length L = 1 and blood

velocity vb = 1 adjacent to tissue cells that consumes oxygen a constant rate κ = 0.2.

Initially the concentrations are zero everywhere until a constant unit amount of oxygen is

suddenlyinfused at the arterial end of the capillary, or Cc(0, t) = 1. Assume the outer tissue

boundary is reflecting, which means that there is no flux from or to the neighboring regions.

Using the discrete compartment model with the sliding element technique, the calculation is

carried with ∆x = 0.01 and ∆t = 0.01. Figure 3.18 shows the numerical results of transient
concentration profiles in both regions at t = 0.5, t = 1, t = 5, and t = 10. As t → ∞, the

unsteady solution approaches

· z,
Cc(z) = 1 − κ
v
Ct(z) = 1 − κ − κ
v

· z,

(3.92)

(3.93)

which is the steady state solution.

86

t=0.5

t=1

t=5

t=10

Figure 3.18: Oxygen convection into a consuming region using a two region model. Red dots
represent oxygen concentration flowing in capillary, blue dots represent oxygen concentration
in adjacent tissue region with constant consumption and solid lines are analytic solutions at
steady state.

As the figure shows, under the given velocity vb = 1, blood slides to the middle of

the capillary at t = 0.5 and slides to the veinous end of the capillary at t = 1. Oxygen

concentration in the tissue region gradually increases with the movement of blood in the

capillary. Oxygen in the capillary and tissues on the arterial end first reaches the steady

state solution with a difference of 0.2, which is the given consumption rate κ. At t = 5,

oxygen concentration in both regions gradually approaches the stead state solution and the

their difference gradually approaches κ = 0.2. At t = 10, oxygen concentration in both

regions reaches the steady state solution. The numerical solution using DCM with the

sliding element method agrees perfectly with the analytic solution.

87

������������������������������������������������������������������������������������������������������������������������3.3.3 Unsteady Convection-Diffusion-Consumption

Now add oxygen diffusion within the tissue and develop the discrete compartmental model.

The partial differential equations for the unsteady convection-diffusion-consumption problem

are derived at the beginning of this section 3.3.

Consider a two space dimensional discrete compartmental space of ordered identical cube

compartments with equal side length ∆x and cross-section area A(Fig. 3.19). Differing

from the diffusion-consumption problem in section 3.2, the choice of ∆x is arbitrary in 2D

convection-diffusion-consumption model with a single capillary. Since we assumed oxygen

to be well mixed inside capillary, one single row of compartments is adequate to represent

the capillary region. Therefore, compartments with side length of any reasonably small

∆x can be used but there lower surfaces have to be matched perfectly with the capillary

compartments. On the other hand, it is advantageous to choose a ∆x that divide the tissue

region exactly to reduce the discretization error. With a fixed ∆x, the corresponding time

step ∆t is determined by equation 2.5.

Figure 3.19: Compartment model for capillary-tissue exchange with diffusion.

Table 3.14 shows a sample Map for this problem, in addition to the problem with two

layers, tissue compartments are further divided. Capillary adjacent tissue compartments are

88

i,ki,k+1i,k-1i+1,kACAPILLARYTISSUEDκ0i-1,k-1vbi-1,ki-1,k+1marked out differently as wells as the different types of boundary compartments. At the

venous end, one extra compartment is added in order to quantify the oxygen concentration

through the outflow. The governing difference equations for a tissue compartment and the

boundary tissue compartments not adjacent to the capillary are identical to the compart-

ments with the same label in section 3.2. Here, we focus on deriving the difference equation

for the capillary compartment and neighboring capillary compartments.

CCA CC

CC

CC

CC CCV

IC

IC

IC

BC11 CTC CTC BC12

BC21

TC

TC BC22

BC60 BC32 BC32 BC70

IC

IC

IC

Table 3.14: Discrete compartmental map of 2D unsteady convection-diffusion-consumption
problem. Label “CC”=capillary compartment;
label “CCA”=capillary compartment on
label “TC”
the arterial end;
= tissue compartment;
label
“BC##”=boundary compartment of type ##; label “IC” = inactivated compartment.

label “CTC” = capillary neighboring tissue compartment;

label “CCV”=capillary compartment on the venous end;

Similar to the two layer model, a given concentration value from blood on the arterial

end of the capillary is initially stored in the source compartment

Label “CCA”: C0

i,k = CA.

(3.94)

1. Difference Equation of Capillary Compartments:

As time begins, oxygen inside the capillary compartments slides instantaneous down

stream at the beginning of time step ∆t+

C

n+
i,k = C

n−
i,k−1.

89

(3.95)

From their to the end of time step ∆t−, oxygen diffuses into the neighboring tissue

regions driven by the partial pressure difference. Balancing the mass or concentration

in a capillary cell at the end of each time step gives

(n+1)−
i,k

C

= C

n+
i,k +

∆t

(∆x)2

n+

i+1,k − C

n+
i,k

C

(cid:16)

(cid:17)

.

(3.96)

Substitute into Eq. (3.95), drop subscripts and simplify give the differential equation

for capillary compartments

Label “CC”: Cn+1

i,k = Cn

i,k−1 +

(cid:104)

∆t

(∆x)2

(cid:105)

.

i+1,k − Cn
Cn

i,k−1

The outflow at the end of each time step is found by

Label “CCV”: Cn+1

i,k = Cn

i,k−1

(3.97)

(3.98)

2. Difference Equation of Capillary Neighboring Tissue Compartments:

Oxygen diffused from the capillary continuously exchange among the tissue compart-

ments due to the partial pressure difference and κ or consumed within the tissue.

Similarly to section 3.2, concentration in a tissue compartment at the end of each time

step can be find as

(n+1)−
i.k

C

= max

(cid:26)

C

n+
i,k +

∆t

(∆x)2

(cid:104)

ΣnbC

n+

i,k − 4C

n+
i,k

(cid:27)
(cid:105) − ∆t · κ, 0

,

(3.99)

where

ΣnbC

n+
i,k = C

n+
i−1,k + C

n+
i,k−1 + C

n+
i+1,k + C

n+
i,k+1.

(3.100)

90

In the capillary, oxygen slide from the (i−1, k−1)th capillary neighbor to the (i−1, k)th

capillary neighbor and diffuse into the tissue during each time interval. Thus, for a

tissue compartment adjacent to a capillary, the infusion of oxygen for time interval

∆t+ can be expressed as

C

n+
i−1,k = C

n−
i−1,k−1.

(3.101)

Substitute the above relation into Eq. (3.99) and drop subscripts on t. The difference

equation for concentration tissue compartment adjacent to a capillary is

(cid:26)

(cid:104)

∆t

(∆x)2

Label “CTC”: Cn+1

i.k = max

Cn

i,k +

Cn
i−1,k−1 + Cn

+Cn

i,k−1 + Cn

i,k+1 − 4Cn

i,k

i+1,k

(cid:105) − ∆t · κ, 0
(cid:111)

(3.102)

,

3. Difference Equation of Boundary Compartments Adjacent to a Capillary

Combining the technique above for tissue compartments adjacent to a capillary and the

compartmental treatment on the reflecting boundary in section 3.2, difference equations

for tissue compartments adjacent to a capillary are

Label “BC11”: Cn+1

i,k = max

(cid:26)

Cn

i,k +

∆t

(∆x)2

Cn
i−1,k−1
i,k+1 − 4Cn

i,k

(cid:105) − ∆t · κ, 0
(cid:111)

,

(3.103)

+Cn

i+1,k + 2Cn

(cid:104)

(cid:104)

and

Label “BC12”: Cn+1

i,k = max

(cid:26)

Cn

i,k +

∆t

(∆x)2

Cn
i−1,k−1
i,k−1 − 4Cn

i,k

(cid:105) − ∆t · κ, 0
(cid:111).

(3.104)

+Cn

i+1,k + 2Cn

91

Example 3.3.2. Consider a capillary-tissue domain as diagrammed in Figure 3.15. As an

illustrative example, a single capillary of length L = 40, radius Rc = 0.1 and blood velocity

vb = 20 is adjacent to a rectangular tissue region of width W = 1. Oxygen concentration on

the arterial end is Ca = 1, the diffusion constant is D = 2 and the consumption coefficient is

κ = 0.1. Assume the outer tissue boundary is reflecting, which means that there is no flux

from or to the neighboring regions.

We aim to compute the steady state solution from the transient state. Assume there are

initially zero oxygen everywhere and at starting time blood with oxygen suddenly flows in.

Our numerical solution is computed under of ∆x = 0.05 and ∆t = 0.0005. The space step

size is determined by the capillary radius, such that the Rc is an integral multiple of ∆x

and the time step size is determined following the convergence theorem (see section 5 for the

theorem and proof).

A transient solution is calculated until reaching steady state. Figure 3.20 shows the

transient oxygen concentration profiles in the tissue region at t=0.1, t=0.25, t=1, t=2.5,

t=5, t=7.5 and t=10. As the figure shows, under the given velocity vb = 20, blood slides

around 40 capillary compartments down stream at t = 0.1. Different from the previous simple

model, where only one layer of tissue compartment is considered, movement of the plasma

inside the capillary slows down while flowing further to the veinous end. This is because the

tissue region is thicker and consumes more oxygen. Hence more oxygen diffuses out to the

tissue and less oxygen lefts inside the capillary. At t = 10, the entire tissue is perfused with

oxygen. Horizontally, the concentration is higher on the arterial end and gradually decreases

towards the veinous end. Vertically, the concentration is higher over tissue regions near the

capillary and slightly lower over tissue regions away from the capillary.

92

Figure 3.20: Concentration curves of a single capillary with blood velocity at t=0.1, t=0.25,
t=1, t=2.5, t=5, t=7.5 and t=10. The vertical and horizontal axes shows the compartment
numbers used in each direction.

93

0.10.20.30.4��������������������������������������0.10.20.30.40.5��������������������������������������0.10.20.30.40.5��������������������������������������0.10.20.40.50.60.70.80.90.3��������������������������������������0.10.20.30.40.50.60.70.80.9��������������������������������������0.10.20.30.40.50.60.70.80.9��������������������������������������0.20.30.40.50.60.70.80.9��������������������������������������Figure 3.21: Concentration profile on the capillary boundary and the tissue boundary, com-
pared with the analytic solution. Circles are numerical solution of oxygen concentration on
the tissue boundary. Triangles are numerical solution of oxygen concentration on the capil-
lary wall. Solid lines represents the oxygen concentration profile obtained from the analytic
solution. Horizontal axis shows the normalized capillary length and vertical axis show the
normalized concentration value.

At t = 15, the unsteady solution approaches the steady state solution, which is obtained

from solving the differential equations over the capillary region and the tissue region. Figure

3.21 shows the comparison of the numerical solution and the analytic solution of oxygen

concentration on the tissue boundary and on the capillary wall. The numerical solution

using DCM with the sliding element method agrees perfectly with the analytic solution.

Discussion

In this section, we constructed 2D discrete compartment model with the convection of oxygen

inside the capillary taking into account. Two examples are demonstrated and the result at

steady state agrees perfectly with the analytical solutions. It is shown by the steady state

solution that the oxygen concentration inside the capillary decreases from the arterial end

to the veinous end. In addition to that, the transient solutions show that the velocity of

oxygen (not plasma) inside the capillary may not be constant or even linear, but depends

on the oxygen concentration of the surrounded tissues. Again, a 2D model can not be used

94

○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△����������������������to explain the real biological phenomena. We will construct the 3D discrete compartmental

model in the following chapter.

95

Chapter 4

The Three Dimensional Modeling

The problem of capillary-tissue supply is unlikely to be limited to 2D tissue cross section

area, 3D tissue cylinder with a single capillary or a 3D steady state. None of these could

represent the real, especially transient physiological phenomena. In this chapter, we extend

the techniques of the previous chapters and develop the 3D discrete compartmental model for

simulating the transient process of multi-capillary oxygen supply to striated muscle tissue.

4.1 Three-Dimensional Unsteady Diffusion-Consumption-

Convection

The problem of three-dimensional unsteady convection, diffusion and consumption problem

describing the oxygen transfer from capillary to tissue can be generally stated as follows.
Consider M straight capillaries of same radius R(cid:48)
of z(cid:48), but different oxygen concentration C(cid:48)
m, where m = 1, 2, ..., M , located irregularly in
a bounded and non-simply connected tissue region Ω(cid:48) ∈ R3. Define individual capillary

c, same blood flow velocity v(cid:48)

b in the direction

domain

Ω(cid:48)
m =

(cid:110)
(x(cid:48), y(cid:48), z(cid:48)) : (x(cid:48), y(cid:48), z(cid:48)) ∈ Ω(cid:48) and (x(cid:48) − x(cid:48)

m)2 + (y(cid:48) − y(cid:48)

m)2 (cid:54) (R(cid:48)

c)2(cid:111)

,

(4.1)

96

where (x(cid:48)

m, y(cid:48)

m, z(cid:48)) are the coordinates of points along the center line of each capillary. Thus,

the total capillary domain is

and the absorbing tissue domain is

M(cid:91)

m=1

Ω(cid:48)

m

Ω(cid:48)
c =

Ω(cid:48)
t = Ω(cid:48) \ Ω(cid:48)
c.

(4.2)

(4.3)

Assume oxygen permeate through the thin capillary membrane wall without resistance and
diffuse into the nearby tissue region at a diffusion rate D(cid:48) and is consumed in the tissue at
a rate of κ(cid:48). Let C(cid:48)
t(x(cid:48), y(cid:48), z(cid:48), t(cid:48)) be the function of non-negative oxygen concentration of the
c(z(cid:48), t(cid:48)) the non-negative oxygen function in the capillaries, the governing
tissue region and C(cid:48)
(cid:32)

partial differential equations to this unsteady convection-diffusion-consumption problem are

(cid:33)

− κ(cid:48),

(x(cid:48), y(cid:48), z(cid:48)) ∈ Ω(cid:48)

t, t(cid:48) ∈ [0,∞),

∂C(cid:48)
∂t(cid:48) = D(cid:48)

t

∂2Ct
∂x(cid:48)2 +

∂2C(cid:48)
∂y(cid:48)2 +

t

∂2C(cid:48)
∂z(cid:48)2

t

(4.4)

(4.5)

(4.6)

and

∂C(cid:48)
∂t(cid:48) = −v(cid:48)

b · ∂C(cid:48)
∂z(cid:48) ,

c

c

(x(cid:48), y(cid:48), z(cid:48)) ∈ Ω(cid:48)

c, t(cid:48) ∈ [0,∞).

On the outer tissue boundary ∂Ω(cid:48) there is no flux or the normal derivative of concentration

is zero

∂C(cid:48)
∂n

(cid:12)(cid:12)(cid:12)(cid:12)(x(cid:48),y(cid:48),z(cid:48))∈∂Ω(cid:48) = 0.

97

On the arterial end of each capillary(z(cid:48) = 0), the concentration is C(cid:48)

m, differing in capillary

strength

C(cid:48)(x(cid:48), y(cid:48), 0, t(cid:48)) = C(cid:48)
m,

(cid:0)x(cid:48), y(cid:48), 0(cid:1) ∈ Ω(cid:48)

m, t(cid:48) ∈ [0,∞) , m = 1,··· , M.

(4.7)

Since the permeability through the capillary wall can be regarded as infinite for gas solutes

like oxygen, the boundary conditions are

t(x(cid:48), y(cid:48), z(cid:48), t(cid:48)) = C(cid:48)
C(cid:48)

c(z(cid:48), t(cid:48)),

(x(cid:48), y(cid:48), z(cid:48)) ∈ ∂Ω(cid:48)

c

(4.8)

Normalize all lengths by ˜L =(cid:112)A(Ω(cid:48) ∩ {z = 0}), the square root of the cross sectional

area of Ω(cid:48) on the xy-plane, concentration by the average on the arterial end C(cid:48)
a, time by
˜L2/D(cid:48) and drop primes. In terms of normalized parameters, the governing equations and

boundary conditions are

(cid:18) ∂2Ct

∂Ct
∂t

∂C
∂t

=

∂x2 +
= −vb · ∂C
∂z

,

(cid:19)

− κ,

∂2Ct
∂z2

∂2Ct
∂y2 +
(x, y, z) ∈ Ωc, t ∈ [0,∞)

(x, y, z) ∈ Ωt, t ∈ [0,∞)

Ct(x, y, z, t) = Cc(z, t),

(x, y, z) ∈ ∂Ωc,

C(x, y, 0, t) = Cm,

(cid:12)(cid:12)(cid:12)(cid:12)(x,y,z)∈∂Ω

∂Ct
∂n

= 0,

(x, y, 0) ∈ Ωm, t ∈ [0,∞) , m = 1,··· , M,

(4.9)

(4.10)

(4.11)

(4.12)

(4.13)

where κ = κ(cid:48) ˜L2/D(cid:48)C(cid:48)

a, vb = ˜Lv(cid:48)

b/D(cid:48) and Rc = R(cid:48)

c/ ˜L.

If oxygen deficiency occurs, the flux and the concentration value equals zero on the

98

moving boundaries ∂On

(cid:12)(cid:12)(cid:12)(cid:12)(x,y,z)∈∂On

∂Ct
∂n

= Ct|(x,y,z)∈∂On = 0.

(4.14)

Oxygen deficit region and its moving boundaries are defined in section 2.3.1. Note the

governing equations are 3D and unsteady.

4.1.1 Formulation

Consider a three-dimensional discrete compartmental space of identical cube compartments

with side length ∆x.

It is of much importance to define the cross-section area of each

capillary. Hence, in three dimensional modeling, the choice of ∆x must satisfy

(cid:115)

∆x =

πR2
c
nCC

(4.15)

where nCC = 1, 4, 5, 9, 13, 21..., which is the same as in 2D modeling. Again, shapes of

capillary compartment arrangements need to be approximately radially symmetric.

The formulation of governing difference equations for 3D modeling follows the same

procedure of the 1D and 2D modeling. The (i, j, k)th cell has six adjacent neighbors

(i − 1, j, k), (i + 1, j, k), (i, j − 1, k), (i, j + 1, k), (i, j, k − 1), (i, j, k + 1).

(4.16)

Let Cn

i,j,k denote the oxygen concentration of the (i, j, k)th compartment at time ∆t· n. The

change in mass inside this closed volume V = (∆x)3 over a short time period ∆t is

(cid:104)

(cid:105) · V.

i,j,k − Cn
Cn+1

i,j,k

99

(4.17)

Figure 4.1: A 3D compartment for deriving the difference equations for DCM.

In the direction of the positive z-axis, the oxygen flux across each side from/to the top,

bottom, left, right, back and front neighboring compartments during a time interval ∆t are

(cid:104)
(cid:104)
(cid:104)

1 = −∆x
J n
∆t
3 = −∆x
J n
∆t
5 = −∆x
J n
∆t

i,j,k − Cn
Cn
i,j,k − Cn
Cn
i,j,k − Cn
Cn

i−1,j,k

i,j−1,k

i,j,k−1

(cid:105)
(cid:105)
(cid:105)

,

,

,

2 = −∆x
J n
∆t
4 = −∆x
J n
∆t
6 = −∆x
J n
∆t

i+1,j,k − Cn
Cn
i,j+1,k − Cn
Cn
i,j,k+1 − Cn
Cn

i,j,k

i,j,k

i,j,k

,

,

.

(4.18)

(cid:105)
(cid:105)
(cid:105)

(cid:104)
(cid:104)
(cid:104)

A schematic diagram is given in Fig. 4.1.

1. Difference Equation of Tissue Compartments:

Consider the interior tissue compartments which constantly consume oxygen at the rate

κ. The concentration change caused by the difference between inflows, outflows and

100

consumption in the (i, j, k)th compartment over the time period ∆t can be described

as

[(J n

1 + J n

3 + J n

5 − J n

2 − J n

4 − J n

6 ) · A − κ · V ] · ∆t,

(4.19)

where A = (∆x)2 is the single face area of a cube compartment and V = (∆x)3 the

volume of a cube compartment. Balancing Equation (4.19) with Equation (4.17) gives

D ·(cid:104)

i,j,k − Cn
Cn+1

i.j,k

(cid:34)(cid:32) 6(cid:88)

(cid:105) · V =

(cid:33)

J n
m

m=1

(cid:35)

· A − κ · V

· ∆t.

(4.20)

Substituting fluxes by concentration differences from (4.18) and simplifying yields

(cid:104)

ΣnbCn

i,j,k − 6Cn

i,j,k

(cid:105) − ∆t · κ,

Cn+1
i.j,k = Cn

i.j,k +

∆t

(∆x)2

where

(4.21)

(4.22)

ΣnbCn

i,j,k = Cn

i−1,j + Cn

i+1,j + Cn

i,j−1 + Cn

i,j+1 + Cn

i,j,k−1 + Cn

i,j,k+1.

Now, consider the situation of a local oxygen deficiency, where

(cid:104)

(cid:105)

Cn

i.j,k +

∆t

(∆x)2

ΣnbCn

i,j,k − 6Cn

i,j,k

< ∆t · κ,

(4.23)

in which case the tissue compartment will consume all the available oxygen over that

time step, and thus

Cn+1

i.j,k = 0.

(4.24)

Combining the oxygen sufficient case with the oxygen deficient case, the governing

101

difference equation of an interior tissue compartment is

Cn+1

i.j,k = max

Cn

i,j,k +

∆t

(∆x)2

ΣnbCn

i,j,k − 6Cn

i,j,k

(cid:26)

(cid:104)

(cid:27)
(cid:105) − ∆t · κ, 0

,

(4.25)

where ∆t/(∆x)2 shall be less than 1/6 in the three dimensional problem for stability

and convergence. See Chapter 5 for the theoretical proof of convergence.

2. Difference Equation of Capillary Compartments:

In the three dimensional problem, oxygen inside the capillary compartments not only

diffuse into the tissue region but also flow downstream at the same time. Dealing

with this complex situation, we shall extend the discrete compartmental technique

introduced in sections 3.2 and 3.3. The fundamental idea is to execute the 3D diffusion

and sliding process by a time-space matching technique.

Consider the case ∆x =(cid:112)πR2

c , where each capillary is represented by a single row of

compartments. Each capillary compartment has four tissue neighbors and two capillary

neighbors. The concentration change in the (i, j)th capillary compartment over the

time period ∆t can be described as

[(J n

1 + J n

3 − J n

2 − J n

4 ) · A] · ∆t.

(4.26)

Balancing of mass flow with concentration change in the (i, j, k)th capillary compart-

ment over the time period ∆t yields

D ·(cid:104)

n+

i,j,k − Cn

i,j,k

C

(cid:34)(cid:32) 4(cid:88)

(cid:105) · V =

(cid:33)

(cid:35)

· A

J n
m

· ∆t,

m=1

(4.27)

102

where n+ represents time at the end of the nth time interval. Substitute in J n

1 , J n
2 ,

J n
3 , J n

4 , from (4.18) and simplify. The oxygen concentration in the (i, j, k)th capillary

compartment at the end of the nth time interval is

C

n+
i.j,k = Cn

i,j,k +

∆t

(∆x)2

i,j,k − 4Cn
Cn

i,j,k

 (cid:88)

(i,j)nb

 ,

(4.28)

(4.29)

(4.30)

(cid:88)

(i,j)nb

where

and

Cn
i,j,k = Cn

i−1,j,k + Cn

i,j−1,k + Cn

i+1,j,k + Cn

i,j+1,k.

∆t =

p
b

,

in which p is the sliding period determined by ∆x/vb and b is any positive integer.

Here ∆t is chosen to that give the desired accuracy and also satisfies the condition

of convergence. During each sliding period, the diffusion-consumption is executed for

every ∆t and oxygen in each capillary compartment slide one segment downstream
for every b · ∆t. Thus, the oxygen concentration of a capillary compartment in the

beginning of the next time interval is

Cn+1

i,j,k = C

n+
i,j,k−1.

(4.31)

For the general case ∆x ≤(cid:112)πR2

c , where each capillary is represented by more then

one row of compartments, an average of the oxygen concentration over each capillary’s

cross-section area shall be taken and reassigned after the diffusion process by Eq.

(4.28).

103

3. Difference Equation of Boundary Compartments:

Consider a cuboid discrete compartmental space, where the boundary compartments

can be divided into three categories: 1. tissue compartment with one inactive neigh-

bor(a side boundary), tissue compartment with two inactive neighbors(an edge bound-

ary), 3. tissue compartment with three inactive neighbors(a corner boundary). The

derivation of the governing difference equation for boundary tissue compartments satis-

fying the reflecting boundary condition follows a similar process as in section 3.2.1. For
example, if a boundary compartment on the corner, whose (i − 1, j, k)th, (i, j − 1, k)th
and (i, j, k − 1)th neighbors are inactive compartments, then the difference equation

updating oxygen concentration of this compartment is

(cid:26)

(cid:104)

Cn+1

i.j,k = max

Cn

i,j,k +

∆t

(∆x)2

2Cn

i+1,j,k + 2Cn

i,j+1,k + 2Cn

i,j,k+1 − 6Cn

i,j,k

(cid:27)
(cid:105) − ∆t · κ, 0

.

(4.32)

104

4.1.2 Three capillaries of uneven strength in a cuboid domain

Consider three capillary of uneven strength and random distribution in a cube region, Ω =
W × H × L = 1 × 1 × 1. Locations and plasma oxygen concentration on the arterial end of

the capillaries are randomly generated and shown in table 4.1. The uniform capillary radius

Rc is set to be 0.056419 for demonstration, the blood velocity vb = 180, and the tissue

consumption coefficient κ = 1.

xm

ym

Ca of the mth capillary

0.679

0.390

0.258

0.299

0.511

0.605

0.311

0.186

0.473

Table 4.1: Locations and plasma oxygen concentration on the arterial end of three randomly
distributed capillaries used in the example.

Assume this cubic multi-capillary tissue domain initially has zero oxygen everywhere.

Blood with oxygen suddenly flows into the capillaries and the oxygen starts to diffuse into

the surrounding tissue. Figure 4.2 shows the computed oxygen concentration profiles at time

t = 0.1, t = 0.2, t = 0.5 and t = 2.2. The domain is gradually filled with oxygen and finally

reaching a steady state, which is determined by comparing the concentration values of each

compartment over a constant time period and see if the difference of every corresponding

compartment is less then a given threshold.

105

t = 0.1

t = 0.2

t = 0.5

t = 2.2

Figure 4.2: Concentration profile of three capillaries of uneven strength in a cube domain at
t = 0.1, t = 0.2, t = 0.5 and t = 2.2 (steady state).

106

Of interest is time when the domain becomes fully perfused, tf , which is found by search-

ing for zero-concentration compartments over each time step. If every compartment has a

nonzero value, we say the domain is fully perfused; if the tissue domain has both zero and

nonzero compartments, we say the domain is partially perfused. The first time when the

domain changed from partially perfused to fully perfused is the domain fully perfused time.

Another subject of interest is the time for the transient solution of the whole domain

to reach a steady state, ts. Since between two adjacent time steps the difference between

the concentration profiles is very small, ts is determined by comparing the concentration

values of each compartment over a longer time period say t = 0.1. If the difference of every
corresponding compartment is less then 1 × 10−6, we say the transient solution has reached

a steady state.

∆x

0.1

0.05

0.02

0.0111111

0.00995803

tf

ts

0.1511

0.1367

0.1965

0.2059

0.2047

2.200

2.200

2.200

2.200

2.200

Table 4.2: Time to fully perfused (tf ) and to steady state (ts) corresponding to different
valid choices of spacial step size (∆x)

Table 4.2 shows the time to full perfusion (tf ) and to steady state (ts) corresponding to

different choices of spacial step size (∆x). As we can see that the full perfusion time is highly

influenced by the compartment size but converges to a certain value. The steady state time

is much more stable and not being influenced by the compartment’s size. This makes sense

because when there are too few compartments, the boundary compartments are thicker and

closer to the center, thus become nonzero earlier while not influencing the stability.

We are also interested about the numerical convergence of this algorithm under transient

107

state. Figure 4.3 shows effect of space step size at time t = 0.1 for two locations, the oxygen

concentration of the compartment at the center and at a random location. The vertical
axis represents the concentration value(×106) of the compartment and the horizontal axis

is the inverse of spacial step sizes (∆x) corresponding to different valid capillary shapes.

For each location, 24 different spacial step sizes are tested. It can be seen from both Fig.

4.3(a) and Table 4.3, that the concentration quickly converges to within a certain range and

oscillate inside that range as the size of ∆x decreases. This oscillation phenomenon may

due to the different cross sections of approximating the capillary compartment by squares.

For example, using 1, 4 or 9 compartments to approximate the capillary cross section the

shape is a square, but using 5 or 13 compartments, the capillary is cross-shaped. To further

analyze the convergence, we divide the result into sub-groups and calculate there averages.

∆x

0.1

0.05

0.02

0.0111111

t = 0.1

28437

40875

42190

41207

t = 0.2

59394

72609

70437

68811

t = 0.5

97947

108828

104477

102693

t = 1

107338

117156

113228

111468

t = 2

108028

117751

113936

112178

t = 4

108031

117754

113940

112182

t = 8

108031

117754

113940

112182

4.3:

Convergence

for

Table
random point
(0.243, 0.525, 0.725) corresponding to an interior tissue compartment for different choices
of valid spacial step size (∆x)

concentration value

on a

(×106)

108

Due to the convergence constraint and the given blood velocity, smaller ∆x requires

larger sliding frequencies corresponding to a smaller sliding period and thus requires extra

sliding steps during each time step. Vertical black lines in the figure divides the plane into

sub-regions corresponding to different sliding frequencies b = 1, 2, 3 and 4 from left to right.

Average of values in each sub-region is represented by the short black lines. Starting from

the smallest ∆x on the right end, convergence of the points is estimated using the uniformly

minimum variance unbiased estimator and represented by the long horizontal black line. We

found that, although individual points vibrate around the line without a clear pattern, the

average values of each sub-region gradually converge to the estimated value as expected.

109

(a) Concentration value(×106) on a random point (0.243, 0.525, 0.725)

(b) Concentration value(×106) on the central point (0.5, 0.5, 0.5)

Figure 4.3: Convergence of the 3D algorithm on interior tissue compartments for different
choices of valid spacial step size (∆x), at time t = 0.1. The vertical axis represents the
concentration value(×106) of the compartment and the horizontal axis is the inverse of
spacial step size (∆x). Vertical black lines divides the plane into sub-regions due to the
sliding period corresponding to the compartment size. Short horizontal lines are the average
of values in each rectangular region and the long horizontal line is calculated using the
uniformly minimum variance unbiased estimator from the end of the smallest ∆x obtained.

110

2040608010012030000350004000045000204060801001201100001200001300001400001500001600001700004.1.3 The Krogh Tissue Cylinder

In this section, we reexamine the clasical work of the oxygen transport from a single capillary

to its surrounding tissues by Krogh [38], taking account of the blood flow through the

capillary as well as the oxygen diffusion in the tissue in the axial direction. We shall compute

the steady state solution from unsteady state, and compare the results with the analytic

solution.

Figure 4.4: The Krogh tissue cylinder

Figure 4.4 shows the schematic diagram of a Krogh capillary-tissue cylinder model. In the

example we use a Krogh cylinder of 110 µm diameter tissue surrounding a 10 µm diameter

capillary of 1000 µm long. The tissue diffusivities are 900 µm2/s, which is assumesd to be

the same in the axial and radial direction. The blood flow inside the capillary is 400 µm/sec.

To determine the normalized concentrations, it is sufficient to only specify the ratio of the

tissue oxygen consumption to the oxygen concentration on the arterial side of the capillary.
Here we used κ(cid:48)/Ca = 0.002 sec−1.

111

Cc(z,t)vbCt(r,z,t)DκFigure 4.5: Concentration curves of the Krogh tissue cylinder at (C = 0.6, 0.7, 0.8, 0.9)
at steady state from the numerical solution. The vertical and horizontal axes shows the
compartment numbers used in each direction.

We aim to compute the steady state solution from the transient state. Assume there

are initially zero oxygen everywhere and at starting time blood with oxygen suddenly flows

in. Our numerical solution is computed under space step of ∆x = 0.0590818, which is

determined by the normalized capillary radius and using four compartments to represent
c = 4 · (∆x)2). A transient solution is calculated until reaching steady
one capillary (π · R2

state. Figure 4.5 shows the oxygen concentration curves for typical values obtained from

the numerical solution. Figure 4.6 shows almost perfect agreement of the numerical and the

analytic solution on the tissue boundary and on the capillary wall. Derivation of the analytic

solution can be find in the Appendix 6.2.

Figure 4.6: Concentration profile of the Krogh tissue cylinder on the capillary boundary and
the tissue boundary. Circles are numerical solution of oxygen concentration on the capillary
wall. Triangles are numerical solution of oxygen concentration on the tissue boundary.
Solid lines represents the oxygen concentration profile obtained from the analytic solution.
Horizontal axis shows the normalized capillary length and vertical axis show the normalized
concentration value.

112

0.60.70.80.905010015020025030051015202530○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○○△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△△����������������������4.2 Simulating 3D Multi-Capillary Oxygen Transport

using DCM

Parameter

c

R(cid:48)
L(cid:48)

FCD
v(cid:48)

b

Ca
D(cid:48)
κ(cid:48)
κ(cid:48)

b

max

Value
2 − 5 µm

150 - 1000 µm

1000 - 1500 capillaries/mm2

130 - 1500 µm/s

3 × 10−3 (ml O2)/(ml blood)
1.21 - 2.41 × 10−5 cm2/s

1.67 × 10−4 (ml O2)(ml tissue)−1 s−1
3 × 10−3 (ml O2)(ml tissue)−1 s−1

Reference

[37] [54]

[58] [23]

[67]

[24]

[50]

[24] [72]

[24]

[4]

Table 4.4: Biophysical parameters used for the simulation of oxygen transport in multi-
capillary tissue bed. R(cid:48)
density; v(cid:48)
diffusivity; κ(cid:48)

c, capillary radius; L(cid:48), capillary length; FCD, functional capillary
b, normal RBC velocity; Ca, arterial oxygen concentration; D(cid:48), tissue oxygen
b, basal tissue oxygen consumption; κ(cid:48)
max, maximal tissue oxygen consumption.

Consider a cube shaped multi-capillary bad with biological parameters are shown in table
4.4. As an illustrative example, let capillary radius R(cid:48)

c = 2.8µm, capillary length L(cid:48) = 240µm
and functional capillary density FCD = 1250 capillaries/mm2. Consider 8 capillaries and a

square cross sectional tissue domain. The dimensions of the tissue region is then calculated
from the functional capillary density, Ω(cid:48) = W(cid:48) × H(cid:48) × L(cid:48) = 80µm × 80µm × 240µm and
Ω = W ×L×H = 1×1×3 after normalization. Table 4.5 shows the capillary locations (from

reference [67]) and the plasma oxygen concentration on the arterial end for each capillary.

All the initial oxygen concentrations on the arterial end were set to be lower than that in the

113

artery and were randomly generated from the range of [0.8, 1]. Note that, the normalized

arterial oxygen concentration equals 1. Based on the convergence analysis result in section

4.1.2, we use 13 compartments to represent each capillary, which determines the space step

∆x = 0.0172057.

The steady state solution was first established from transient solution of zero oxygen ev-

erywhere and then blood suddenly flows into the capillaries. The tissue diffusivity 24.1µm/s

is assumesd to be the same in the axial and radial direction. Under the basal tissue oxygen

consumption rate and a blood velocity of 1150µm/s, oxygen concentration in the modeled

capillary-tissue system reaches steady state at normalized time t = 0.6.

m xm

ym

Ca of the mth capillary

1

2

3

4

5

6

7

8

0.385

0.879

0.183

0.673

0.837

0.855

0.602

0.690

0.138

0.365

0.691

0.428

0.310

0.161

0.760

0.203

0.960

0.935

0.978

0.928

0.931

0.995

0.976

0.962

Table 4.5: Locations and plasma oxygen concentration on the arterial end of eight capillaries.

Figure 4.7 shows the concentration profile of the capillary-tissue region at steady state.

Capillary-tissue region with oxygen concentration lower than 0.0001 is considered as anoxia.

Note that there is a tiny oxygen deficient region on the top corner of the top layer. Oxygen

114

diffused to that part was depleted during the same time step and so it is at the borderline

to anoxia.

We shall use this steady state result as the initial condition to simulate the development

of anoxia and its recovery in the following examples. We shall use the DCM to develop

tissue anoxia from two different ways: 1. by blocking capillaries and 2. by increasing the

tissue oxygen consumption rate. In the first case, the recovery process will be simulated by

unclogging the blocked capillary. In both cases, the recovery process will be simulated by

decreasing the consumption rate, by increasing capillary blood velocity and by increasing

oxygen concentration in the incoming plasma.

115

Figure 4.7: Oxygen concentration profile of eight capillaries of uneven strength and distri-
bution in a cuboid domain at steady state. Biophysical parameters collected in Table 4.4
is used to simulate the actual biological process. Under capillary blood velocity 1150µm/s
and basal tissue oxygen consumption rate 1.67 × 10−4 (ml O2)(ml tissue)−1 s−1, oxygen
concentration in this capillary-tissue system reaches the steady state at a normalized time
of t = 0.6.

116

4.2.1 Anoxia Following Capillary Occlusion

In the previous subsection, the steady state oxygen concentration in a capillary-tissue region

with 8 capillaries of uneven strength and experimental parameters are established. Now,

consider one capillary is suddenly blocked, which means blood in that capillary suddenly

stopped flowing. That means, new oxygen will stop entering the blocked capillary and only

the remaining oxygen in the blocked capillary is available for the surrounding tissue.

Figure 4.8 shows the change of oxygen concentration after blocking a capillary (m = 4)

which is closer to the middle of the tissue domain at time t = 0.01, 0.05, 0.1 and 0.2. At

t = 0.1, oxygen level around the blocked capillary decreased dramatically. However, the first

oxygen deficient region appears on the top margin (at t = 0.05), which is far away from the

blocked capillary and close to two healthy capillaries. This happens because oxygen always

diffuse from high concentration region to low concentration region, which enables capillaries

to support each other and balancing the oxygen distribution in the tissue region. This

surprising result is different from the traditional view, that tissue surrounding the blocked

capillary will be anoxic first. At t = 0.2 several oxygen deficient areas occur in separate

locations. The distribution of oxygen concentration reaches steady state again at t = 0.6.

Since the response time should be short in the reality, we shall use this result at t = 0.2 as

the initial condition on the next simulation.

The present result shows that, while occlusion happens, oxygen is redistributed to reach

a new balance with minimum anoxic area and this area may not occur near the blocked

capillary. Afterwards, anoxia may occur in several places.

117

t = 0.01

t = 0.05

t = 0.1

t = 0.2

Figure 4.8: The development of oxygen deficit by blocking one capillary closer to the middle
of the tissue domain, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic areas are shown in white
for contrast.

118

Figure 4.9 shows the results of oxygen distribution from blocking a capillary near the

outer tissue boundary (m = 1), after time t = 0.01, 0.05, 0.1 and 0.2. At t = 0.01,

oxygen concentration around the blocked capillary decreased greatly but unevenly due to

the distribution of other healthy capillaries. At t = 0.05, a significant anoxic area appears on

the veinous end. At t = 0.1, the zero oxygen region expands and extends to the arterial end.

Besides that, additional small anoxia regions starts to occur elsewhere near the veinous end.

The oxygen deficient regions keep expending from t = 0.1 to t = 0.2. Since the boundary

is defined to be reflected, this example demonstrates the result of blocking a bunch of four

neighboring capillaries.

This section is important for the study of heart attack and its recovery. A heart attack

or Myocardial infarction (MI) is the death of a segment of heart muscle usually caused by

a loss of blood supply. The blood is usually cut off when an artery supplying the heart

muscle is blocked by a blood clot. During a heart attack, the heart muscle loses oxygen

supply and is damaged. If impaired blood flow to the heart lasts long enough, the heart

cells die, chiefly through necrosis, and do not grow back. Treatment of an MI is time-

critical. Common treatments include nitroprusside, theophyline, coronary artery bypass

surgery, and thrombolysis, where the blood flow is restored; supplemental oxygen, where

oxygen concentration in the artery blood is increased and so on. In the following subsections,

DCM is applied to simulate different types of the recovery process.

For convenience, note that t = 0.01 is corresponding to 2.6556 seconds, t = 0.05 is

corresponding to 13.278 seconds and t = 0.2 is corresponding to 53.112 seconds.

119

t = 0.01

t = 0.05

t = 0.1

t = 0.2

Figure 4.9: The development of oxygen deficit by blocking one capillary near the tissue
domain boundary, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic areas are shown in white for
contrast.

120

4.2.1.1 Recovery Through Eliminating Capillary Blockage

In the previous simulation, one of the eight capillaries of uneven strength and distribution

was blocked, causing oxygen deficient in several separate regions (Figure 4.8). Now, consider

the blocked capillary was unblocked and resumes flow. We are interested in the transient

recovery process of tissue oxygen redistribution. Using the results from blocking the capillary

(m = 4) for a duration of t = 0.2 as the initial condition here, we restore the blood flow

of that capillary. Thus, original rate of oxygen flows into the recovered capillary again and

diffuse into the surrounding tissue region.

Figure 4.10 shows the oxygen concentration profiles after the blocked capillary is un-

clogged after time t = 0.001, 0.05, 0.1 and 0.2. At t = 0.001, the recovered capillary was

infused again with oxygenated blood and re-nourish the surrounding tissue. At t = 0.05,

oxygen deficient region in the middle has fully recovered but other anoxic regions remain.

From t = 0.05 to t = 0.1, the rest zero oxygen region shrunk greatly but still remain existed.

The concentration profile reaches a fully nourished steady state again at t = 0.5.

The result shows that, after eliminating the capillary blockage of a single capillary, oxygen

deficient regions near the blocked capillary will be recovered first. Anoxia tissue regions away

from the blocked capillary will need to wait a longer time for oxygen to diffuse there, and

thus have a greater chance to encounter tissue death.

121

t = 0.001

t = 0.05

t = 0.1

t = 0.2

Figure 4.10: Recovery process of tissue oxygen distribution through unclog the blocked
capillary, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown in white for contrast.

122

4.2.1.2 Recovery Through Decreasing Consumption Rate

In section 4.2.1, anoxia in capillary-tissue region with 8 capillaries of uneven strength and

experimental parameters was simulated by suddenly stopping the blood flow velocity inside

a capillary from steady state. Now, assume that capillary was still blocked and the living

system responds to anoxia by decreasing oxygen consumption. That means, the amount of

oxygen entering the system per time period remains the same as after the capillary occlusion,

but the less oxygen was consumed during each time period. So that extra oxygen will flow

towards anoxic regions. In reality, the basal tissue consumption rate can be further reduced

by decreasing the temperature, such as putting tissues into ice.

Consider the oxygen consumption rate is reduced from 1.67 × 10−4 (ml O2)(ml tissue)−1
s−1 to 1.55 × 10−4 (ml O2)(ml tissue)−1 s−1 at time t = 0.2 after the capillary blockage

happens. Figure 4.11 shows the recovery process of the tissue oxygen distribution after

reducing the oxygen consumption rate at time t = 0.001, 0.05, 0.1 and 0.2. At t = 0.001,

there is no significant change on the concentration profile. At t = 0.05, every anoxic region

is slightly reduced but none fully recovered. The zero oxygen regions continue to shrink from

t = 0.05 to t = 0.2.

The oxygen distribution reaches steady state t = 0.6 after the reduction of the consump-

tion rate. Under the selected new consumption rate, a significant oxygen deficient region

still exists at the middle of the top margin on the top layer. Tissues in this region are cer-

tainly going to die if nothing else happens, and tissues surround the blocked capillary on the

veinous end are in great danger. In addition, the top corner on the top level, which used to

have the lowest concentration in the steady state in the beginning (4.2), was no longer the

most anoxic.

123

t = 0.001

t = 0.05

t = 0.1

t = 0.2

Figure 4.11: Recovery process of tissue oxygen distribution following a reduction of the
tissue consumption, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown in white
for contrast.

124

4.2.1.3 Recovery Through Increasing Capillary Blood Velocity

In section 4.2.1, we developed anoxia in a capillary-tissue region with 8 capillaries of uneven

strength from steady state by blocking a single capillary and reduced the total oxygen supply.

Now, assume the living system is attempting to avoid tissue death by increasing the blood

velocity in the other healthy capillaries. That means, more oxygen will enter the region in

the same amount of time, resulting an increase on the amount to oxygen supply to tissue.

Physically, increasing of the heartbeat due to many reasons will lead to an increase of the

blood flow.

Assume the blood velocity of the seven healthy capillaries is suddenly increased from

1150µm/s to 1500µm/s, at time t = 0.2 after the capillary blockage happens. The resulting

oxygen concentration redistribution process at time t = 0.001, 0.05, 0.1 and 0.2 is shown

in Figure 4.12. At t = 0.001, there is no significant change on the concentration profiles.

However, at t = 0.05 most of the anoxic region has already recovered, leaving only tissues

near the blocked capillary and in the middle of two blocked capillaries (the boundary is

reflecting). Although the majority of the anoxia region recovered quickly, a zero oxygen

region exists at t = 0.2. Actually, the oxygen distribution reaches steady state at t = 0.5

with that region remain anoxic.

Increasing the blood velocity by 50µm every time, we found the anoxic region can be

fully recovered at vb = 1750µm. It takes t = 0.1922 for the entire tissue region to be perfused

and t = 0.5 to reach the steady state. The result shows us, by increasing the blood velocity,

the majority of the oxygen deficient region can be recovered in a short time. However, if

the velocity is not fast enough, some anoxic regions may exist permanently, leading to tissue

death.

125

t = 0.001

t = 0.05

t = 0.1

t = 0.2

Figure 4.12: Recovery process following an increment in the capillary blood velocity, at time
t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown in white for contrast.

126

4.2.1.4 Recovery Through Increasing Oxygen Concentration in Blood

In this simulation, we started with the oxygen shortage developed in section 4.2.1 by blocking

a single capillary. Assume the oxygen concentration in the incoming blood in the remaining

capillaries is suddenly increased by a certain amount, which means the total amount of

oxygen entering the capillary-tissue region will increase during each time period. Here, we

increase the arterial oxygen concentration Ca by 10%, at time t = 0.2 after the single capillary

occlusion occurs. Oxygen therapy, also known as supplemental oxygen, is a commonly used

medical treatment which increases the oxygen concentration in the artery blood.

Figure 4.13 shows the oxygen concentration profiles following increasing oxygen concen-

tration in blood, at time t = 0.001, 0.05, 0.1 and 0.2. The recovery process is very similar

to the result from increasing the blood velocity in the previous part 4.2.1.3. There is no

immediate change on the concentration profile at t = 0.001 but the anoxia region shrunk

greatly at t = 0.05. With the applied 10% increment on the blood oxygen concentration,

the oxygen distribution reaches steady state at t = 0.6 with a small zero oxygen area at the

top margin on the veinous end.

Note that, the anoxic regions could fully recovers with a greater increment on the oxygen

concentration in the artery blood. Here, we are aiming to explore and predict the location for

possible tissue death under severe cases. As the result shows, tissue death caused by capillary

occlusion may not happen on tissues neighboring the blocked capillary nor on tissues initially

having the lowest oxygen concentration but elsewhere.

127

t = 0.001

t = 0.05

t = 0.1

t = 0.2

Figure 4.13: Recovery process of tissue oxygen distribution following increasing oxygen con-
centration in the arterial blood, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are shown
in white for contrast.

128

4.2.2 Anoxia Following Increased Tissue Oxygen Consumption

In section 4.2, the steady state of oxygen concentration in capillary-tissue region with 8

capillaries of uneven strength and experimental parameters are calculated. Now, consider

the tissue oxygen consumption is suddenly increased. Increments on oxygen consumption

caused by muscle movement are usually accompanied by the increase of the blood velocity

or blood oxygen concentration. Here, we consider the case of acute exertion without proper

warm-up and focus on the occurrence of the resulting oxygen shortage and its recovery

processes.

Consider the tissue oxygen consumption rate is suddenly increased by 50%, from the
basal rate 1.67 × 10−4 (ml O2)(ml tissue)−1 s−1 to 2.505 × 10−4 (ml O2)(ml tissue)−1 s−1

after the steady state in section 4.2. Figure 4.14 shows the change of oxygen concentration

at time t = 0.01, 0.05, 0.1 and 0.2. Comparing to the result from blocking capillaries in

section 4.2.1, the oxygen shortage occurs more separately in the tissue domain at t = 0.01

and affects not only the veinous end at all levels. The oxygen deficient area expands to the

arterial end over time. The distribution of oxygen concentration reaches steady state again

at t = 0.4. Since the response time should be short in the reality, we shall use this result at

t = 0.2 as the initial condition of the following subsections.

129

t = 0.01

t = 0.05

t = 0.1

t = 0.2

Figure 4.14: The development of oxygen deficit by increasing the tissue oxygen consumption
rate, at time t = 0.01, 0.05, 0.1 and 0.2. Anoxic areas are shown in white for contrast.

130

4.2.2.1 Recovery Through Decreasing the Consumption Rate

In the previous subsection 4.2.2, anoxia in capillary-tissue region with 8 capillaries of un-

even strength and experimental parameters was developed by suddenly increasing the tissue

oxygen consumption rate. Now, assume the tissue oxygen consumption rate has decreased
back to the original 1.67 × 10−4 (ml O2)(ml tissue)−1 s−1 from 2.505 × 10−4 (ml O2)(ml
tissue)−1 s−1, at time t = 0.2. That means, although the amount of oxygen entering the

system remains the same, the oxygen concentration will be redistributed since every tissue

cell consumes less oxygen per time step allowing extra oxygen to flow back to anoxia regions.

Figure 4.15 shows the recovery process of the tissue oxygen distribution after reducing

the oxygen consumption rate at time t = 0.001, 0.01, 0.05 and 0.2. During the recovery

process, each of the oxygen deficient region gradually shrink and eventually disappear. The

oxygen distribution reaches steady state again t = 0.6, which is longer than the time taken

during the anoxia development process in the previous section 4.2.2.

Comparing the recovered area between different time periods, we observed that, even a

full recovery at the end is guaranteed, it took less time for the oxygen perfused region to

expand to tissues close to the capillaries than to tissues far from the capillaries.

131

t = 0.001

t = 0.05

t = 0.1

t = 0.2

Figure 4.15: Recovery process of tissue oxygen distribution following reducing the tissue
consumption rate to original state, at time t = 0.001, 0.05, 0.1 and 0.2. Anoxic areas are
shown in white for contrast.

132

4.2.2.2 Recovery Through Increasing Capillary Blood Velocity

In section 4.2.2, anoxia is developed by suddenly increasing the tissue oxygen consumption

rate at steady state. Muscle exercise always followed by an increase in the blood velocity.

From rest to exercise, blood flow velocity in the muscle capillaries can increase more then

tenfold [43]. Here, we consider two cases. Case 1: the capillary blood flow velocity is suddenly

increased to 1500µm/s t = 0.2 after the rise of tissue oxygen consumption rate, and case 2:

the capillary blood flow velocity is suddenly doubled. We are interested in how soon all the

anoxia regions are fully recovered.

Figure 4.16 shows the oxygen redistribution process at time t = 0.001, 0.01, 0.05 and 0.1

after suddenly increasing the capillary blood velocity to 1500µm/s. The oxygen deficient

regions were fully perfused with oxygen at t = 0.0569. Figure 4.17 shows the oxygen con-

centration profiles t = 0.001, 0.01, 0.05 and 0.1 after the capillary blood velocity is suddenly

doubled. At t = 0.0298, the anoxia regions have fully recovered. Comparing the two results,

we can see that there are significant zero oxygen regions remains at t = 0.01 in both cases,

However, a faster blood velocity lead to a much smaller fully-perfusion time, which means

faster further expansion of perfused region towards tissues located away from the capillaries.

133

t = 0.001

t = 0.01

t = 0.05

t = 0.1

Figure 4.16: Recovery process of tissue oxygen distribution following increasing the capillary
blood velocity to 1500µm/s, at time t = 0.001, 0.01, 0.05 and 0.1. Anoxic areas are shown
in white for contrast.

134

t = 0.001

t = 0.01

t = 0.05

t = 0.1

Figure 4.17: Recovery process of tissue oxygen distribution following doubling the capillary
blood velocity, at time t = 0.001, 0.01, 0.05 and 0.1. Anoxic areas are shown in white for
contrast.

135

4.2.2.3 Recovery Through Increasing Oxygen Concentration in Blood

In this simulation, we started with the oxygen shortage developed in section 4.2.2 by increas-

ing the tissue oxygen consumption rate. Assume the oxygen concentration in the incoming

blood has suddenly increased at time t = 0.2 after the consumption rate has increased. Here,

we increase the arterial oxygen concentration Ca by 10%, which means the total amount of

oxygen entering the capillary-tissue region is increased by 10% during each following time

period.

Figure 4.18 shows the oxygen concentration profiles following increasing oxygen concen-

tration in blood, at time t = 0.001, 0.01, 0.05 and 0.1. The recovery process is very similar

to the result from increasing the blood velocity in the previous part 4.2.2.2. There is no

immediate change on the concentration profile at t = 0.001 but the anoxia region shrunk

greatly at t = 0.01. With the applied 10% increment on the blood oxygen concentration,

the anoxic tissue region is fully recovered at t = 0.0668.

Note that both increasing the blood velocity and increasing the oxygen concentration

in the incoming blood increases the amount of oxygen entering the capillary-tissue region.

Hence, the resulting recovery process are similar, but increasing the arterial oxygen concen-

tration will lead to a higher oxygen concentration in tissues closer to the artery end.

136

t = 0.001

t = 0.01

t = 0.05

t = 0.1

Figure 4.18: Recovery process of tissue oxygen distribution following increasing oxygen con-
centration in the arterial blood, at time t = 0.001, 0.01, 0.05 and 0.1. Anoxic areas are
shown in white for contrast.

137

Chapter 5

The Analytic Theory of the Discrete

Compartmental Method

The proposed Discrete Compartment Method (DCM) is a numerical method that approx-

imates the solution of a transient diffusion-consumption-convection system by solving cor-

responding difference equations on a discrete domain with finite number of compartments.

In the examples from previous sections, we have validated the accuracy of steady state so-

lutions obtained using DCM by comparing them with analytic solutions of the differential

equation. However, for transient problems we generally do not have the exact solution to

compare against, nor do we have any highly accurate solution computed by other means or

experimental results (for example, transient problems in 2D and 3D; transient problems in

1D with anoxic region.) In such cases, we mush rely on some combination of the following

techniques to gain confidence in our numerical results:

ˆ Validation on test problems. The model with particular implementation should

be tested on simple problems for which the true solution is known, or on problems

for which a highly accurate comparison solution can be computed by other methods.

Experimental results may also be available for comparison in some cases.

ˆ Theoretical analysis of convergence and accuracy.

Ideally one would like to

prove that the numerical result converges to the correct solution as the compartment

138

is refined and obtain reasonable error estimates for the numerical error that will be

observed on any particular discrete domain.

In this chapter we concentrate on the theoretical analysis. Here we consider only the

diffusion consumption problem on a bounded cubic domain. We will also assume the problem

is well-posed, which means there exist a unique solution and it depends continuously upon

its initial value. In addition, we assume there is no moving boundary in the tissue region for

all t.

5.1 Numerical Error

In order to talk about accuracy or convergence, we first need to quantify the error and obtain

reasonable error estimates. Numerical solution using DCM is influenced by two sources of

error:

ˆ Discretization error, the difference between the exact analytical solution of the

partial differential equation and the exact solution of the compartmental difference

equation.

ˆ Round-off error, the difference between the exact solution of the compartmental

difference equation and the numerical result using finite-precision.

It is also called

rounding error, which is incurred after a repetitive number of calculations due to in-

exactness in the representation of real numbers on a computer.

In one dimensional problem we have an approximation Cn
i

for each compartment using

the discrete compartment method. For comparison, we let u(x, t) represent the analytical

solution of the partial differential equation and un

i represent the exact value we are hoping

139

to approximate well. For a finite difference method we may choose to compare with the

pointwise value

un
i = u(i∆x, n∆t),

(5.1)

but for a discrete compartment method we may want to instead use

(cid:90) xi+1/2

xi−1/2

un
i =

1
∆x

u(x, n∆t)dx.

(5.2)

However, the pointwise value (5.1) for the compartment centered at xi agrees with the
compartment average (5.2) to O(∆x) if the function u(x, t) is sufficiently smooth. Since

the present discrete compartment method is second-order accurate (will prove in the next

section), comparison with the pointwise value can be used for simplicity. The global error

at time T = N ∆t can then be represented as

EN = CN − uN .

(5.3)

Note that for time dependent problems, errors generally grow with time, and the number of

time steps to reach time level n + 1 will grow to infinity (in the limit that must be considered

in proving convergence) as we refine the grid. So we should not expect that any finite grid

could yield good solutions at arbitrarily large times.

We say that the method is convergent at time T in the norm || · || if

||EN|| → 0, ∆x, ∆t → 0,

(5.4)

and the method is pointwise convergent as the compartment is refined if the method is

140

convergent in the ∞−norm

||E||∞ = max

i∈Z∗ |Ei|.

We say that the method is accurate of order (˜p, ˜q) if

||EN|| = O(∆x˜p) + O(∆t˜q)

(5.5)

(5.6)

(cid:16)

∆x˜p + ∆t˜q(cid:17)

.

(5.7)

which means there exist a constant B, such that

||EN|| (cid:54) B

5.2 Convergence

Consider the Crank-Gupta problem



∂C

∂t = ∂2C

∂x2 − 1,

0 (cid:54) x (cid:54) s(t),

t (cid:62) 0

∂C
∂x = 0,

at x = 0,

C = ∂C

∂x = 0,

at x = s(t),

C(x, 0) = 1

2(1 − x)2,

0 (cid:54) x (cid:54) 1,

(5.8)

141

where s(0) = 1, and the discrete compartment method

(cid:3) − ∆t · κ, 0(cid:9) , n (cid:62) 0

i+1

i + Cn

(cid:3) − ∆t · κ, 0(cid:9) , n (cid:62) 0,
(cid:105) − ∆t · κ, 0
(cid:111)

, n (cid:62) 0,

(5.9)



Cn+1

i

Cn+1

0

Cn+1

I

= max(cid:8)Cn
= max(cid:8)Cn
(cid:110)

= max

Cn

i + α ·(cid:2)Cn
0 + 2α ·(cid:2)Cn
I + 2α ·(cid:104)

i−1 − 2Cn
1 − Cn
I−1 − Cn
Cn
i = 0, 1,··· , I,

0

i

i = 1
C0

2(1 − ∆x · i)2,

where α = ∆t/(∆x)2.

Theorem 5.2.1. If 0 < α (cid:54) 1/2, we have the solution of (5.9) pointwise converges to the

solution of (5.8), i.e.

||En||∞ → 0,

as ∆x, ∆t → 0,

∀n ∈ N.

(5.10)

Proof. First, we show that Cn
i

is decreasing w.r.t i in the region where its values are positive.

We have

i+1 − C0
C0

i =

1
2

[1 − (i + 1)∆x]2 − 1
2

(1 − i∆x)2

(5.11)

= (i +

)∆x2 − ∆x

1
2

Since ∆x = 1/I, for i = 1, ..., I − 1 we have

(i +

1
2

)∆x2 − ∆x =

i + 1
2

I

∆x − ∆x = −I − i − 1

2

I2

< 0.

(5.12)

142

Combining with (5.11) we have

i+1 − C0
C0

i < 0,

for i = 1, ..., I − 1.

(5.13)

Now consider the region where Cn+1

i+1 and Cn+1

i

> 0,

i+1 − Cn+1
Cn+1

i

=(1 − 2α)Cn

i + Cn

i+2) − ∆t · κ

i+1 + α(Cn
i − α(Cn
i+1 − Cn

− (1 − 2α)Cn
=(1 − 2α)(Cn

i−1 + Cn

i+1) − ∆t · κ
i+2 − Cn

i+1 + Cn

i ) + α(Cn

i − Cn

i−1)

Assume Cn

i+1 − Cn

i < 0, we have by (5.14) that

i+1 − Cn+1
Cn+1

i

< 0,

(5.14)

(5.15)

because α (cid:62) 0 and 1 − 2α (cid:62) 0. Recall by (5.13), this is also true for n = 0. Thus by

induction, we have

i+1 − Cn
Cn

i < 0

for all i = 1, ..., I − 1 such that Cn

i > 0, ∀n ∈ N.

On the left boundary, since Cn+1

0

< Cn

0 , by (5.9) we have

Cn+1

0

= (1 − 2α)Cn

0 + 2α · Cn

1 − ∆t · κ < Cn
0 .

(5.16)

(5.17)

143

Simplify and rewrite gives

0 > Cn
Cn

1 − ∆t · κ

2α

1 − ∆x2 · κ

= Cn

2

→ Cn

1

as ∆x → 0,

(5.18)

∀n ∈ N such that Cn

i > 0.

Thus, the solution Cn
i

(5.14) and (5.11) that

is decreasing strictly before it reaches 0. Moreover, we have by

Let(cid:101)in be such that Cn(cid:101)in

i

i+1 − Cn+1
|Cn+1
and Cn(cid:101)in+1
|Cn(cid:101)in

| (cid:54) |Cn

i+1 − Cn

i | → 0 as ∆x → 0.

= 0, by (5.19) we have ∀ 
 > 0, ∃ δ > 0 s.t.

| < 
 when ∆x < δ,

(5.19)

(5.20)

(5.21)

(5.22)

(5.23)

i.e. the numerical solution is asymptotically continuous as ∆x → 0.

Next, we show that the numerical solution converges to the original solution in its positive

values. Let

i = Cn
En

i − u(i∆x, n∆t),

then for i = 1, ..., I − 1 such that Cn

i > 0

En+1

i

− u(i∆x, (n + 1)∆t)

i

=Cn+1
=(1 − 2α)Cn
=(1 − 2α) [En

i+1) − ∆t · κ − u(i∆x, (n + 1)∆t)

i + α(Cn

i−1 + Cn

i + u(i∆x, n∆t)] + α(cid:2)En

i−1 + u((i − 1)∆x, n∆t) + En

i+1

u((i + 1)∆x, n∆t)] − ∆t · κ − u(i∆x, (n + 1)∆t),

144

∀n ∈ N. Consider the Taylor expansion of the analytical solution terms in equation (5.23)

at (i∆x, n∆t)

u((i − 1)∆x, n∆t) = u(i∆x, n∆t) − ∆x · ∂u
∂x
u((i + 1)∆x, n∆t) = u(i∆x, n∆t) + ∆x · ∂u
∂x
u(i∆x, (n + 1)∆t) = u(i∆x, n∆t)) + ∆t · ∂u
∂t

(cid:12)(cid:12)(cid:12)(i∆x,n∆t)
(cid:12)(cid:12)(cid:12)(i∆x,n∆t)
(cid:12)(cid:12)(cid:12)(i∆x,n∆t)

+ O(∆x2),

+ O(∆x2),

+ O(∆t2).

(5.24)

(5.25)

(5.26)

Put these back into (5.23), we have

En+1

i

=(1 − 2α)En

(cid:20)

+ α

i + α(cid:2)En

i−1 + En
u(i∆x, n∆t) − ∆x · ∂u
∂x

i+1

+u(i∆x, n∆t) + ∆x · ∂u
∂x

(cid:20)
i + α(cid:2)En

− ∆t · κ −

u(i∆x, n∆t)) + ∆t · ∂u
∂t

=(1 − 2α)En

i−1 + En

i+1

+ O(∆x2)

(cid:3) + (1 − 2α) · u(i∆x, n∆t)
(cid:12)(cid:12)(cid:12)(i∆x,n∆t)
(cid:12)(cid:12)(cid:12)(i∆x,n∆t)
(cid:12)(cid:12)(cid:12)(i∆x,n∆t)
(cid:3) + O(∆x2) + O(∆t).

+ O(∆x2)

+ O(∆t2)

(cid:21)

Since α (cid:62) 0 and 1 − 2α (cid:62) 0, by the triangle inequality,

|En+1

i

| (cid:54) (1 − 2α)|En

i | + α|En

i−1| + α|En

i+1| + B · (∆x2 + ∆t)

145

(cid:21)

(5.27)

(5.28)

for some constant B > 0. Take supremum over positive region, we have

sup

1(cid:54)i(cid:54)(cid:101)in+1

|En+1

i

1(cid:54)i(cid:54)(cid:101)in+1
| (cid:54) sup
1(cid:54)i(cid:54)(cid:101)i0
|En
(cid:54) sup

|En

i | + B · (∆x2 + ∆t)

i | + nB · (∆x2 + ∆t)

Observe that

We have

E0

i =

1
2

(1 − i∆x)2) − 1
2

(1 − i∆x)2) = 0.

(5.29)

(5.30)

(5.31)

sup

1(cid:54)i(cid:54)(cid:101)in+1

|En+1

i

| (cid:54) nA · (∆x2 + ∆t) → 0 as ∆x, ∆t → 0.

(5.32)

Therefore, for each (i, n), where i = 1, ..., I − 1 and n ∈ N, such that Cn

i > 0, we have

|En

i | = |Cn

i − u(i∆x, n∆t)| → 0 as ∆x, ∆t → 0.

On the right boundary, obviously we have

|En
I | = 0 ∀n ∈ N.

On the left boundary, we have

En+1

0

0 − u(0, (n + 1)∆t)

=Cn
=(1 − 2α)Cn
=(1 − 2α) [En

0 + 2αCn

1 − ∆t · κ − C(0, (n + 1)∆t)

0 + u(0, n∆t)] + 2α [En

1 + u(∆x, n∆t)]

− ∆t · κ − u(0, (n + 1)∆t)

146

(5.33)

(5.34)

(5.35)

(5.36)

(5.37)

(5.38)

∀n ∈ N. Consider the Taylor expansion of the analytical solution terms in equation (5.35)

at (i0, n∆t)

u(∆x, n∆t) = u(0, n∆t) + ∆x · ∂u
∂x

(cid:12)(cid:12)(cid:12)(0,n∆t)
(cid:12)(cid:12)(cid:12)(0,n∆t)

+ O(∆x2),

u(0, (n + 1)∆t) = u(0, n∆t)) + ∆t · ∂u
∂t

+ O(∆t2).

Put these back into (5.35), we have

En+1

0

=(1 − 2α)En

(cid:20)

0 + (1 − 2α)u(0, n∆t) + 2αEn
(cid:20)

(cid:12)(cid:12)(cid:12)(0,n∆t)

u(0, n∆t) + ∆x · ∂u
∂x

1

u(0, n∆t)) + ∆t · ∂u
∂t

(cid:12)(cid:12)(cid:12)(0,n∆t)

+ O(∆x2)

+ 2α

− ∆t · κ −

(cid:21)

+ O(∆t2)

(cid:21)

=(1 − 2α)En

0 + 2αEn

1 + O(∆x2) + O(∆t)

where (5.41) is because ∂u/∂x = 0 at x = 0. Since 0 < α (cid:54) 1/2

|En+1

0

0 | + |En

1 | + B1 · (∆x2 + ∆t)

| + |En

1 | + 2B1 · (∆x2 + ∆t)

| (cid:54) |En
(cid:54) |En−1
0
(cid:54) ... (cid:54)

1

| + |En−1
n(cid:88)
0| +
n(cid:88)

|Ek

k=0

(cid:54) |E0

= 0 +

1| + nB1 · (∆x2 + ∆t)

|Ek

1| + nB1 · (∆x2 + ∆t),

k=0

147

(5.39)

(5.40)

(5.41)

(5.42)

for some constant B1 > 0. Together with (5.33), we have

|En

0 | = |Cn

0 − u(0, n∆t)| → 0 as ∆x, ∆t → 0.

(5.43)

To handle the zero value part, recall by (5.20) that the numerical solution Cn
i

is asymp-

totically continuous. Thus ∀
 > 0, ∃δ s.t. (by (5.33))

− u((cid:101)in∆x, n∆t)| < 


|Cn(cid:101)in

if max(∆x, ∆t) < δ.

(5.44)

By (5.20), ∃δ1 > 0 s.t.

|Cn(cid:101)in

| < 


for ∆x < δ1.

(5.45)

Combining (5.44) and (5.45), we have

|u((cid:101)in∆x, n∆t)| = |u((cid:101)in∆x, n∆t) − Cn(cid:101)in
(cid:54) |u((cid:101)in∆x, n∆t) − Cn(cid:101)in

+ Cn(cid:101)in
|
| + |Cn(cid:101)in

|

< 2
.

(5.46)

Since u(x, t) is non-increasing and positive in x, we have

0 < u(((cid:101)in + 1)∆x, n∆t) < u((cid:101)in∆x, n∆t).

Recall by definition of(cid:101)in we have Cn(cid:101)in+1

= 0, i.e.

− u(((cid:101)in + 1)∆x, n∆t)| (cid:54) 2
.

|Cn(cid:101)in+1

148

(5.47)

(5.48)

Therefore, Cn
i

is also convergent to u(x, t) in the zero value region. Combining (5.33) and

(5.48), we have

i → u(i∆x, n∆t) as ∆x, ∆t → 0,
Cn

(5.49)

∀i = 1, ..., I and ∀n ∈ N, i.e.

||En||∞ = ||Cn − un||∞ → 0,

as ∆x, ∆t → 0,

∀n ∈ N.

(5.50)

Remark 5.2.1. The result can be easily extended to two-dimensional cases with 0 < α (cid:54) 1/4
and three-dimensional cases with 0 < α (cid:54) 1/6.

When the sliding element method (SEM) is introduced to deal with the convection,

the time step size is dependent on the blood velocity. Since our the original goal is to

simulate and analyze real biological process, biological parameters are used through out the

computation. For some combination of the parameters, we will need to adjust the sliding

frequency b accordingly to maintain the method’s stability. Based on the stability condition

of the discrete compartment method and the definition of the sliding period p, we have the

following corollary.

Corollary 5.2.1. When the sliding element method (SEM) is used along with the discrete

compartment method (DCM), for the convergence of an m-dimensional problem, we require

b (cid:62) 2m
∆x · vb

.

(5.51)

where vb is the normalized blood velocity and b ∈ Z+ is the sliding frequency.

149

Proof. For method convergence, we require

By the definition of ∆t (4.30), we have

∆t

(∆x)2

(cid:54) 1
2m

.

p
b

(cid:54) (∆x)2
2m

.

Finally, apply the definition of p gives

∆x
b · vb

(cid:54) (∆x)2
2m

.

(5.52)

(5.53)

(5.54)

Rewriting and simplifying the above equation yields the desired inequality (5.51).

The above corollary is crucial to the development of the discrete compartment model.

It helps determining the smallest b under selected blood velocity and space step size, which

satisfies the stability condition with the best accuracy.

Remark 5.2.2. In original coefficients with units and the way of normalization used in this

thesis, an m-dimensional DCM+SEM scheme is convergent, if

b (cid:62) 2m · D(cid:48)
∆x · ˜L · v(cid:48)

(5.55)

where ˜L =(cid:112)A(Ω(cid:48) ∩ {z = 0}) is the square root of the cross sectional area of the capillary-

b

tissue region Ω(cid:48) on the xy-plane in the unit of length, D(cid:48) is the diffusion coefficient in the
unit of area per time, v(cid:48)

b is blood velocity in the unit of length per time and b ∈ Z+.

150

Chapter 6

Conclusions

6.1 Discussion

The delivery of oxygen through capillaries and blood flow in the microcirculation is critical

in human physiology. Due to the advance of technology, it has been shown that many of the

macro-cycle system of major diseases, such as all kinds of stroke, ischemic heart disease, high

blood pressure, vascular occlusive diseases, diabetes, malignant tumor metastasis, etc, has

its occurrence and development not in the macro-circulation but in the micro-circulation.

A three-dimensional discrete compartmental model is developed to estimate the oxygen

concentration provided by multiple straight capillaries with different strengths among stri-

ated muscles. A reflecting boundary condition is assumed and blood flow velocity inside

the capillary is considered. Supply of oxygen from capillaries and tissue and the tissue oxy-

gen concentration level is critical in the function of all living creatures. Using the present

method, multiple capillaries with arbitrary characteristics in a three dimensional absorb-

ing tissue domain were first modeled in this paper. The proposed discrete compartmental

method is advantageous when compared to other methods due to the following four major

reasons:

1. The solution is transient. Previous works has obtained steady-state solutions in many

different setups, but the transient solution is more important has not been studied.

151

Actually, a true steady state may never exist in reality since the oxygen partial pres-

sure inside the capillary changes with breathing and cardiac pulse. Also, the tissue

consumption rate changes from rest to normal movement, to physical exercise.

2. The model is three-dimensional. Previous works only studied two-dimensional transient

problems with multiple capillaries, but these results cannot reveal the physiology of

the three-dimension case.

In two-dimensional models, the axial diffusion in tissue

was neglected. However, recent experiments on the tissue surrounding a capillary

clearly indicated that axial diffusion have significant effect on the oxygen concentration

profiles.

3. The model deals with anoxia among multiple capillaries. Some authors studied anoxia

using the Krogh tissue cylinder with a single capillary, but such case may never exist in

reality due to the concurrent flow. It is evident that capillaries have different potential

strengths and the distribution is uneven. Also the previous models do not predict the

location of anoxia which actually occurs somewhere away from a blocked capillary.

4. The model is constructed directly from Fick’s Law of diffusion. Therefore, the model is

easier to apply to three-dimensional case than previous models. At the same time, the

resulting difference equation is easy to compute using modern computer technology.

The model itself and its construction method can be applied to future studies of more

complicated biological systems.

The three major results are 1. The three dimensional transient oxygen concentration can

be solved in tissue supplied by a bunch of capillaries 2. If one capillary is blocked, tissue

anoxia dose not occur in the vicinity of the blocked capillary but occurs in some affected

152

region which is susceptible to low oxygen concentration 3. If neighboring capillaries are well

oxygenated, tissue anoxia may not occur at any time.

The present model is a significant contribution to the research field of both numerical

methods and micro-circulation. It reveals and simulates many important biological processes

which are not measurable or observable by experiments. Also, it demonstrates a successful

and novel way of developing efficient computational model for complicated biological phe-

nomena. Lastly and most importantly, it provides a new research tool with many possible

extensions on similar transport problems, especially for the study of heart attack and stroke

and their recovery.

6.2 Future Works

The three dimensional discrete compartmental modeling along with the time and space

matching technique is first developed in this paper. The discrete compartmental model is

designed especially for the oxygen transport problem in muscle tissue where the governing

difference equations are derived directly from the Ficks’ law of diffusion. Therefore, it is

much simpler than any other traditional numerical methods such as finite difference method

or finite element method and it can be used as a powerful tool to study the oxygen transport

among cellular units and oxygen concentration among muscle tissue.

Possible extensions of the discrete compartmental model are:

1. consider diverse type of cells with different consumption rates,

2. apply the model to other substrates such as carbon dioxide or substrates with finite

permeability,

3. instead of sudden blockage, consider different reduced blood velocities,

153

4. treatment of different boundary shapes,

5. use other compartment shapes such as hexagonal cube or tetrahedron,

6. instead of straight capillary, consider slightly curved capillaries,

7. instead of the zero oxygen area, study also the low oxygen areas,

8. simulation of large scaled tissue bed,

9. explore ways study low blood velocity without violate the convergence constraint,

10. study the change of individual and total outflow versus any abnormal oxygen distri-

bution such as anoxia,

11. consider capillaries of different radius,

12. consider oxygen axial diffusion inside the capillary,

13. deal with capillary branching,

14. calculate the anoxia time of individual zero oxygen regions and predict tissue death

time and location.

Theoretical proofs of consistency, convergence and stability of the discrete compartmental

method with moving boundary property remain to be improved. The model can be also be

applied to pharmacological problems such as the analysis of the drug delivery to tissues.

154

APPENDIX

155

APPENDIX

Analytic Solution of Krogh Tissue

Cylinder

A derivation of the steady state analytic solution of oxygen distribution in a Krogh cylinder

is given below. The model consists of a single capillary from which oxygen diffuses into a

surrounding cylinder of tissue. The equations governing this model were derived by Blum

[6] in 1960, who also attempted to obtain their solution but end up with incorrect results. A

correct mathematical analysis of a model for substrate concentration in tissue surrounding

single capillries is given by Salathe and Wang [57] in 1980. We shall apply their method to

the substrate of oxygen.

Frist, we will review and define the geometry and nature of the idealized cylindrical

capillary bad.

1. A single capillary of radius R(cid:48)

c and length L(cid:48) is surrounded by a concentric cylinder of

tissue of radius R(cid:48)

t. Both ends of the tissue cylinder is bounded by planes normal to the

capillary, and it is assumed that no oxygen crosses the boundary of the tissue cylinder

including the ends and the cylindrical sides. The origin of a cylindrical coordinate
system (x(cid:48), z(cid:48)) is place on the symmetry axis at the arterial end.

2. Blood containing oxygen molecules at a concentration C(cid:48)

a enters the capillary at the

arterial end and flows down to the venous end at a uniform velocity v(cid:48)

b. Inside the

156

capillary, we assume the hemoglobin was uniformly dissolved in the plasma and use

the average oxygen concentration at any cross section of the capillary. We also assume

that diffusion of oxygen inside the capillary are negligible compared to convection due

to blood flow through capillary.

3. As the oxygen is convected through the capillary at any longitudinal distance z(cid:48) from

the entrance, it diffuses through the membrane into the surrounding tissue. The per-

meability, P , of the capillary endothelium to the oxygen molecules is unlimited.

4. In the surrounding tissue the diffusivities of the oxygen in the radial and axial direction,

D(cid:48), is assumed to be the same, i.e. isotropic diffusion.

5. The surrounding tissue cells are assumed to be similar and has uniformly distributed

enzyme system that consumes the oxygen at a constant rate κ(cid:48).

6. Under the assumptions above, the concentration of oxygen in the capillary, C(cid:48)

c is a
founction only of the axial distance z(cid:48), while the concentration in the surrounding
tissue, C(cid:48)

t, depends on both z and the radial diatance r(cid:48).

Normalize all length by the tissue radius R(cid:48)
end C(cid:48)
z = z(cid:48)/R(cid:48)

t, l = L(cid:48)/R(cid:48)

c/R(cid:48)

a. In terms of the nondimensional variables Cc = C(cid:48)

t and oxygen concentration by that on the arterial
t/Ca, r = r(cid:48)/R(cid:48)
t,
t, the governing equation in the tissue arround the

a, Ct = C(cid:48)

c/C(cid:48)

capillary is

t and R = R(cid:48)
(cid:18)

∂2Ct
∂z2 +

1
r

∂
∂r

r

∂Ct
∂r

(cid:19)

= κ, 0 (cid:54) z (cid:54) l, R (cid:54) r (cid:54) 1,

(.1)

where κ = κ(cid:48)(R(cid:48)

t)2/D(cid:48)C(cid:48)
a.

Inside the capillary, the amount of oxygen passing through a cross section at a distance

z(cid:48) along the capillary can be denoted as πR(cid:48)2
c v(cid:48)

bC(cid:48)

c, and the amount leaving the capillary at

157

c v(cid:2)C(cid:48)

(cid:3). The excess, which is the amount of oxygen diffused

the distance z(cid:48) +dz(cid:48) is πR(cid:48)2
away to the surrounding tissues in a distance of dz(cid:48) is −πR(cid:48)2
c v(cid:48)

c + dC(cid:48)

c

b(dC(cid:48)

c/dz(cid:48))dz(cid:48). Meanwhile, the

amount of oxygen permeates across through the surface of the capillary membrane enclosed
by planes z(cid:48) and z(cid:48) + dz(cid:48) is given by the product of that area, 2πR(cid:48)
unit area, C(cid:48)

. These two flows have to be equal in the steady state which gives

cdz(cid:48), times the flow per

c − C(cid:48)

t|r(cid:48)=R(cid:48)

c

the expressions of the changes of the oxygen concentration inside and on the boundary of

the capillary cylinder.

(cid:104)

dC(cid:48)
dz(cid:48) = − 2P
cv(cid:48)
R(cid:48)

c

b

(cid:105)

, 0 (cid:54) z(cid:48) (cid:54) L(cid:48)

C(cid:48)
c − C(cid:48)

t|r(cid:48)=R(cid:48)

c

After the normalization we have

β

dCc
dz

= − [Cc − Ct|r=R] , 0 (cid:54) z (cid:54) l

cv(cid:48)
b/4R(cid:48)
where β = R(cid:48)

tP . The oxygen flux in the radial direction is

∂C(cid:48)
∂r(cid:48) |r=R = −P

t

D

and after normalization is

(cid:104)
C(cid:48)
c − C(cid:48)

(cid:105)

t|r(cid:48)=R(cid:48)

c

, 0 (cid:54) z(cid:48) (cid:54) L(cid:48),

∂Ct
∂r

|r=R = −η [Cc − Ct|r=R] , 0 (cid:54) z (cid:54) l,

(.2)

(.3)

(.4)

(.5)

where η = 2R(cid:48)

tP/D(cid:48).The inner and outer oxygen flows on the capillary membrane can then

be related by

γ

dCc
dz

=

∂Ct
∂r

|r=R, 0 (cid:54) z (cid:54) l

(.6)

158

cv(cid:48)
where γ = R(cid:48)
b/2D(cid:48). Since we assumed that permeability is unlimited, the boundary

condition on the capillary wall is given by taking the limit on equation (.3) with P going to

infinty

Cc(z) − Ct(R, z) = 0,

0 (cid:54) z (cid:54) l

(.7)

The other boundary conditions in terms of normalized veriables for this pair of differential

equations are

Cc = 1, z = 0

= 0, r = 1, 0 (cid:54) z (cid:54) l

∂Ct
∂r

= 0, z = 0, l, R (cid:54) r (cid:54) 1.

∂Ct
∂z

(.8)

(.9)

(.10)

We shall solve (.1) for Ct(r, z) in terms of Cc(z) by separation of variables, then substitute

it into (.6) to obtain the equation of Cc(z). Assume the corresponding homogeneous PDE

(cid:18)

(cid:19)

∂2Ct
∂z2 +

1
r

∂
∂r

r

∂Ct
∂r

= 0, 0 (cid:54) z (cid:54) l, R (cid:54) r (cid:54) 1,

(.11)

(.12)

(.13)

has a solution of the form

Substituting into (.11) yields

Ct(r, z) = F (r)G(z).

F (r)G(cid:48)(cid:48)(z) + F(cid:48)(cid:48)(r)G(z) +

F(cid:48)(r)G(z) = 0.

1
r

159

Separation of variables yields two ODE’s

F(cid:48)(cid:48)(r) +

1
r

F(cid:48)(r) − λ2F (r) = 0 and G(cid:48)(cid:48)(z) + λ2G(z) = 0

(.14)

with boundary conditions

F(cid:48)(1) = 0, and G(cid:48)(0) = G(cid:48)(l) = 0.

(.15)

Combine the solutions to (.14) that satisfies these boundary consitions express the solution

to (.11) by the eigenfunction expansion

∞(cid:88)

n=1

Ct(r, z) =

an [K1(λn)I0(λnr) + I1(λn)K0(λnr)] cos(λnz) +

κ
4

r2 − κ
2

ln r + A0

(.16)

where λn = nπ/l, for n = 1, 2, 3,··· . Substituting into the remaining boundary condition

(cid:16) nπz

(cid:17)

l

+

κ
2

angn cos

∞(cid:88)

n=1

Cc(z) =

(.7) yields

where

(cid:18) R2

2

(cid:19)

− ln R

+ A0,

(.17)

From the orthogonality of the eigenfunctions cos( nπz

that

and

gn = K1(λn)I0(λnR) + I1(λn)K0(λnR).

(.18)

(cid:17)

l ) on the interval 0 (cid:54) z (cid:54) l it follows
(cid:16)nπz
(cid:18) R2

(cid:19)

(.19)

dz,

l

(.20)

− ln R

.

2

(cid:90) l

0

an =

2
l gn

Cc(z) cos

(cid:90) l

0

1
l

A0 =

Cc(z)dz − κ
2

160

Substituting the solution of Ct(r, z) with constants given explicitly in Cc(z) into equation

(.6) and together with the boundary condition (.8) yields

(cid:16)nπz

(cid:17)(cid:90) l

l

0

2fn

γnπgn

sin

∞(cid:88)

n=1

(cid:16)nπω

(cid:17)

l

dω +

κ
2γ

(cid:19)

(cid:18)

R − 1
R

z + 1.

(.21)

Cc(ω) cos

Cc(z) =

where

fn = λn [K1(λn)I1(λnR) − I1(λn)K1(λnR)] .

Multiply equation (.21) by cos(cid:0) mπz

(cid:1) and integrate over dz from z = 0 to z = l yields a

(.22)

l

Fredholm integral equation of the second kind with separable kernel

cm = bm +

amncn, m = 1, 2, 3,··· ,

(.23)

∞(cid:88)

n=1

where

amn =

=

m = n or m + n even

(cid:16) mπz

(cid:17)

(cid:90) l

0

(cid:16)nπz

(cid:17)

l

2fn

dz

l

sin

cos

γnπgn

γπ2gn

4(l − 1)fn

0,
(cid:0)n2 − m2(cid:1) , m + n odd,
(cid:20) κ
(cid:18)
(cid:16) mπz
0,

R − 1
R

, m odd,

(cid:90) l

− κ l2
m2π2γ

R − 1
R

z + 1

cos

m even

2γ

0

(cid:18)

(cid:19)

(cid:19)

(cid:21)

l

bm =

=

(cid:17)

dz

(.24)

(.25)

161

and

(cid:90) l

0

cm =

(cid:16) mπz

(cid:17)

l

dz.

Cc(z) cos

(.26)

Equation (.23) represents an infinite set of linear algebric equations for the infinite set of

unknown constants cm, which can be truncated into the form

N(cid:88)

cm = bm +

amncn, m = 1, 2, 3,··· ,

Writting in matrix form

So that

n=1

Cm×1 = Bm×1 + Am×nCn×1.

(I − A) C = B.

This ”reduced” system can be solved as

C = (I − A)−1 B,

(.27)

(.28)

(.29)

(.30)

where is can be shown that cm converges to cm as N → ∞. By choosing a sufficiently large

N, one can obtain a solution of cm of any desired accuracy. Therefore, the solution to Cc(z)

and Ct(r, z) are found as

Ct(r, z) =

∞(cid:88)

n=1
κ
2

+

2cn
l gn

(cid:20)r2

2

[K1(λn)I0(λnr) + I1(λn)K0(λnr)] cos(

nπz

l

)

− ln r − R2
2

+ ln R

+

Cc(z)dz,

(.31)

(cid:21)

(cid:90) l

0

1
l

162

and

where

(cid:90) l

0

1
l

∞(cid:88)

n=1

2cnfn
γnπgn

sin

Cc(z) =

(cid:16)nπz

(cid:17)

l

+

κ
2γ

(cid:19)

(cid:18)

R − 1
R

z + 1,

(.32)

∞(cid:88)

n=1

Cc(z)dz =

2cnfn

γ(nπ)2gn

[1 − cos(nπ)] +

κ l
4γ

(cid:19)

(cid:18)

R − 1
R

+ 1.

(.33)

This result is used for the comparison with the numerical solution obtained using DCM

in example 4.1.3.

Note that, since the solution is obtained as an infinite series, it must be truncated after a

finite number of terms for numerical computations. An error analysis of the applied method

is given in [57]. In the example shown in this thesis, 100 terms were used in the expansion,

and the error was exceedingly small.

163

BIBLIOGRAPHY

164

BIBLIOGRAPHY

[1] Keumseong Bang and JaeGwi Go. The simulation of oxygen diffusion following occlu-

sion. Applied mathematics and computation, 183(2):1301–1309, 2006.

[2] James B Bassingthwaighte, Erik Butterworth, Bartholomew Jardine, and Gary M Ray-
mond. Compartmental modeling in the analysis of biological systems. In Computational
toxicology, pages 391–438. Springer, 2012.

[3] James B Bassingthwaighte, IS Joseph Chan, and CY Wang. Computationally efficient
algorithms for convection-permeation-diffusion models for blood-tissue exchange. An-
nals of biomedical engineering, 20(6):687–725, 1992.

[4] Daniel A Beard and James B Bassingthwaighte. Modeling advection and diffusion of
oxygen in complex vascular networks. Annals of biomedical engineering, 29(4):298–310,
2001.

[5] Alan E. Berger, Melvyn Ciment, and Joel C. W. Rogers. Numerical solution of a
diffusion consumption problem with a free boundary. SIAM Journal on Numerical
Analysis, 12(4):646–672, 1975.

[6] Jacob J Blum. Concentration profiles in and around capillaries. American Journal of

Physiology-Legacy Content, 198(5):991–998, 1960.

[7] Walter F Boron and Emile L Boulpaep. Medical physiology E-book. Elsevier Health

Sciences, 2016.

[8] Alessandra Bosutti, Stuart Egginton, Yoann Barnouin, Bergita Ganse, J¨orn Rittweger,
and Hans Degens. Local capillary supply in muscle is not determined by local oxidative
capacity. Journal of experimental biology, 218(21):3377–3380, 2015.

[9] John Crank. Free and moving boundary problems. Clarendon press Oxford, 1984.

[10] John Crank et al. The mathematics of diffusion. Oxford university press, 1979.

[11] John Crank and Radhey S Gupta. A method for solving moving boundary problems
in heat flow using cubic splines or polynomials. IMA Journal of Applied Mathematics,
10(3):296–304, 1972.

[12] John Crank and Radhey S Gupta. A moving boundary problem arising from the diffu-
sion of oxygen in absorbing tissue. IMA Journal of Applied Mathematics, 10(1):19–33,
1972.

165

[13] HO Dahmardah and DF Mayers. A fourier-series solution of the crank—gupta equation.

IMA Journal of Numerical Analysis, 3(1):81–86, 1983.

[14] Stuart Egginton and Eamonn Gaffney. Experimental physiology–review article: Tis-
sue capillary supply–it’s quality not quantity that counts! Experimental physiology,
95(10):971–979, 2010.

[15] H Fabel. Normal and critical o2 supply of the heart. In Oxygen transport in blood and

tissue, pages 159–171. George Thieme Stuttgart, 1968.

[16] John E Fletcher. A mathematical model of the unsteady transport of oxygen to tissues
in the microcirculation. In Oxygen Transport to Tissue, pages 819–825. Springer, 1973.

[17] John E Fletcher. A model describing the unsteady transport of substrate to tissue from

the microcirculation. SIAM Journal on Applied Mathematics, 29(3):449–480, 1975.

[18] John E Fletcher. Mathematical modeling of the microcirculation. Mathematical bio-

sciences, 38(3-4):159–202, 1978.

[19] Graham M Fraser, Daniel Goldman, and Christopher G Ellis. Comparison of generated
parallel capillary arrays to three-dimensional reconstructed capillary networks in model-
ing oxygen transport in discrete microvascular volumes. Microcirculation, 20(8):748–763,
2013.

[20] JaeGwi Go. Oxygen delivery through capillaries. Mathematical biosciences, 208(1):166–

176, 2007.

[21] D Goldman and Aleksander S Popel. Computational modeling of oxygen transport
from complex capillary networks. In Oxygen Transport to Tissue XXI, pages 555–563.
Springer, 1999.

[22] Daniel Goldman. Simulations of capillary network oxygen transport during transient
In Oxygen Transport To

ischemia in the presence and absence of tissue myoglobin.
Tissue XXIII, pages 355–359. Springer, 2003.

[23] Daniel Goldman. Theoretical models of microvascular oxygen transport to tissue. Mi-

crocirculation, 15(8):795–811, 2008.

[24] Daniel Goldman, Ryon M Bateman, and Christopher G Ellis. Effect of sepsis on skele-
tal muscle oxygen consumption and tissue oxygenation:
interpreting capillary oxygen
transport data using a mathematical model. American Journal of Physiology-Heart and
Circulatory Physiology, 287(6):H2535–H2544, 2004.

[25] Daniel Goldman, Ryon M Bateman, and Christopher G Ellis. Effect of decreased o2
supply on skeletal muscle oxygenation and o2 consumption during sepsis: role of hetero-

166

geneous capillary spacing and blood flow. American Journal of Physiology-Heart and
Circulatory Physiology, 290(6):H2277–H2285, 2006.

[26] Jos´e M Gonzalez-Fernandez and Susie E Atta. Transport and consumption of oxygen

in capillary-tissue structures. Mathematical Biosciences, 2(3-4):225–262, 1968.

[27] Carl A Goresky and Harry L Goldsmith. Capillary-tissue exchange kinetics: diffusional
interactions between adjacent capillaries. In Oxygen Transport to Tissue, pages 773–781.
Springer, 1973.

[28] Radhey S Gupta and Dhirendra Kumar. Complete numerical solution of the oxygen dif-
fusion problem involving a moving boundary. Computer Methods in Applied Mechanics
and Engineering, 29(2):233–239, 1981.

[29] E Hansen and P Hougaard. On a moving boundary problem from biomechanics. J.

Inst. Math. Appl., 13, 06 1974.

[30] Archibald Vivian Hill. The diffusion of oxygen and lactic acid through tissues. Pro-
ceedings of the Royal Society of London. Series B, Containing Papers of a Biological
Character, 104(728):39–96, 1928.

[31] G Hutten, G Thews, and P Vaupel. Some special problems concerning the oxygen supply
to tissue, as studied by an analogue computer. Oxygen Supply, Urban & Schwarzenberg,
Munchen, page 25, 1973.

[32] John A Jacquez. Compartmental analysis in biology and medicine, biomedware. Ann

Arbor, MI, 512, 1996.

[33] John Alfred Jacquez et al. Compartmental analysis in biology and medicine. 1972.

[34] John A Johnson and Theodore A Wilson. A model for capillary exchange. American

Journal of Physiology-Legacy Content, 210(6):1299–1303, 1966.

[35] Leonard R Johnson. Essential medical physiology. Elsevier, 2003.

[36] Seymour S Kety. Determinants of tissue oxygen tension.

In Federation proceedings,

volume 16, page 666, 1957.

[37] Ronald J Korthuis. Skeletal muscle circulation.

In Colloquium Series on Integrated
Systems Physiology: From Molecule to Function, volume 3, pages 1–144. Morgan &
Claypool Life Sciences, 2011.

[38] August Krogh. The number and distribution of capillaries in muscles with calculations of
the oxygen pressure head necessary for supplying the tissue. The Journal of physiology,
52(6):409–415, 1919.

167

[39] August Krogh. The anatomy and physiology of capillaries, volume 18. Yale University

Press, 1922.

[40] Randall J LeVeque et al. Finite volume methods for hyperbolic problems, volume 31.

Cambridge university press, 2002.

[41] David G Levitt. Theoretical model of capillary exchange incorporating interactions
between capillaries. American Journal of Physiology-Legacy Content, 220(1):250–255,
1971.

[42] Edwin N Lightfoot et al. Transport phenomena and living systems. 1974.

[43] Ted Lin and Chris Mowatt. Physiology of the circulation, page 315–343. Cambridge

University Press, 4 edition, 2016.

[44] Stanley Middleman. Transport phenomena in the cardiovascular system. John Wiley &

Sons, 1972.

[45] JV Miller, KW Morton, and MJ Baines. A finite element moving boundary computation

with an adaptive mesh. IMA Journal of Applied Mathematics, 22(4):467–477, 1978.

[46] S.L. Mitchell. An accurate application of the integral method applied to the diffusion of
oxygen in absorbing tissue. Applied Mathematical Modelling, 38(17):4396 – 4408, 2014.

[47] SL Mitchell and Michael Vynnycky. The oxygen diffusion problem: Analysis and nu-

merical solution. Applied Mathematical Modelling, 39(9):2763–2776, 2015.

[48] David L Nelson, Albert L Lehninger, and Michael M Cox. Lehninger principles of

biochemistry. Macmillan, 2008.

[49] WHO. “The top 10 causes of death.” World Health Organization, 2018. http://www.

who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death

[50] Roland N Pittman. Regulation of tissue oxygenation. In Colloquium series on inte-
grated systems physiology: from molecule to function, volume 3, pages 1–100. Morgan
& Claypool Life Sciences, 2011.

[51] MJ Plyley, Gloria J Sutherland, and AC Groom. Geometry of the capillary network in

skeletal muscle. Microvascular research, 11(2):161–173, 1976.

[52] Aleksander S Popel. Theory of oxygen transport to tissue. Critical reviews in biomedical

engineering, 17(3):257, 1989.

[53] Aleksander S Popel and Paul C Johnson. Microcirculation and hemorheology. Annu.

Rev. Fluid Mech., 37:43–69, 2005.

168

[54] RF Potter and AC Groom. Capillary diameter and geometry in cardiac and skeletal
muscle studied by means of corrosion casts. Microvascular research, 25(1):68–84, 1983.

[55] Daniel D Reneau, Duane F Bruley, and Melvin H Knisely. A mathematical simulation
of oxygen release, diffusion, and consumption in the capillaries and tissue of the human
brain. In Chemical Engineering in Medicine and Biology, pages 135–241. Springer, 1967.

[56] Francis John Worsley Roughton. Diffusion and chemical reaction velocity in cylindrical
and spherical systems of physiological interest. Proceedings of the Royal Society of
London. Series B-Biological Sciences, 140(899):203–229, 1952.

[57] Eric P Salath´e and Tseng-Chan Wang. Substrate concentrations in tissue surrounding

single capillaries. Mathematical Biosciences, 49(3-4):235–247, 1980.

[58] Eric P Salath´e, Tseng-Chan Wang, and Joseph F Gross. Mathematical analysis of

oxygen transport to tissue. Mathematical Biosciences, 51(1-2):89–115, 1980.

[59] Eric P Salath´e and Tseng-Chang Wang. The development of anoxia following occlusion.

Bulletin of mathematical biology, 44(6):851–877, 1982.

[60] Lauralee Sherwood. Human physiology: from cells to systems. Cengage learning, 2015.

[61] IA Silver. Some observations on the cerebral cortex with an ultramicro, membrane-
covered, oxygen electrode. Medical electronics and biological engineering, 3(4):377–387,
1965.

[62] Liang Sun. Substrate diffusion and consumption in rectangular capillary-tissue bed.

Letters in Biomathematics, 1(1):67–76, 2014.

[63] Liang Sun, Junkoo Park, and Alessandra L. Barrera. A compartmental model for

capillary supply. Letters in Biomathematics, 4(1):133–147, 2017.

[64] Gerhard Thews. ¨Uber die mathematische behandlung physiologischer diffusionsprozesse

in zylinderf¨ormigen objekten. Acta Biotheoretica, 10(3-4):105–138, 1953.

[65] Eleuterio F Toro. Riemann solvers and numerical methods for fluid dynamics: a prac-

tical introduction. Springer Science & Business Media, 2013.

[66] Gerard J Tortora and Bryan Derrickson. Principles of anatomy & physiology. John

Wiley & Sons, Incorporated, 2017.

[67] Karel Tyml, Odile Mathieu-Costello, Linong Cheng, and Earl G Noble. Differential
microvascular response to disuse in rat hindlimb skeletal muscles. Journal of Applied
Physiology, 87(4):1496–1505, 1999.

169

[68] CY Wang and JB Bassingthwaighte. Capillary supply regions. Mathematical biosciences,

173(2):103–114, 2001.

[69] Harvey R Weiss. Control of myocardial oxygenation—effect of atrial pacing. Microvas-

cular research, 8(3):362–376, 1974.

[70] WJ Whalen, D Buerk, and CA Thuning. Blood flow-limited oxygen consumption in
resting cat skeletal muscle. American Journal of Physiology-Legacy Content, 224(4):763–
768, 1973.

[71] Jonathan P Whiteley, David J Gavaghan, and Clive EW Hahn. Mathematical modelling

of oxygen transport to tissue. Journal of mathematical biology, 44(6):503–522, 2002.

[72] Jonathan P Whiteley, DJ Gavaghan, and Clive EW Hahn. Oxygen transport to muscle
tissue where regions of low oxygen tension exist. Mathematical and computer modelling,
42(9-10), 2005.

[73] Mohamed Zerroukat and Chris R Chatwin. Computational moving boundary problems,

volume 8. Research Studies Press, 1994.

170