A fast numerical method for oxygen supply in tissue with complex blood vessel network

Abstract

About the Authors: Yuankai Lu Roles Methodology, Software, Writing – original draft, Writing – review & editing Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China ORCID logo https://orcid.org/0000-0003-3778-4198 Dan Hu Roles Conceptualization, Methodology, Supervision, Writing – review & editing * E-mail: [email protected] Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China Wenjun Ying Roles Resources, Software Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China Abstract Angiogenesis plays an essential role in many pathological processes such as tumor growth, wound healing, and keloid development. [...]existing techniques can hardly offer a complete spatial-temporal picture of oxygen field on microscopic scales. [...]theoretic modeling [18–20] and numerical simulation [21–26] is widely utilized in evaluating oxygen level and studying angiogenesis. [...]in most real applications, it is important to incorporate the complex vessel network structures into the model and numerical algorithms to evaluate the oxygen field in tissue. [...]the corresponding computational cost is huge for tissue with large-scale complex vessel network structure.

About the Authors:

Yuankai Lu

Roles Methodology, Software, Writing – original draft, Writing – review & editing

Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China

ORCID logo https://orcid.org/0000-0003-3778-4198

Dan Hu

Roles Conceptualization, Methodology, Supervision, Writing – review & editing

* E-mail: [email protected]

Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China

Wenjun Ying

Roles Resources, Software

Affiliation: School of Mathematical Sciences, Institute of Natural Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai, China

Abstract

Angiogenesis plays an essential role in many pathological processes such as tumor growth, wound healing, and keloid development. Low oxygen level is the main driving stimulus for angiogenesis. In an animal tissue, the oxygen level is mainly determined by three effects—the oxygen delivery through blood flow in a refined vessel network, the oxygen diffusion from blood to tissue, and the oxygen consumption in cells. Evaluation of the oxygen field is usually the bottleneck in large scale modeling and simulation of angiogenesis and related physiological processes. In this work, a fast numerical method is developed for the simulation of oxygen supply in tissue with a large-scale complex vessel network. This method employs an implicit finite-difference scheme to compute the oxygen field. By virtue of an oxygen source distribution technique from vessel center lines to mesh points and a corresponding post-processing technique that eliminate the local numerical error induced by source distribution, square mesh with relatively large mesh sizes can be applied while sufficient numerical accuracy is maintained. The new method has computational complexity which is slightly higher than linear with respect to the number of mesh points and has a convergence order which is slightly lower than second order with respect to the mesh size. With this new method, accurate evaluation of the oxygen field in a fully vascularized tissue on the scale of centimeter becomes possible.

Figures

Fig 12

Table 2

Table 1

Fig 1

Fig 2

Fig 3

Fig 4

Fig 5

Fig 6

Fig 7

Fig 8

Fig 9

Fig 10

Fig 11

Fig 12

Table 2

Table 1

Fig 1

Fig 2

Citation: Lu Y, Hu D, Ying W (2021) A fast numerical method for oxygen supply in tissue with complex blood vessel network. PLoS ONE 16(2): e0247641. https://doi.org/10.1371/journal.pone.0247641

Editor: Adélia Sequeira, Universidade de Lisboa Instituto Superior Tecnico, PORTUGAL

Received: October 26, 2020; Accepted: February 10, 2021; Published: February 26, 2021

Copyright: © 2021 Lu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the paper and its Supporting information files.

Funding: This work is supported by National Key R&D Program of China (2019YFA0709503), the National Natural Science Foundation of China (Contract no. 11971312, 11771290, and 91630208). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

1 Introduction

Oxygen plays a key role in animal metabolism. Oxygen supply to tissue is mainly achieved by the circulation system of animals. In particular, the efficiency of oxygen delivery is mainly determined by the microcirculation structure. In order to improve the efficiency of their microcirculation structure, animals have developed different physiological processes to modify the geometry and topology of vessel networks, including blood vessel adaptation and angiogenesis [1–3]. In these physiological processes, the oxygen level is the key driving stimulus. For example, poor microcirculation structure results in local tissue hypoxia (low concentration of oxygen), which leads to the production of growth factors for angiogenesis, such as vascular endothelial growth factor (VEGF) [4, 5]. Angiogenesis plays a critical role in many pathological processes, such as wound healing [6], and keloid development [7], and tumor growth [8, 9]. Modeling studies of angiogenesis are important in understanding the progression of tumor [10] and its therapy [11–13]. Accurate evaluation of oxygen supply is also important in generating artificial vascular networks [14, 15] and non-invasively investigating the impact of diseases and therapeutic procedures on the vascular bed [16].

There have been a few techniques for measurement of oxygen level, such as two-photon phosphorescence lifetime microscopy that can be applied to measure oxygen level in vivo [17]. However, existing techniques can hardly offer a complete spatial-temporal picture of oxygen field on microscopic scales. Therefore, theoretic modeling [18–20] and numerical simulation [21–26] is widely utilized in evaluating oxygen level and studying angiogenesis.

The classical Krogh cylinder model roughly describes the oxygen transport from blood vessels to tissues [18, 27]. In this model, evenly spaced capillaries are assumed to be parallel and supply oxygen to a cylindrical tissue domain. Following Krogh’s model, more detailed oxygen consumption mechanisms are taken into account [28]. Coupled models for oxygen delivery including a reasonable oxygen consumption mechanism, a relatively complex vessel network structure, and detailed blood flow in the vessel network are also proposed in later works [19, 29].

Experimental results have indicated disordered spatial distribution of blood vessels [30]. In retinal vascular network, the arterioles and venules reach out from the center to the periphery, forming a roughly radial skeleton, while capillaries form a vessel network that links the arterioles and venules [31]. It is also observed that the microcirculation in tumor tissue presents a chaotic geometry [32]. Therefore, in most real applications, it is important to incorporate the complex vessel network structures into the model and numerical algorithms to evaluate the oxygen field in tissue.

The chaotic vessel geometry brings great challenges in designing effective numerical algorithms. Although existing simulations based on finite difference method and rectangular grids can provide useful insight for tissue oxygen supply, they introduce unrealistic requirement for the geometry of blood vessel networks [2, 22]. A finite element method is also introduced for more general vessel network topology and geometry [28]. This method requires a particular local processing near blood vessels. Meanwhile, the mesh grid is generated according to the vessel geometry. As a result, the corresponding computational cost is huge for tissue with large-scale complex vessel network structure. In more realistic applications, Secomb and Hsu developed a numerical method based on Green’s function, in which each vessel is regarded as a line source of oxygen [26, 32]. Their numerical method can be used to deal with complex vessel geometry. However, their method suffers from the high computational cost due to the all-to-all interaction between all elements. In order to handle large-scale simulations, Welter et al. introduced a fast numerical method based on finite element method [33]. With a uniform grid, they achieved the simulation in a fully vascularized tissue with a domain size of about 0.5 cm3, which was used to study the pathological characteristics of a tumor and its surrounding tissues [34].

