A compartmental model for capillary supply

Abstract Oxygen diffusion for time-dependent diffusion and consumption can be measured for small tissue regions containing a single capillary. An all or none model is reflected by myocardial infarction where necrotic regions are clearly demarcated. However if there is more than one capillary, the problem becomes very difficult; since the boundary of the ischemic area is no longer circular and is not known a priori. A geometric compartmental model using the Fick’s method will be presented for multi-capillary supply. Our method is to approach the steady state by a transient process, which paradoxically may be more efficient than the steady-state problem.


Introduction
Capillaries are tiny vessels connecting arterioles with venules and forming networks through the body. These thin-walled blood vessels are pervasive throughout the organism forming networks called capillary beds. Access to the gases and other nutrients provided by the capillary ensures that cells typically are in close proximity. Oxygen that is bound to erythrocytes releases and diffuses into the nearby tissues in a passive process. The high oxygen content in the blood promotes the oxygen gas to diffuse into the low oxygen tissues.
Oxygen diffusion from the microcirculation through capillaries, and its consumption by tissue cells, is basic in all higher organisms. Krogh (1919) introduced the first basic mathematical model to estimate the oxygen concentration profiles for highly regular capillary bed of skeletal muscle. Extensions of the Krogh cylinder model have been made, such as polygonal regions of supply, capillary effects, axial dependence and oxygen pressure.
Krogh cylinder along with many other models approximated oxygen distribution generated from homogeneous capillaries within regular capillary boundaries. Impermeable barriers and zero interactions among capillaries were also assumed, but in real cases the capillaries may not be parallel and may be arranged more or less randomly in normal tissues. The oxygen concentration, along with its time-wise distribution, may be uneven. Therefore, heterogeneity is considered much more common (Degens, Deveci, Bemden, Hoofd, & Egginton, 2006). An important case of heterogeneity is abnormal perfusion. This may be caused by local ischemia due to blockage such as in heart attack, ischemic stroke or local collapse due to pressure and trauma. Obstruction of these capillaries by a clot or by injury prevents the downstream cells from access to vital resources such as oxygen. Oxygen starvation in tissues is life threatening in organs with high-energy demands such as the heart and brain. Tissue necrosis is irreversible and subsequently results in permanent damage to the organ and potential death of the organism. The impact to these tissues is dependent upon the extent of the area affected by the obstructed capillary. Where does such hypoxia occur? Would neighbouring capillaries provide enough oxygen to the anoxic regions? How soon and how much? These are the questions to be answered. Note that we are concerned about the scale of distance between capillaries and not global scale problems such as increase in demand due to exercise or reactive hyperaemia after a tourquinet is released.
Consider a cross-sectional region perfused with N parallel capillaries at arbitrary locations with arbitrary diffusion fluxes q i (1 ≤ i ≤ N). The axial variation is small and can be neglected. The unsteady diffusion-consumption equation suffices where C is the normalized with respect to oxygen concentration of the arterial blood. The oxygen consumption κ(x, y, t) is the consumption rate at location (x, y) and time t, and at the steady state is equal to the sum of the fluxes out of all the capillaries. Due to uneven capillary characteristics, the capillary domains could be highly contorted and differ from either the Krogh circular cylinder (Krogh, 1919) or Voronoi polygonal cylinders (Gonzalez-Fernandez & Atta, 1972). The boundary of the perfused region p(x, y, t) = 0 is unknown. On this boundary where n is the normal to the boundary. The problem is a difficult moving boundary problem or known as a Stefan problem (Crank, 1975), studied in the context of freezing and melting of materials and seepage flow in porous media. There exist some analytical solutions to Equations (1) and (2) in one dimension including axisymmetric. For two-dimensional problems, such as the case for multiple capillaries, even direct numerical interpolation fails due to the unknown moving boundary. To locate the boundary for numerical computation, several methods are suggested. These include fixing a partial grid on the moving boundary, a flexible grid that include the boundary, grid lines that coincide with time step and exchange of dependent and independent variables. All are impossible to implement in two dimensions. The question is similar to the blockage of one capillary in Figure 1 where we have an ischemic area of unknown boundary. This area may be decreased by an increase in flux from the neighbouring perfused capillaries.
The 2D ischemic quasi-steady problem was attempted by Salathé and Wang (1982) using a Krogh cylinder and a perturbation from steady. The best method to treat truly 2D unsteady problem with moving boundary is the truncation scheme proposed by Berger, Ciment, and Rogers (1975). They used a large, fixed domain, and set or truncate any negative values of concentration to zero. Since the domain is fixed, the usual finite difference numerical schemes can be used. The moving front is then the zero concentration curve. The method is applied to a 2D oxygen diffusion problem by Evans and Gourlay (1977) who used a hopscotch finite difference method to increase stability.
In what follows we shall propose a compartmental diffusion consumption model, which is much simpler than the truncation-finite difference method of Berger et al. (1975).

