The tissue that we are simulating is 100 × 100 × 100 μm3
in size and contains a single blood vessel with diameter of d
= 10 μm. The vessel is passing through the center of the tissue and parallel to the x
-axis. Unless otherwise stated, for this simple simulation the velocity of blood inside the vessel is 400 μm/s, hematocrit is set to 0.1, and PO2
= 90 mmHg at the input boundary of the vessel. The vessel wall permeability is Kw
μmol/(μm s mmHg) and wall thickness is 1 μm. The oxygen diffusion coefficient is DO2
/s. Oxygen consumption is 41.7×10-16
s), or 10% of the baseline value for the rat cortex [31
]. A mesh with 25554 nodes and 120942 tetrahedral elements was created using our customized mesh generator. The meshing process took about 14 seconds on a 1.8GHz Intel Xeon (64bit) PC and building the matrices in Eqs. (9)
took about 11 seconds. Subsequently, we solved Eqs. (9)
for about 1.3 seconds per time step.
In we show a contour plot of the steady state PO2, obtained for the above simulation parameters in X-Y plane that contains the vessel axis. We see a high PO2 for the inflowing blood on the left (starting at 90 mmHg), diminishing as the blood flows across the vessel to a minimum of 64 mmHg indicated by the last contour on the right. In , we plot the PO2 along the interior of the vessel for various values of the blood velocity, hematocrit, and tissue oxygen consumption, respectively. As we increase velocity from 0.4 mm/s to 0.8 mm/s and 1.6 mm/s, we see that the PO2 is increasing as more oxygen is being delivered by the blood per unit time, but the oxygen extraction by the tissue is constant, and thus the oxygen extraction fraction is expected to decrease. Similarly, as we increase the hematocrit from 0 to 0.2 we see that PO2 increases as more oxygen is delivered by the blood. When tissue oxygen consumption is increased, the PO2 in the blood decreases, as expected. The difference between the oxygen flowing into the volume per unit time and the oxygen flowing out should be equal to the tissue volume integral of the oxygen consumption. Assuming the uniform oxygen consumption (OC, in μM/s), the oxygen conservation requires
Model validations with a single vessel: (a) PO2 distribution across an X-Y plane along the vessel; (b-d) PO2 along the vessel for different values of (b) blood velocity, (c) hematocrit, and (d) oxygen consumption.
is the velocity of blood, CT,in
are the total concentrations of oxygen at the input and output boundaries of the blood vessel, respectively, and Vtis
is the tissue volume. The left-hand-side of Eq. (19)
denotes the total oxygen extracted from the vessel per unit time while the right-hand-side denotes the quantity of oxygen consumed by the tissue.
To validate the steady-state solution obtained by our finite element model, we calculated the PO2 distributions under various velocities, hematocrit rates and OC values, and the solutions are shown in . The expected conservation relationship, i.e. Eq. (19)
, holds for all the steady-state simulations.
The second step is to validate the model for dynamic processes. In these simulations, we separately explored the dynamics of oxygen advection within the vessel, flux across the vessel membrane and diffusion in the tissue and compared the calculations with simple analytic solutions. These results are shown in and detailed below.
Fig. 3 Validation of dynamic processes with a single vessel. (a) PO2 profiles along the vessel at 20 ms intervals for isolated advection. The modeled PO2 profile is given by the solid line. The star represents the expected advancement of oxygen pressure profile (more ...)
Oxygen advection within the vessel was isolated from oxygen diffusion by setting the vessel wall permeability Kw
= 0. We started with no oxygen in the vessel and introduced 90 mmHg of oxygen pressure at the input of the vessel ramping down to 0 mmHg at 40 μm into the vessel. shows the oxygen pressure profiles along the vessel at 20 ms time intervals. We see that the oxygen wavefront is advancing almost 8 μm every 20 ms, consistent with the velocity of 0.4 mm/s. It is known that a two-level implicit advection solver produces solutions that underestimate the true velocity at high spatial frequencies and accurately estimate the velocity at low spatial frequencies [26
]. This is evident in with smoothing of the edges at 90 mmHg and 0 mmHg as the oxygen advects along the vessel. We note that this is an intrinsic limitation of the discretization scheme, and although we have implemented the most accurate solver, it places limits on the accuracy of our solution for oxygen content varying with high spatial frequency along a vessel. As noted in Ref. 
, this limitation can be reduced by using a refined time step or smaller grid spacing.
Flux of oxygen across the vessel wall was validated by fixing the vessel PO2
at 90 mmHg, initiating the tissue PO2
at 0 mmHg, setting the tissue oxygen consumption to 0, and increasing the oxygen diffusion coefficient within the tissue by a factor of 1000 such that the system dynamics would be limited by the vessel wall permeability Kw
. In , we plot the tissue temporal evolution of PO2
, which is increasing to 90 mmHg due to oxygen flux across the permeable vessel wall. When we reduce Kw
by a factor of 2, we see that the rate of increase for the tissue PO2
is reduced. The rate of change in the total oxygen content of the tissue volume is given by Fick's law [28
is the number of moles of oxygen in the vessel or tissue, d
is the diameter of the vessel, l
is the vessel length, w
is the vessel wall thickness, and we assume that any oxygen that crosses the vessel wall distributes uniformly throughout the tissue volume effectively immediately (i.e., the diffusion time constant is much shorter than the vessel wall flux time constant). Under these conditions, the solution of Eq. (20)
is a simple exponential. Note that
indicates either tissue (tis
) or vessel (ves
). The solution of Eq. (21)
is shown in and agrees well with the solution of our finite element model.
Oxygen diffusion was validated by setting oxygen consumption to 0 and initiating the tissue PO2
to 0 mmHg except for a single volume element far from the boundaries of the tissue volume that was given a unit quantity of oxygen. This oxygen was then diffused with our finite element model. The early responses of diffusion after 20 (blue), 25 (green), and 30 ms (red) are plotted in and they roughly agree with the analytic solution of the diffusion equation for a homogeneous infinite medium. We also computed the finite element solution reducing the diffusion coefficient by a factor of 2 and plot the results at 40 (blue), 50 (green), and 60 ms (red) and observed the expected DO2t
scaling of the diffusion equation. That is, if we half the diffusion coefficient and double the time, the analytic solution [32
] predicts the same solution.
Tissue oxygen consumption was validated by setting PO2,tis to 90 mmHg throughout the tissue volume, inhibiting oxygen flux across the vessel wall by setting Kw = 0, and increasing oxygen consumption. Our simulation result (not shown) demonstrates a linear decrease in oxygen concentration in the tissue with time, as we expected.
3.2 Application to a microvessel network
After validating the steady state and dynamic solutions of our oxygen advection-diffusion solver, we proceed in solving for the delivery of oxygen to tissue by a realistic microvessel network. We obtain a volumetric image of the microvessel network from the somatosensory cortex of a rat using two-photon microscopy as described in the methods. The volumetric image spanning 230 × 230 × 450 μm3 is shown in . In this image, we identify a 25 μm in diameter surface pial arteriole in the upper right delivering oxygenated blood to the tissue and a draining ascending venule in the upper left of the image, connected by the intervening capillary network.
Fig. 4 Anatomical vessel network extracted from two-photon scans of a rodent somatosensory cortex. (a) Raw image (Media 1, View 1). (b-d) Top view projections of the 3D maps of (b) vessel types: arterioles - red, venules - blue, and capillaries - green, (c) (more ...)
We manually graphed this vessel network in under 2 hours using custom software to accelerate the process, by placing node points and connecting edges along the vessels. The diameter of the vessel at each node point, the blood pressure, and velocity of blood was calculated as described in the methods. The velocity at each node in the vascular graph was used in solving the oxygen advection diffusion equation to obtain the partial pressure of oxygen at every node point in the vessels and throughout the tissue. The vessel in-flow and out-flow points (i.e., the boundary points as mentioned in Section 2.1) were manually labeled and used for the subsequent advection modeling. The generated 3D mesh for this case contains 138685 nodes and 607406 elements. The total time for mesh generation, calculating the matrices and solving for each time step are 136, 62 and 6 seconds, respectively.
In we plot the graphed vascular network identifying the arterioles (red), capillaries (green), and venules (blue). We see that the surface pial arteriole branches twice. The left-hand branch propagates along the surface to the lower left of the image and then descends 450 μm to the bottom of the image without branching. The right-hand branch descends into the tissue at the lower right. This descending arteriole quickly branches a few times leading to a capillary network that connects with an ascending venule on the upper left of the image. The velocity of blood along the vascular network is shown in . The velocity has a maximum of 10 mm/s in the pial arteriole and reaches a minimum of ~1 mm/s in the capillaries and draining venules. The partial pressure of oxygen along the network is graphed in , showing the expected decrease from feeding arteriole to draining venule.
The median velocity, blood pressure, and partial pressure of oxygen versus vessel diameter are shown in . Our pressure distribution compares well with a survey of experimental data [33
]. Our velocity and partial pressure of oxygen, however, drop below experimental data [33
] in the capillaries and venules. This is likely a result of our incomplete vascular graph not accounting for the surrounding microvessels that are providing more blood into the draining vein and delivering more oxygen to the tissue. We will correct this in the future by considering a vascular graph over a larger volume of tissue in order to obtain a more completely connected graph in the central regions.
Fig. 5 The (a) median velocity, (b) blood pressure and (c) partial pressure of oxygen versus vessel diameter compared with experimental data from the literature. The experimental data in panel (a-b) were derived from  and (c) from .
The partial pressure of oxygen within the tissue is shown in in slices at different depths and in a movie of the volumetric data. From these data, we see the high tissue PO2
around the arterioles dropping off rapidly with distance, while the lower PO2
around the capillaries is much more spatially uniform. In we plot the drop in PO2
versus distance from selected arterioles, capillaries, and venules. These results qualitatively agree with experimental in-vivo
The 3D rendering (left, Media 3, View 2) and cross-sections (right) of the PO2 distributions in the vessel and tissue.
The tissue PO2 versus distance from a vessel wall of a few selected (a) arterioles, (b) capillaries and (c) venules. Different colors represent different vessels selected for this test.