In this work, we develop a fast numerical method to evaluate the coupled system for oxygen delivery in tissue, with particular attention paid to the large-scale complex blood vessel network structures. In our model, the blood vessel is also regarded as a line source for oxygen field. A source distribution technique is used so that we can solve the oxygen field on a square mesh. This technique ensures that blood vessels can be freely embedded in the tissue without the need to fit the grid mesh. A post-processing technique is employed to remove the numerical error induced by source distribution, which allows us to use a relatively large mesh size for the square mesh, while sufficient numerical accuracy is maintained. Furthermore, an efficient iteration method is used to deal with the nonlinearity of the coupled system. The convergence order of our new method is slightly smaller than second order with respect to the mesh size and the computational cost is slightly larger than linear with respect to the total number of mesh points. Due to the advantages of this new method, we can accurately evaluate the oxygen field of a three dimensional fully vascularized tissue on the scale of centimeter within 20h.

2 Modeling

The model for oxygen delivery in blood vessels and tissues is well established in previous works [19, 29]. The model mainly includes a partial differential equation (PDE) for the oxygen field, a linear system for blood flow and blood pressure, and a system of ordinary differential equations (ODE) for the blood oxygen.

2.1 Oxygen diffusion in tissue

The oxygen-consuming tissue is a mixture of cells, extracellular matrix and extracellular fluid. Oxygen transport in tissue depends mostly on diffusion. The diffusion constant and solubility of oxygen may vary slightly in different tissue. Here we assume a uniform oxygen diffusivity D and uniform solubility α. Oxygen diffusion in tissue at steady state satisfies the reaction-diffusion equation(1)where PO is the partial pressure of oxygen and M(⋅) is the oxygen-consuming rate in tissue. The oxygen consumption in cells consists of various biochemical processes, which can be described by the Michaelis-Menten equation in general(2)where M0 represents the maximum consumption under infinite oxygen supply and P0 represents the partial pressure of oxygen at half-maximal consumption.

2.2 Blood flow and blood pressure

The blood flow and blood pressure can be computed with Ohm’s law and Kirchhoff’s circuit law that form a system of linear equations. For a small blood vessel, which can be regarded as a cylinder, the blood flow in the vessel can be well approximated by the Poiseuille flow. The conductance of the vessel is(3)where i and j are indices of the two end nodes of the vessel, μ is the viscosity of the blood, and Rij and Lij are the radius and length of the vessel, respectively. It is worth noting that, the effective viscosity μ can be dependent on the vessel radius and blood oxygen level [26, 32]. In this case, the system becomes a nonlinear system that is coupled to the whole system for oxygen delivery.

The blood pressure and blood flow rate from node i to node j, Qij, satisfies the Ohm’s law

The Kirchhoff’s circuit law describes the conservation of mass(4)where si gives the sources and sinks at the inlets and outlets of the vessel network. In most cases, fluid exchange between the blood and tissue may only be a small percent of the total flow (e.g., <0.5% as evaluated from Ref [35]). In this case, we can simply set si = 0 for all inner nodes including the junction points. When this fluid exchange becomes large, such as in the kidney or in particular pathological states, the Starling’s law can be used by introducing continuous flow sinks si along the vessels [35]. In real applications, the boundary conditions at the inlets and outlets can also be replaced by other conditions, such as fixed pressure condition.

2.3 Oxygen flux in vessels

Oxygen carried by blood includes two parts, the minor of which is directly dissolved in the plasma while the major of which is associated with hemoglobin. The relation between oxygen saturation Sa and blood oxygen partial pressure Pb satisfies the Oxygen-hemoglobin dissociation relation(5)where n = 2 is used in this work and P50 is the half-saturated oxygen pressure that may depend on the pH-value of blood. Correspondingly, the oxygen flux f through a cross-section of the vessel includes two parts(6)where Q is the blood flow rate, αb is the oxygen solubility in blood plasma, HD is the discharge hematocrit, and C0 is the concentration of hemoglobin-bound oxygen in a fully saturated red blood cell (RBC).

2.4 Oxygen exchange on vessel walls

Let s be the arc-length parameter along the centerline of a vessel and x(s) be the coordinate of the centerline. Assume the cross-section of the blood vessel is a circle with radius R, then the total oxygen flux q(s) through the blood vessel wall per unit length is(7)where r is the polar radius and θ is the polar angle on the cross-sectional plane (the pole is at x(s)). In this work, we assume that the diffusion constant of oxygen in the blood vessel wall is the same as that in the tissue, in which the vessel wall can be regarded as part of the tissue. Therefore, the oxygen flux can be evaluated from the average PO gradient on the vessel wall. The oxygen flux can also be described by the Kedem-Katalchsky’s law [36],(8)where Lp is the hydraulic permeability of the vessel wall, Pb(s) is the blood oxygen pressure, PO,wall(s) is the circumferential average oxygen partial pressure on the outer surface of the vessel wall, and σ and ΔΠ are the osmotic reflection coefficient and the osmotic pressure difference, respectively. In the Kedem-Katalchsky’s law, precise experimental measurements of multiple parameters (e.g., the thickness of blood vessel wall and the vascular permeability) are used to evaluate the conductivity of oxygen through the blood vessel wall.

Ignoring the oxygen consumption in blood, we obtain the conservation of oxygen(9)

Given the oxygen flux q(s), the above equation becomes a first-order ordinary differential equation (ODE). A boundary condition is required to solve the equation for each vessel. In real applications, we may have Pb at the inlet of each vessel: (1) At the inlets of the vessel system, Pb is given; (2) At the bifurcation points, Pb at the inlet of the downstream vessels is inherited from the parent vessel; (3) At collecting junctions, Pb at the inlet of the downstream vessel is obtained from its parent vessels as the mixed value by conservation of oxygen.

2.5 Model simplification

Since the oxygen partial pressure should be continuous in the whole domain, i.e., the tissue domain and the vessel domain, Pb provides a Dirichlet boundary condition on vessel walls for PO in Eq (1). However, the disordered structure of blood vessel walls brings great difficulties in meshing and numerical simulations to solve the above coupled system. In previous studies [37], the oxygen flux is represented by oxygen sources on the centerlines of vessels. Under this simplification, the governing equation for oxygen supply is defined in the whole domain. In this case, the governing equation becomes(10)where S(x) = ∫q(s)δ(x − x(s))ds is the oxygen source supplied by the vessel.