Essence of the compartmental model
Consider the whole plane partitioned into contiguous finite tissue bins. Each bin would consume one unit of oxygen per unit of time, if sufficient oxygen supply is available. If oxygen is less than the amount to fully suffice the plane, the bin is hypoxic and therefore consumes all available oxygen. Some of the bins are designated as capillaries, which produce relatively large amount of oxygen per unit of time. By Fick's diffusion law through membrane separations, those bins with high substrate/oxygen concentration would diffuse into adjacent bins where substrate/oxygen concentration is low. Within each bin, the oxygen content is even or well mixed. Therefore, the whole system is a network of physical compartments. The resulting equations are similar to the discretized diffusionconsumption equations but are formulated from completely different principles.
We will use 1D first to derive the governing equations for such compartment model, as shown next.

1D unsteady diffusion
Consider the one dimensional problem C(x, t), x ≥ 0, where the concentration at x = 0 when t = 0 is suddenly raised and kept at a normalized value of 1. The governing equation is with the boundary and initial conditions We shall compare it with the exact solution where erfc is the complementary error function (Crank, 1975). Instead of the partial differential equations approach, consider a domain x ≥ 0, which is composed of ordered compartments of width x in one dimension. By Fick's diffusion law, each compartment is assumed to exchange oxygen with adjacent compartments and the rate of change of oxygen moving is proportional to transfer coefficient, surface area and difference of concentration between adjacent compartments. Then from the mass balance, the following equation can be obtained where C i is the concentration of the ith compartment. Thus where is the area of diffusion and k is the transfer coefficient. Equation (7) can be rewritten as In the limit of x → 0, we see the last bracket is ∂ 2 C/∂x 2 ; thus, (k x) must be the diffusivity, which we have normalized to unity in Equation (3). Now consider very small time step t and replace dC i /dt by (C i,n+1 − C i,n )/ t where n denotes the time. Equation (8) now becomes the algebraic equation where α = t/( x) 2 , which is in the order of o(1). This is exactly the finite difference explicit discretization of Equation (3).
In the beginning, all compartments are empty except the concentration of the zeroth compartment is maintained at 1.

1D unsteady diffusion-consumption
In a one-dimensional unsteady diffusion-consumption problem, a normalized constant consumption rate of 1 is considered for simplicity. Let C(x, t) be the oxygen concentration at time t and at location x and let be the perfused region of tissue. Then the governing equation is obtained by with the boundary and initial conditions The solution for the final steady state is As shown, the region 0 ≤ x ≤ √ 2 is fully perfused and ischemic elsewhere. We shall use the compartmental method on the transient problem. Balancing mass on the ith compartment, we obtain where is the area of diffusion and k is the transfer coefficient. Thus If C i,n+1 is found to be negative, it is set to be zero and the compartment is labelled ischemic. Figure 3 shows the spread of the oxygen diffusion from a plane into a consuming region. Note when t = 1 the computed results have reached the steady-state exact solution in Equation (16).
For stability, similar to finite difference explicit methods, α shall be less than .5. Otherwise implicit or semi-implicit schemes would improve stability, but introduce systems of equations to be solved.

2D compartmental model
For one-dimensional unsteady diffusion, the compartmental model leads to similar equations as the finite difference method. For 2D, the compartmental model has decisive advantages.
Here hexagonal compartments are utilized to cover the whole tissue region. Firstly, the compartmental model is quite simple to implement. In comparison, the 2D unsteady finite difference schemes are quite complicated. Secondly, finite differences are normally used in a square grid, while the compartments used here can be of any shape, i.e. one can use tissue cell boundaries or Voronoi domains (Gonzalez-Fernandez & Atta, 1972). We shall use hexagonal compartments shown in Figure 4(b), which are more natural than squares.
Notice in Figure 4(a) the spread of a square grid, having only four neighbours, is more distorted or uneven than the spread of a hexagonal grid, which has six neighbours. The neighbours of the (i, j)th cell for hexagonal compartments from the grid lines in Figure 4(c) are (i − 1, j + 1), (i, j + 1), (i − 1, j), (i + 1, j), (i, j − 1), (i + 1, j − 1).
Mass balance gives Here C ij is the concentration of the (i, j)th cell, κ is the transfer coefficient, or permeability through the cell boundary and κ is the consumption as in mass per volume per time.
Comparing the above with the diffusion equation we find diffusivity Normalizing all lengths by R, which is a measure of the size of the region, the concentration by κR 2 /D and the time by R 2 /D, Equation (22) becomes To discretize the time using explicit method, if time is discretized by t = n · t C i,j,n+1 = C i,j,n + α neighbours − 6C i,j,n − t; here α = 2 t/3( x) 2 and The coordinates (x, y) are related to (i, j) by For this free boundary problem, one can choose a large region that contains the area where oxygen is well sufficed and consumed. For example, if the outer boundary is not oxygenated, we can use a rhombic region where i = −N, N, j = −N, N large enough to contain any perfused regions. On the rhombic boundary, one can prescribe C = 0, which is equivalent to the zero flux condition ∂C/∂n = 0. In the case the zero flux condition is truly needed, the concentrations of adjacent cells in the normal direction can be prescribed equal. For example alone j = −N, or even better, remove the two outside neighbouring bins or remove the two neighbouring bins C −N,j , C −N,j+1 . Similar relations can be obtained for Dirichlet boundary conditions y = constant or any curved boundary. Now consider one of the cells is a capillary or oxygen source. Suppose the flux is Q, nourishing a region of πr 2 in the steady state. This region contains πr 2 / − 1 consuming cells. If each cell consumes one unit of substrate per time, the total flux needed from the capillary is Let r = J x, then from Equation (21) For a radius of r = 1, x = .1, Q = 271, a radius of .5, Q = 67. Now by Fick's diffusion law, the oxygen flux is evenly distributed through the membrane wall between each two adjacent cells, therefore Q/6 amount of oxygen is supplied to each of the capillary neighbouring cells in time of t. The governing equation of oxygen concentration for capillary neighbouring compartments (i, j) is