Note that the vessel radius is relatively small (2 ∼ 3μm) compared to the distance between capillaries (∼100μm). Therefore, this simplification can be a good approximation at least for far field. The oxygen partial pressure may be overestimated near the vessels, which can lead to a slight overestimate of the oxygen source M(PO). Further improvement on such an approximation has been discussed in the work of Ref. [26, 38]. Instead of concentrated oxygen sources on the centerlines of blood vessels, distributed sources using smooth kernel functions are also used in a recent study [39]. This can be effectively used to avoid the weak singularities of the oxygen field.

Now we have obtained a coupled system for oxygen delivery, which includes a PDE (10) on the whole domain, a set of nonlinear ODEs (9) on the vessel centerlines, and a system of linear equations for blood pressures and blood flows in all vessels.

3 Numerical method

In order to develop a fast numerical method to solve the above coupled system, we are left with two main tasks: (1) find an efficient iteration method to deal with the nonlinearity of the coupled system and (2) develop a fast numerical solver for the PDE (10) with complex sources on blood vessel centerlines.

3.1 Nonlinear iteration

Mainly due to the nonlinearity of the ODEs (9) and the unusual form of coupling, general iteration methods, such as the full Newton-like iterative methods, are both hard in code implementation and lack of convergence guarantee for the coupled system. Here we propose an iterative method by decoupling the system and alternatively updating PO and Pb. In particular, we introduce pseudo time step in order to reach the steady state solution. The rest of the section then explains how each (pseudo) time step is solved, namely with a semi-implicit Euler method.

Given and at step n, by assuming that the blood oxygen partial pressure PO = Pb on vessel walls, Eq (7) allows us to evaluate the flux qn(s) numerically. In the 3-D case, we expand the Laplacian operator in series near the vessel first, and then use particular fitting to evaluate the derivative in Eq (7). In the 2-D case, linear fitting on both sides of the vessel can be directly utilized to evaluate the flux qn(s). Detailed discussion on the fitting method is included in Appendix A. Similar treatment can also be used to evaluate the oxygen flux utilizing the Kedem-Katalchsky’s law. In general, we denote the evaluation of qn as(11)which is linearly dependent on and . According to our numerical tests, it is not efficient for convergence to update the flux directly using qn*. Instead, a weighted sum(12)is more efficient, where the weight is used in our simulation. Once the flux is obtained, Eq (9) is used to update Pb in all vessels(13)

The fourth-order Runge-Kutta method is used to solve the ODE (13) to obtain in each vessel segment. In practice, qn* can also be updated in a Gauss-Seidel fashion. Namely, the values of on the updated nodes can be used in Eq (11) to evaluate the flux.

The flux qn+1 is also used to calculate the oxygen source Sn+1, which is required in computing . It is possible to update PO by fully solving the PDE (10) with given oxygen source Sn+1. However, since we still need the nonlinear iteration, it is neither necessary nor efficient to find the accurate solution at each iteration step. Instead, we use the following scheme to update PO (14)where △h is the numerical Laplacian operator and Δt is the temporal step size. This implicit numerical scheme solves the time-dependent reaction-diffusion equation for one time step. When the iteration converges, the solution satisfies the steady state PDE (10).

Note that an increase of the oxygen partial pressure in tissue can lead to a decrease of the oxygen flux qn+1 (thus the oxygen source Sn+1). This implies at least one negative eigenvalue of the coupled system in the sense of linearization. As a result, it may even be not stable to update PO by fully solving the steady state PDE (10) for each iteration step (namely, Δt = ∞). The instability has been observed in our numerical test. From this point of view, a suitable time step size Δt should be selected so that it is both good for iteration stability (Δt is not too big) and good for fast convergence (Δt is sufficiently big).

3.2 PDE solver

Based on the iteration framework above, our remaining task is to numerically solve the PDE (10) efficiently. For the simplified coupled system, we do not consider the detailed geometry of vessel walls. Thus, a square mesh can be easily used in our simulation. We simply use the central difference as the numerical Laplacian △h, where h is the mesh size. This brings great convenience in developing a fast solver. The system in this paper contains no fluid convection and maintains a uniform diffusion coefficient. Under this situation, the central differential scheme is equivalent to the finite volume method, which ensures the local mass conservation. When convection becomes significant, we can use the mass-conservative finite volume method instead.

There are two problems left for us to solve: First, the oxygen sources are located on vessel centerlines, we need to distribute the sources onto the square mesh points; Second, in order to conduct a large scale simulation, we need to use a relatively large spatial mesh size while maintaining sufficient numerical accuracy.

In order to distribute the oxygen sources onto the square mesh points, first we discretize the sources to be point sources on the center lines, where k is the index of the point source at xk = (xk,1, xk,2, …, xk,d) (d is the space dimension) and hk is the step size on the centerline; Then, we use Peskin’s numerical δ-function (see Appendix B) to distribute all point sources onto their neighboring mesh points [40]. Namely, we have(15)where is the numerical δ-function and i = (i1, i2, …, id) is the index of the mesh point.

Due to the nonlinearity in the consumption function M(⋅), we use the standard multigrid algorithm combined with the Newton’s method to solve Eq (14). According to our numerical tests, only a few steps of Newton-iteration are sufficient to make the numerical error small enough.

Notably, the above redistribution of oxygen-sources can lead to numerical error in the solution. As a result of using Peskin’s numerical δ-function, the numerical error induced by source redistribution is local—mainly on the local mesh points to which the oxygen sources are distributed (see Appendix C). Nevertheless, the local oxygen field around the blood vessels must be sufficiently accurate for evaluating the oxygen flux using Eq (11). Therefore, without further improvement, we can only use a small spatial mesh size h to perform simulations.

Next, we introduce a post-processing technique to reduce the local error induced by the redistribution of oxygen-sources. With this post-processing step, we are able to use a relatively large mesh size h (e.g., be comparable to or even larger than vessel diameters) while maintaining a sufficiently small numerical error.

3.3 Post-processing

The post-processing is designed to reduce the error introduced by oxygen-source redistribution. This error can be defined as the difference , where and satisfy the following equationsrespectively. Therefore, δP satisfies(16)where .

Eq (16) appears to be nonlinear due to the nonlinearity in δM. However, in real applications, this term is always very small and we can neglect it. The reason that δM is small is double fold: For mesh points near a vessel, the oxygen partial pressure is relatively high (much bigger than P0), thus is very small; whereas for mesh points away from all vessels, δP becomes very small thanks to the good property of the numerical δ-function. Therefore, we can evaluate δP by neglecting the nonlinear term(17)

Due to the linearity of Eq (17), δP can be regarded as the linear combination , where δPk satisfies(18)

The boundary condition for Eq (18) is that δPk vanishes when x → ∞. In fact, as shown in Appendix C, only a few mesh-sizes away from xk, the error δPk is already negligible. This observation allows us to solve δPk only on local mesh near to xk (e.g., an 8 × 8 mesh in the 2-D case).

The error δPk can be further decomposed into two parts , where the first part is the fundamental solution corresponding to the point source δ(xk − x) and the second part is corresponding to the distributed sources(19) (20)where ∂Ωk is the boundary of the local domain Ωk. Note that the coefficient matrix to numerically solve the linear system (19) on the local mesh is independent of k and the size of the coefficient matrix is small (e.g., 64 × 64). We can find and save the inverse of the coefficient matrix numerically. Then a matrix-vector multiplication is used to find .

In summary, for each source point xk, δPk is evaluated on a local mesh at the beginning. At each iteration step, the linear combination is used to evaluate the error on mesh points close to vessels introduced by oxygen-source redistribution. This error is then removed from .

4 Parameters for simulation

The model parameters used in this work are listed in Table 1. The data are adapted from previous experiments on different animal retinas [41–45].

[Figure omitted. See PDF.]

Table 1. Model parameters.

https://doi.org/10.1371/journal.pone.0247641.t001

5 Numerical results

5.1 Blood vessel network structures

Four model structures of blood vessel networks are used in our simulation (see Fig 1). The simple cobweb structure (Fig 1(a)) in 2D and the single vessel (Fig 1(c)) in 3D are used for model validation and numerical convergence analysis. The 2D and 3D refined structures (Fig 1(b) and 1(d)) are adapted from real experimental measurements on a mouse retina. They are further used for analyzing the efficiency of oxygen supply. In this study, the no-flux boundary condition is used in solving the oxygen field with Eq (1).

[Figure omitted. See PDF.]

Fig 1. Model vessel structures.

(a) A cobweb vessel structure mimicking the main branches of retinal vessel networks. There are six inlets and six outlets for blood flow near the center of the network. The inlets and outlets are in a spaced arrangement. The arrows show the direction of blood flow. (b) A refined vessel network with 4306 blood vessels. The network structure is adapted from a real retinal vessel network measured in experiment by stretching and symmetrical extensions. There are four flow inlets and flow outlets near the center. (c)A single vessel embedded in a tissue cube. (d) A refined 3D vessel network with 7815 blood vessels. The intermediate and deep capillary plexi layers in retina are reconstructed by projecting the above 2D network onto two spherical shells. The diameter of the two spheres for the deep capillary plexi layer, the intermediate layer, and the choroid layer is 2 mm, 2.1 mm, and 2.2 mm, respectively. The widths of the lines in each figure show the radii of blood vessels. The unit for the axes is micrometer for (c), and millimeter for others.

https://doi.org/10.1371/journal.pone.0247641.g001

5.2 Numerical solution for the cobweb vessel network

In Fig 2(a), we show the numerical solution of the oxygen partial pressure obtained with a 1024 × 1024 square mesh for the cobweb vessel structure. Since the distance between neighboring vessels is much bigger than that of real microvasculature, tissue away from the vessels is in a hypoxic state, i.e., the oxygen partial pressure is very low. The width of the well irrigated region for each vessel depends on the blood flow in it, being roughly on the scale of 100 μm. Note that the gap between neighboring capillaries in normal tissue is also roughly on the scale of 100 μm, which is much smaller than that in the cobweb structure. With such a high capillary density, normal tissue can be well irrigated.

[Figure omitted. See PDF.]

Fig 2. Oxygen partial pressure and its numerical error.

(a) Numerical solution of the oxygen partial pressure obtained with a 1024 × 1024 square mesh. (b) Numerical error of the oxygen partial pressure calculated by the difference between the solution obtained with a 1024 × 1024 mesh and a 2048 × 2048 mesh. The unit for oxygen partial pressure is mmHg.

https://doi.org/10.1371/journal.pone.0247641.g002

The numerical error of the oxygen partial pressure is shown in Fig 2(b). The error is calculated by comparing the solution obtained with a 1024 × 1024 mesh and a 2048 × 2048 mesh. The error is smaller than 1 mmHg. This suggests that with a mesh size of about 5 μm, the numerical error can be smaller than one percent.

In Fig 3, we show the blood oxygen partial pressure Pb and oxygen flux q on all vessels. Because the blood vessels of the most inner loop are very small, they have a big resistance to blood flow. As a result, the blood flow in these vessels is very small, which leads to a fast decay in blood oxygen partial pressure along the flow direction (see Fig 3(a)). From Fig 3(c), we can see that the blood oxygen partial pressure Pb decreases along each vessel segment as the oxygen diffuses from vessel to tissue. Similarly, in Fig 3(b) and 3(d), we show the oxygen flux q. From Fig 3(b), we can see that the flux q reaches minimums at most cross-links of the vessel segments. This is because that more than one vessels share the oxygen supply to the local tissue at the cross-links. On the contrary, from Fig 3(d), we can see that the oxygen flux in vessel a and d reaches the maximums at their one or two ends. This is due to the vacancy of vessels at the center and the outer domain. At the end of the vessel segment b, the blood oxygen partial pressure becomes a constant and the flux becomes zero, because the oxygen partial pressure is higher in the surrounding tissue than the blood inside the vessel. Note that this unphysiological behavior is a consequence of the cobweb vessel structure. Here we set the negative flux to be zero. This setting does not significantly change the tissue oxygen level and blood oxygen concentration in downstream vessels since this piece of vessel is always short.

[Figure omitted. See PDF.]

Fig 3. Blood oxygen partial pressure and oxygen fluxes from the vessels.

(a) Blood oxygen partial pressure Pb on the vessels. (b) Oxygen fluxes q from the vessels to tissue. (c-d) Pb and q on four marked vessel segments in (a). The x-axis denotes the distance from the inlet for each vessel (i.e., arc-length coordinate). The unit of Pb is mmHg. The unit of q is 10−9cm3O2/cm/s.

https://doi.org/10.1371/journal.pone.0247641.g003

5.3 Numerical solution for the 2D refined vessel network

The numerical solutions of the oxygen partial pressure PO obtained with the refined vessel network with a 1024 × 1024 mesh are shown in Fig 4. The oxygen partial pressure field is obtained with a normal tissue consumption (M0 = 1.7 × 10−4cm3O2/cm3/s) (see Fig 4(a)), whereas the partial pressure field is obtained with a reduced tissue consumption (M0 = 2.89 × 10−5cm3O2/cm3/s) (see Fig 4(c)). In both cases, we can see that the oxygen partial pressure directly reflects the refined vessel structure. In Fig 4(b) and 4(d), we statistically analyze the oxygen fields inside the circles shown in Fig 4(a) and 4(c), respectively. The distance between the circle and the outer vessels is greater than 200 μm. Hence the boundary effects is insignificant. For the normal consumption case as shown in Fig 4(a) and 4(b), tissue in the outer domain has a very low oxygen supply. This is due to the large stretch in the outer domain when we generate the refined vessel network, which significantly affects the oxygen supply in two folds: First, it increases the irrigated area of the corresponding blood vessels in the outer domain; Second, it increases the resistance of these vessels, thus decreases the blood flow in these vessels. For the reduced consumption case as shown in Fig 4(c) and 4(d), all tissue inside the circle is well irrigated. This statistical results of oxygen field under various tissue consumption are consistent with that reported in previous works [21, 46].