Examples and discussion
I. To illustrate this model, we start first with the case of a single capillary in an infinite domain. For an unbounded domain with one single capillary, it does not supply the needs of the entire tissue area. The region is ischemic outside the circular supply boundary where r ≤ R. Let C denote the oxygen concentration by perfusion from a single capillary with oxygen strength q at the centre. The governing equation is ) and the oxygen flux is zero at r = R: The solution to Equations (34)-(36) is where the constants R and p are found by Equation (35): The oxygen concentration decreases from the centre of the source through the tissue area at √ q. This region enclosed by the circular boundary is also referred as the capillary supply region (Wang & Bassingthwaighte, 2001). Figure 5 shows the result.  First, we look at the case with one single capillary in the centre and flux strength sufficiently enough for a circular supply region with radius r = 1.
As t → ∞, the concentration profile should approach  which is from the steady-state solution for Equation (37) when q = 1. Here we use α = .1, x = .01 and t = .5. The result is displayed as in Figure 6, which shows the oxygen concentration profiles at t = 1.
When time is small, the spread is approximately hexagonal, showing the effect of the grid. For larger time, the spread becomes circular. Table (1) shows the radius of the supply region against time. At t = .5 the radius reaches .972, which is close to the steady-state estimate of 1. Figure 7 shows the radial oxygen distribution error as results are compared to Equation (39). As the compartment partition becomes smaller, we can adjust the size of our capillary source by adding a fixed number of compartments so that the capillary wall is approximated by the outer region of the source compartment boundary.
II. Using the compartmental model, we can construct boundaries with different shapes, as well as for studying insulated closed regions. Due to the natural properties of hexagonal boundaries, the insulated region boundary can be hexagonal or rhombic. We can also    (26) and make flux zero on each boundary cell walls.
Note that for Neumann boundary conditions such as the fully perfused insulated region, the steady concentrations are unique up to an additive constant, while for Dirichlet or mixed boundary conditions such as partially ischemic cases, the steady concentrations are unique. This means the concentrations are independent of history, or of how steady state is achieved.
For multiple capillaries with insulated regions, we approximate the rectangular region by the hexagonal compartmental model. The following sources are given in Table (2). Use x = .004, α = .1. The total flux strength from three capillaries is calculated to satisfy the metabolic needs of the whole rectangular region. Figures 9 and 10 show the oxygen concentration at t = .5 and 3.5, respectively. The difference between the profiles at t = 3.5 and 4 is very small (< 10 −4 ) and we assume a steady state is well achieved. Figure 11 shows the oxygen concentration comparison between the compartmental model and the analytical results using the series solution (Sun, 2014). Minimum oxygen concentration can be obtained on the boundary using the exact solution (Sun, 2014). The two results differ by a constant, and this constant can be eliminated from consideration due to Neumann type of boundary conditions. One can, from the concentration distribution, determine the capillary domains of each source. The method is to find local minimums along the three grid lines i = constant, j = constant, and i + j= constant. The boundary of the capillary domains must contain three local minimums. However, due to the approximate nature of the course grid, it may be difficult to discern a true minimum especially at the outlying areas where the concentrations are almost uniform. Other geometric types of compartmental shapes can also be applied, with the formulation similar to the one proposed in the previous section.
For 2D supply regions, the finite difference method (Berger et al., 1975) used relies heavily on the grid network chosen. For oxygen diffusion problems where multiple capillary sources reside, these capillaries create singular points. In Berger's and almost all the other methods, due to lack of mesh points and high concentration, the solutions near the singular points either require smaller time steps or special formulation. Gupta and Kumar (1981) used a variable time step grid network to avoid the large number of time steps generally required. The compartmental model we propose only requires uniform time steps and uniform spatial discretization, and therefore it is much easier to implement. In the traditional numerical methods for Stefan problems, it is always necessary to evaluate the time to reach the oxygen front. Our method does not require such extrapolation.
The Compartmental Model, which is also a diffusion-consumption model, points out a direction when the boundary is constantly moving (i.e. a Stefan problem). Using Fick's law of diffusion, the hexagonal compartmental model can be used to describe two-dimensional unsteady diffusion-consumption for multi-capillary supply regions. The method models transport among cellular units and is much simpler than traditional finite differences.

Disclosure statement
No potential conflict of interest was reported by the authors.