[Figure omitted. See PDF.]

Fig 4. Oxygen partial pressure obtained with the refined vessel network.

The partial pressure for normal tissue consumption and reduced tissue consumption are shown in (a) and (c), respectively. The tissue inside the white circles shown in (a) and (c) are used to statistically evaluate the area percentages of tissue with particular partial pressure of oxygen. The statistical results are shown in Figure (b) and (c), respectively. The unit of the x− and y− axes in (a) and (c) is mm. The unit of oxygen partial pressure is mmHg.

https://doi.org/10.1371/journal.pone.0247641.g004

5.4 Numerical solution for the 3D refined vessel network

The 3D refined structure with 7815 vessels is adapted from the intermediate and deep capillary plexi of a mouse retina and used for simulation of the oxygen field with a 512 × 512 × 512 mesh. The tissue consumption rate M0 = 2 × 10−3cm3O2/cm3/s and total blood flow Q = 10nl/s are adapted from experimental data to simulate the real situation [47, 48]. The simulated oxygen field at different layers is displayed in Fig 5. While low oxygen partial pressure can be observed in the layers away from the capillary plexus, the capillary-rich areas are well irrigated, which agrees well with the experimental data. Note that the irrigation efficiency depends on the total blood flow rate Q and the vessel density.

[Figure omitted. See PDF.]

Fig 5. The oxygen partial pressure at different layers of the retina.

The oxygen profiles on the layers corresponding to 0, 25, and 50 percent of the retina depth are shown in the Figure (a), (b), and (c), respectively. The unit of the axes is mm and the unit of oxygen partial pressure is mmHg.

https://doi.org/10.1371/journal.pone.0247641.g005

5.5 Model validation

The comparison between our results and previous experimental and simulation results are shown in Fig 6. The horizontal axis represents the normalized retinal depth and the vertical axis represents the oxygen partial pressure. Despite the large differences in the independent experimental profiles, the profiles indeed share similar features, e.g., relatively high oxygen partial pressure is observed near the choroid surface and the deep capillary plexi, which is attributed to the high local vessel density. Meanwhile, a ‘peak’ can be observed in the middle, which corresponds to the intermediate capillary plexi layer. Despite the differences induced by particular parameter settings and experimental conditions, the simulation results are relatively consistent with the experimental data.

[Figure omitted. See PDF.]

Fig 6. Comparison between retina oxygen partial pressure profiles simulated by our model for the mouse retina and those experimentally measured in mice by [47], in macaques by [41], in rats by [48] and in a human retina model [31].

DCP: deep capillary plexi, ICP: intermediate capillary plexi.

https://doi.org/10.1371/journal.pone.0247641.g006

The comparison between our simulation results and that obtained with the numerical method developed in Ref. [26] are shown in Fig 7. Both simulations are performed on the same single-vessel domain as shown in Fig 1(C). The edge length of the cube is 640 μm and the vessel radius is 10 μm. The oxygen profile in Fig 7(a) was obtained from a line through the cube center perpendicular to the vessel. All curves exhibit a similar diffusion distance that agrees well with the data under physiological conditions. Our results show good consistency with that obtained with the method developed in Ref. [26]. The blood oxygen partial pressure Pb obtained with different mesh size is shown in Fig 7(b). A more detailed illustration of the error in oxygen field is shown in Fig 12. According to these results, a mesh size of 10μm is sufficient to achieve a relatively low numerical error (1%).

[Figure omitted. See PDF.]

Fig 7. Model validation and convergence analysis.

(a) The profiles of oxygen partial pressure on a line perpendicular to the vessel, where x = 0 represents cube center. The oxygen partial pressure inside the vessel (−10, 10) are set to be equal to the blood oxygen partial pressure Pb. The line “Secomb” is obtained from the method developed in Ref. [26] with a mesh size of 10μm. All simulations are performed with tissue consumption M0 = 2.0 × 10−3cm3O2/cm3/s and total blood inflow Q0 = 0.05nl/s. (b) The blood oxygen partial pressure Pb along the blood vessel.

https://doi.org/10.1371/journal.pone.0247641.g007

5.6 Convergence analysis and efficiency analysis

We use the cobweb vessel network in 2D and the simple single-vessel network in 3D to numerically study the convergence of our new method. The results for mesh refinement test are shown in Fig 8(a). The relative error En is computed by the normalized L2-norm of the difference between the solutions obtained with an Nk mesh and a 2Nk mesh, where k is the dimension of the model. The numerical convergence order is about 1.7 for both the 2-D and 3-D cases. For the 3-D case, the relative error reaches at about 1% with a mesh size of 10μm, which is consistent with the result shown in and Fig 7.

[Figure omitted. See PDF.]

Fig 8. Numerical convergence analysis.

(a)Mesh refinement test. The stars show the numerical error. The fitted lines show a convergence order of 1.73 for the 2D case and 1.71 for the 3D case with respect to mesh size. The relative error is computed by the normalized L2-norm of PO. (b) The decay of the relative difference between two iteration steps. The relative difference is computed by the normalized L2-norm of q. (c) Iteration numbers for different mesh size.

https://doi.org/10.1371/journal.pone.0247641.g008

The iteration convergence is shown in Fig 8(b), where the relative difference is shown by the normalized L2-norm of the difference of the corresponding functions, q and PO, between two iteration steps. We can see that the relative differences have a fast decay at the beginning, followed by a relatively slower linear convergence. In Fig 8(c), we show the iteration numbers required to make the relative difference smaller than 1 × 10−3. We can see that the iteration numbers are comparable for the simple and refined vessel network structures, while both of them increase slightly with the mesh size N. The increase in iteration numbers is mainly due to the choice of the temporal step size Δt in Eq (14). As we have discussed above, a large Δt is helpful for convergence while a too large Δt can induce instability in the iteration. According to our numerical test, the optimal Δt decreases with N. For example, Δt ≈ 5000 is used for N = 1024 in our simulation, while Δt ≈ 3500 is used for N = 2048. The behavior in the 3-D case is similar.

The average number of Newton iterations required for each full iteration step is about 5 ∼ 6. Therefore, the time cost is about 40 seconds to find the solution with a relative error smaller than 1% on a 1024 × 1024 mesh with a standard 3.7 GHz personal computer, whereas about 170 seconds on a 2048 × 2048 mesh. The iteration process for the 3-D simulation shares the similar convergence rate to that for the 2-D case. In our simulation, the complete 3-D simulation requires about 726 seconds on a 256 × 256 × 256 mesh and about 8200 seconds on a 512 × 512 × 512 mesh. This time increase from 2-D to 3-D comes mainly from the increase in solving the PDE.

The average time cost for each Newton iteration is analyzed in Fig 9. At the beginning of each iteration, three to six Newton iterations are required for each full iteration step. In each Newton iteration step, two to three V-cycles of multigrid iterations are required to solve the oxygen partial pressure field. After a few full iteration steps, only two to three Newton iteration steps, each with only one V-cycle, are sufficient to make the error small enough for each iteration step. In Fig 9(a) and 9(b), we show the average time cost for solving the oxygen field in tissue (PDE-time), solving the blood oxygen partial pressure (ODE-time), and post-processing (P-time) in each Newton iteration step, tested on a standard 3.7 GHz personal computer.

[Figure omitted. See PDF.]

Fig 9. Numerical efficiency analysis.

(a) and (b) Average time cost for each iteration step in the 2-D and 3-D cases, respectively. The black solid lines, the blue dotted lines, and the red dashed lines show the average time cost for solving the PDE, the ODEs, and the post-processing in each Newton iteration step, respectively. For the 2-D case, the circled and squared lines show the time cost for the refined vessel network and the simple cobweb network, respectively; whereas for the 3-D case, the circled and squared lines show the time cost for the refined retina network and the single vessel system, respectively. (c) The time ratio between the OP-time and PDE-time for the 2-D and 3-D cases. (d) The normalized time ratio for different systems.

https://doi.org/10.1371/journal.pone.0247641.g009

In our simulations, the step size along the blood vessels is set to be , where L is the domain size and n is the mesh size on one direction. For a given vessel network structure, the total number of discrete nodes on the blood vessels is proportional to n. Hence the sum of the ODE-time and P-time (OP-time) is also proportional to n. Note that for the standard multigrid method, the computation complexity and the PDE-time are O(n2logn) for 2-D cases and O(n3logn) for 3-D cases, which increase much faster than the OP-time (see Fig 9(a) and 9(b)). As shown in Fig 9(c), for large mesh size n, the OP-time is only a small fraction of the PDE-time for the refined retina vessel networks.

In many real applications, simulations are required to perform in fully-vascularized 3-D tissue domains. Thus it is of interest to estimate the OP-time required in such cases. Obviously, the OP-time increases linearly with the total length of blood vessels. In principle, the OP-time is proportional to the total number of discrete grids on the blood vessels , where Lv is the total length of blood vessels. In order for a fair comparison for different simulations, we define the normalized time ration by

As shown in Fig 9(d), Rn is relatively invariant with the mesh size n and the total vessel length Lv. This suggests that the OP-time is really proportional to the total vessel length for given mesh size n.

For a fully vascularized tissue with fixed vessel density, the total vessel length is proportional to the volume V of the tissue domain. As a result, for given mesh step size h, the ratio between the OP-time and the PDE-time is almost invariant with the domain size. Using the vessel density suggested in Ref. [33], the total vessel length in a cube with the edge length of 5.12 mm is about 13422 mm, which is 13.3 times of that of the 3-D retina network (about 1012 mm) in this work. Therefore, for h ≤ 20μm, the ratio between the OP-time and PDE-time is about smaller than 1 for fully vascularized tissue (see Fig 9(c)). In particular, for a fully vascularized (5.12mm)3 cube, the OP-time is about 789 seconds for h = 10μm, and the total time cost is about 8000 seconds; For a fully vascularized (10.24mm)3 cube, which contains about 106 vessel segments, the estimated total time is about 70390 seconds for h = 10μm.

6 Conclusions and discussions

Oxygen delivery in tissue plays an important role in many physiological processes such as angiogenesis, blood flow regulation, and blood vessel adaptation. Repeatedly evaluating the oxygen field in tissue is the key bottleneck that limits the large scale modeling and simulation of these important processes. In this work, a fast numerical method is developed for the computation of the nonlinear coupled system of oxygen consumption, oxygen diffusion, and oxygen delivery in blood vessels with complex network structures.

Our fast numerical method involves an implicit finite-difference method for solving the partial differential equation of oxygen partial pressure with a square mesh. The key techniques we used include (1) the Peskin’s numerical δ-function to distribute the oxygen sources onto mesh points and (2) the post-processing to remove the numerical error induced by the distribution of oxygen sources. With these techniques, relatively large spacial mesh size can be used while sufficient numerical accuracy is maintained. The computational complexity is slightly bigger than linear with respect to the number of mesh points, taking into account the increase in iteration steps for refined meshes. The convergence of numerical error is slightly less than second order with respect to the mesh size. For three-dimensional simulations, the numerical error can be controlled to be about 1% with a mesh step size of h = 10μm.

Although we have not performed a 3-D simulation of fully vascularized tissue because we do not have a suitable 3-D vascular structure, our numerical tests show that for given step size h ≤ 20μm, the simulation time is mainly costed in the multigrid method for solving the oxygen field. A large scale simulation in a fully vascularized (10.24mm)3 cube can be achieved within 20h for h = 10μm. Moreover, the natural extension of our method from serial to parallel also leaves abundant possibilities for further applications on various large organs.

Nevertheless, we can see that the total number of iterations used in our current simulations is still large. Hence, to pave the way for more realistic three-dimensional simulations with complex blood vessel networks, better iteration strategies should be explored to reduce the time cost in our method.

7 Appendix

A Oxygen flux through blood vessel walls

For the two-dimensional case, oxygen flux is evaluated on both side of the vessel walls separately. The oxygen partial pressure PO is assumed to be equal to Pb on vessel walls. In order to evaluate the oxygen flux at a discrete node x = (x, y) on the vessel center line, we first numerically calculate the unit normal direction n by central difference. Then, we define two corresponding points x± = x±r n on the vessel wall (see Fig 10), where r is the vessel radius. Next, the oxygen partial pressures PO on four nearest mesh points to x+ outside the vessel are used to fit a linear function PO(x, y)≈Pb + G ⋅ (x − x+) in the sense of least squares, where G is the gradient to be fitted. Finally, the oxygen flux is given by q = Dα G ⋅ n.

[Figure omitted. See PDF.]

Fig 10.

(a) Evaluation of oxygen flux in 2D. The red solid lines represent the vessel wall. The red dashed line shows the center line of the blood vessel, the crosses on which show the discrete nodes. The blue (green) solid circle shows x+ (x−), whereas the blue (green) solid squares show the nearest mesh points used to fit the linear function and evaluate the flux. (b) Evaluation of oxygen flux in 3D. The red dashed line shows the center line of the blood vessel. The black circle shows the cross-section of blood vessel, whereas the blue (green) solid squares show the nearest mesh points used to fit the solution and evaluate the flux.

https://doi.org/10.1371/journal.pone.0247641.g010

For the three-dimensional case, the oxygen field near the blood vessel can be described by the Laplace equation under cylindrical coordinate system.where r0 denotes the radius of the blood vessel, z0 and Pb are z coordinate of one discrete node on vessel center line and the corresponding blood oxygen partial pressure, respectively. Then the series expansion of the solution can be writen as,

By keeping the first term and applying the boundary condition, the solution can be written as

In our numerical tests, a few nearest mesh points are used fit the function in the sense of least squares. The distance between the fitting point and the point on center line is less than r0 + h. Finally, the oxygen flux is obtained by q = 2πDαc1.

B Peskin’s numerical δ-function

Peskin’s numerical δ-function is defined aswhere d is the dimension of x, h is the mesh size, and xi are the Cartesian components of x. The one-dimensional continuous function ϕ(x) is given by ϕ(x) satisfies the following constraints

C Localization of δPk

An example of δPk defined in Eq (18) on the local mesh is illustrated in Fig 11, where the mesh size is h = 0.002, xk = (0.5 + h/2, 0.5 + h/2) is set at the center of the local mesh. It can be clearly seen that δPk decays very rapidly. When |x − xk|>2h, the error is already negligible.

[Figure omitted. See PDF.]

Fig 11. The relative error δPk on a local mesh.

https://doi.org/10.1371/journal.pone.0247641.g011

D Detail error analysis for single vessel

A detailed error analysis for this classical Krogh model is shown in Fig 12. In order to estimate the numerical error induced by the spatial discretization, we performed simulations with mesh sizes of 5,10,20 and 40 μm, respectively. The result of 5-μm mesh was used as a reference to determine the errors of those obtained with other mesh sizes. It can be observed that the errors are mainly distributed in the vessel and surrounding tissue. A space step of 40 μm will bring about 10% numerical error in the blood vessel, where the case of 10 μm only has an error less than 1%.

[Figure omitted. See PDF.]

Fig 12. Detailed oxygen field error analysis for the single vessel model.

(a) The oxygen profile on the plane of the vessel under the fine grid (h = 5μm). The numerical error for the h = 10, h = 20 and h = 40 are shown in (b-d), respectively, while the result of 5-μm mesh were used as a standard. The unit of the axes is μm and the unit of oxygen partial pressure is mmHg.

https://doi.org/10.1371/journal.pone.0247641.g012

E Inflow and outflow conditions

The inflow and outflow conditions of the four vascular systems used in in this work are listed in Table 2. When there are multiple inlets and outlets, they share the same blood flow rate.

[Figure omitted. See PDF.]

Table 2.

https://doi.org/10.1371/journal.pone.0247641.t002

Supporting information

S1 File.

https://doi.org/10.1371/journal.pone.0247641.s001

(ZIP)

Citation: Lu Y, Hu D, Ying W (2021) A fast numerical method for oxygen supply in tissue with complex blood vessel network. PLoS ONE 16(2): e0247641. https://doi.org/10.1371/journal.pone.0247641

References

1. Secomb TW, Alberding JP, Hsu R, Dewhirst MW, Pries AR. Angiogenesis: an adaptive dynamic biological patterning problem. PLoS computational biology. 2013;9(3):e1002983. pmid:23555218

2. Bartha K, Rieger H. Vascular network remodeling via vessel cooption, regression and growth in tumors. Journal of Theoretical Biology. 2006;241(4):903–918. pmid:16545398

3. Gorodnova N, Kolobov A, Mynbaev O, Simakov S. Mathematical modeling of blood flow alteration in microcirculatory network due to angiogenesis. Lobachevskii Journal of Mathematics. 2016;37(5):541–549.

4. Dewhirst MW, Cao Y, Moeller B. Cycling hypoxia and free radicals regulate angiogenesis and radiotherapy response. Nature Reviews Cancer. 2008;8(6):425. pmid:18500244

5. Stone J, Itin A, Alon T, Pe’Er J, Gnessin H, Chan-Ling T, et al. Development of retinal vasculature is mediated by hypoxia-induced vascular endothelial growth factor (VEGF) expression by neuroglia. Journal of Neuroscience. 1995;15(7):4738–4747. pmid:7623107

6. Knighton D, Silver I, Hunt T. Regulation of wound-healing angiogenesis-effect of oxygen gradients and inspired oxygen concentration. Surgery. 1981;90(2):262–270. pmid:6166996

7. Beer TW, Baldwin HC, Goddard JR, Gallagher PJ, Wright DH. Angiogenesis in pathological and surgical scars. Human pathology. 1998;29(11):1273–1278. pmid:9824106

8. Kerbel RS. Tumor angiogenesis. New England Journal of Medicine. 2008;358(19):2039–2049. pmid:18463380

9. Kuznetsov M, Kolobov A. The Influence of Chemotherapy on the Progression of a Biclonal Tumor: Analysis Using Mathematical Modeling. Biophysics. 2019;64(2):279–292.

10. Kuznetsov M, Gorodnova N, Simakov S, Kolobov A. Multiscale modeling of angiogenic tumor growth, progression, and therapy. Biophysics. 2016;61(6):1042–1051.

11. Kuznetsov M, Kolobov A. Optimization of Dose Fractionation for Radiotherapy of a Solid Tumor with Account of Oxygen Effect and Proliferative Heterogeneity. Mathematics. 2020;8(8):1204.

12. Kuznetsov M. Mathematical Modeling Shows That the Responseof a Solid Tumor to Antiangiogenic TherapyDepends on the Type of Growth. Mathematics. 2020;8(5):760.

13. Kuznetsov MB, Gubernov VV, Kolobov AV. Analysis of anticancer efficiency of combined fractionated radiotherapy and antiangiogenic therapy via mathematical modelling. Russian Journal of Numerical Analysis and Mathematical Modelling. 2018;33(4):225–242.

14. Köppl T, Vidotto E, Wohlmuth B. A 3D-1D coupled blood flow and oxygen transport model to generate microvascular networks. International journal for numerical methods in biomedical engineering. 2020;36(10):e3386. pmid:32659047

15. Schneider M, Reichold J, Weber B, Székely G, Hirsch S. Tissue metabolism driven arterial tree generation. Medical image analysis. 2012;16(7):1397–1414. pmid:22705288

16. Dewhirst MW, Secomb TW. Transport of drugs from blood vessels to tumour tissue. Nature Reviews Cancer. 2017;17(12):738–750. pmid:29123246

17. Lecoq J, Parpaleix A, Roussakis E, Ducros M, Houssen YG, Vinogradov SA, et al. Simultaneous two-photon imaging of oxygen and blood flow in deep cerebral vessels. Nature medicine. 2011;17(7):893. pmid:21642977

18. Krogh A. The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue. The Journal of physiology. 1919;52(6):409–415. pmid:16993405

19. Popel AS. Theory of oxygen transport to tissue. Critical reviews in biomedical engineering. 1989;17(3):257. pmid:2673661

20. Lücker A, Weber B, Jenny P. A dynamic model of oxygen transport from capillaries to tissue with moving red blood cells. American Journal of Physiology-Heart and Circulatory Physiology. 2014;. pmid:25398979

21. Secomb TW, Hsu R, Dewhirst M, Klitzman B, Gross J. Analysis of oxygen transport to tumor tissue by microvascular networks. International Journal of Radiation Oncology* Biology* Physics. 1993;25(3):481–489. pmid:8436527

22. Goldman D, Popel AS. A computational study of the effect of capillary network anastomoses and tortuosity on oxygen transport. Journal of Theoretical Biology. 2000;206(2):181–194. pmid:10966756

23. Goldman D, Popel AS. A computational study of the effect of vasomotion on oxygen transport from capillary networks. Journal of Theoretical Biology. 2001;209(2):189–199. pmid:11401461

24. Zhu TC, Liu B, Penjweini R. Study of tissue oxygen supply rate in a macroscopic photodynamic therapy singlet oxygen model. Journal of biomedical optics. 2015;20(3):038001. pmid:25741665

25. Liu D, Wood N, Witt N, Hughes A, Thom S, Xu X. Computational analysis of oxygen transport in the retinal arterial network. Current eye research. 2009;34(11):945–956. pmid:19958111

26. Secomb TW, Hsu R, Park EY, Dewhirst MW. Green’s function methods for analysis of oxygen delivery to tissue by microvascular networks. Annals of biomedical engineering. 2004;32(11):1519–1529. pmid:15636112

27. Kreuzer F. Oxygen supply to tissues: the Krogh model and its assumptions. Experientia. 1982;38(12):1415–1426. pmid:7151956

28. D’Angelo C. Multiscale modelling of metabolism and transport phenomena in living tissues. EPFL; 2007.

29. Goldman D. Theoretical models of microvascular oxygen transport to tissue. Microcirculation. 2008;15(8):795–811. pmid:18608981

30. Cassot F, Lauwers F, Fouard C, Prohaska S, LAUWERS-CANCES V. A novel three-dimensional computer-assisted method for a quantitative study of microvascular networks of the human cerebral cortex. Microcirculation. 2006;13(1):1–18. pmid:16393942

31. Causin P, Guidoboni G, Malgaroli F, Sacco R, Harris A. Blood flow mechanics and oxygen transport and delivery in the retinal microcirculation: multiscale mathematical modeling and numerical simulation. Biomechanics and modeling in mechanobiology. 2016;15(3):525–542. pmid:26232093

32. Grimes DR, Kannan P, Warren DR, Markelc B, Bates R, Muschel R, et al. Estimating oxygen distribution from vasculature in three-dimensional tumour tissue. Journal of The Royal Society Interface. 2016;13(116):20160070. pmid:26935806

33. Welter M, Bartha K, Rieger H. Emergent vascular network inhomogeneities and resulting blood flow patterns in a growing tumor. Journal of theoretical biology. 2008;250(2):257–280. pmid:17996256

34. Welter M, Fredrich T, Rinneberg H, Rieger H. Computational model for tumor oxygenation applied to clinical data on breast tumor hemoglobin concentrations suggests vascular dilatation and compression. PloS one. 2016;11(8):e0161267. pmid:27547939

35. Scallan J, Huxley VH, Korthuis RJ. Capillary fluid exchange: regulation, functions, and pathology. In: Colloquium Lectures on Integrated Systems Physiology- From Molecules to Function. vol. 2. Morgan & Claypool Publishers; 2010. p. 1–94.

36. Cattaneo L, Zunino P. A computational model of drug delivery through microcirculation to compare different tumor treatments. International journal for numerical methods in biomedical engineering. 2014;30(11):1347–1371. pmid:25044965

37. Hsu R, Secomb TW. A Green’s function method for analysis of oxygen delivery to tissue by microvascular networks. Mathematical biosciences. 1989;96(1):61–78. pmid:2520192

38. Pozrikidis C, Farrow D. A model of fluid flow in solid tumors. Annals of biomedical engineering. 2003;31(2):181–194. pmid:12627826

39. Koch T, Schneider M, Helmig R, Jenny P. Modeling tissue perfusion in terms of 1d-3d embedded mixed-dimension coupled problems with distributed sources. Journal of Computational Physics. 2020;410:109370.

40. Peskin CS. The immersed boundary method. Acta numerica. 2002;11:479–517.

41. Birol G, Wang S, Budzynski E, Wangsa-Wirawan ND, Linsenmeier RA. Oxygen distribution and consumption in the macaque retina. American Journal of Physiology-Heart and Circulatory Physiology. 2007;. pmid:17557923

42. Hardarson SH, Harris A, Karlsson RA, Halldorsson GH, Kagemann L, Rechtman E, et al. Automatic retinal oximetry. Investigative ophthalmology & visual science. 2006;47(11):5011–5016. pmid:17065521

43. Williamson TH, Rumley A, Lowe G. Blood viscosity, coagulation, and activated protein C resistance in central retinal vein occlusion: a population controlled study. British Journal of Ophthalmology. 1996;80(3):203–208. pmid:8703856

44. Riva CE, Grunwald JE, Sinclair SH, Petrig B. Blood velocity and volumetric flow rate in human retinal vessels. Investigative ophthalmology & visual science. 1985;26(8):1124–1132. pmid:4019103

45. Funk R. Blood supply of the retina. Ophthalmic research. 1997;29(5):320–325. pmid:9323723

46. Klitzman B, Popel AS, Duling BR. Oxygen transport in resting and contracting hamster cremaster muscles: experimental and theoretical microvascular studies. Microvascular research. 1983;25(1):108–131. pmid:6835096

47. Yu DY, Cringle SJ. Oxygen distribution in the mouse retina. Investigative ophthalmology & visual science. 2006;47(3):1109–1112. pmid:16505048

48. Shahidi M, Wanek J, Blair NP, Little DM, Wu T. Retinal tissue oxygen tension imaging in the rat. Investigative ophthalmology & visual science. 2010;51(9):4766–4770. pmid:20375336