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

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.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 [8,9]. Modeling studies of angiogenesis are important in understanding the progression of tumor [10] and its therapy [11][12][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][19][20] and numerical simulation [21][22][23][24][25][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 allto-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 cm 3 , 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.

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.

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 where P O 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 where M 0 represents the maximum consumption under infinite oxygen supply and P 0 represents the partial pressure of oxygen at half-maximal consumption.

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 where i and j are indices of the two end nodes of the vessel, μ is the viscosity of the blood, and R ij and L ij 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, Q ij , satisfies the Ohm's law The Kirchhoff's circuit law describes the conservation of mass X where s i 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 s i = 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 s i 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.

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 S a and blood oxygen partial pressure P b satisfies the Oxygen-hemoglobin dissociation relation where n = 2 is used in this work and P 50 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 where Q is the blood flow rate, α b is the oxygen solubility in blood plasma, H D is the discharge hematocrit, and C 0 is the concentration of hemoglobin-bound oxygen in a fully saturated red blood cell (RBC).

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 @P O ðr; y; sÞ @r j r¼R Rdy 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 P O gradient on the vessel wall. The oxygen flux can also be described by the Kedem-Katalchsky's law [36], where L p is the hydraulic permeability of the vessel wall, P b (s) is the blood oxygen pressure, P O,wall (s) is the circumferential average oxygen partial pressure on the outer surface of the vessel wall, and σ and ΔP 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 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 P b at the inlet of each vessel: (1) At the inlets of the vessel system, P b is given; (2) At the bifurcation points, P b at the inlet of the downstream vessels is inherited from the parent vessel; (3) At collecting junctions, P b at the inlet of the downstream vessel is obtained from its parent vessels as the mixed value by conservation of oxygen.

Model simplification
Since the oxygen partial pressure should be continuous in the whole domain, i.e., the tissue domain and the vessel domain, P b provides a Dirichlet boundary condition on vessel walls for P O 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 where S(x) = R 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(P O ). 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.

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.

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 P O and P b . 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 P n O and P n b at step n, by assuming that the blood oxygen partial pressure P O = P b on vessel walls, Eq (7) allows us to evaluate the flux q n (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 q n (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 q n as which is linearly dependent on P n O and P n b . According to our numerical tests, it is not efficient for convergence to update the flux directly using q n � . Instead, a weighted sum is more efficient, where the weight l ¼ 1 2 is used in our simulation. Once the flux is obtained, Eq (9) is used to update P b in all vessels The fourth-order Runge-Kutta method is used to solve the ODE (13) to obtain P nþ1 b in each vessel segment. In practice, q n � can also be updated in a Gauss-Seidel fashion. Namely, the values of P nþ1 b on the updated nodes can be used in Eq (11) to evaluate the flux. The flux q n+1 is also used to calculate the oxygen source S n+1 , which is required in computing P nþ1 O . It is possible to update P O by fully solving the PDE (10) with given oxygen source S n+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 where 4 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 P n O in tissue can lead to a decrease of the oxygen flux q n+1 (thus the oxygen source S n+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 P O by fully solving the steady state PDE (10) for each iteration step (namely, Δt = 1). 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).

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 4 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 S n k ðxÞ ¼ q n ðx k Þh k dðx À x k Þ on the center lines, where k is the index of the point source at and h k is the step size on the centerline; Then, we use Peskin's numerical δ-function � dð�Þ (see Appendix B) to distribute all point sources onto their neighboring mesh points [40]. Namely, we have where � dð�Þ is the numerical δ-function and i = (i 1 , i 2 , . . ., i d ) 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.

Post-processing
The post-processing is designed to reduce the error introduced by oxygen-source redistribution. This error can be defined as the difference dP ¼ P nþ1 O satisfy the following equations respectively. Therefore, δP satisfies 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 P 0 ), thus M 0 ðP nþ1 O Þ 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 Due to the linearity of Eq (17), δP can be regarded as the linear combination The boundary condition for Eq (18) is that δP k vanishes when x ! 1. In fact, as shown in Appendix C, only a few mesh-sizes away from x k , the error δP k is already negligible. This observation allows us to solve δP k only on local mesh near to x k (e.g., an 8 × 8 mesh in the 2-D case).
The error δP k can be further decomposed into two parts dP k ¼ P k À � P k , where the first part ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ds is the fundamental solution corresponding to the point source δ(x k − x) and the second part � P k is corresponding to the distributed sources where @O k is the boundary of the local domain O 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 � P k . In summary, for each source point x k , δP k is evaluated on a local mesh at the beginning. At each iteration step, the linear combination dP ¼ P k q nþ1 k h k dP k is used to evaluate the error on mesh points close to vessels introduced by oxygen-source redistribution. This error is then removed from � P nþ1 O .

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][42][43][44][45].

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

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.
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 P b 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 P b 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.

PLOS ONE
Fast numerical method for tissue oxygen supply

Numerical solution for the 2D refined vessel network
The numerical solutions of the oxygen partial pressure P O obtained with the refined vessel network with a 1024 × 1024 mesh are shown in Fig 4. The oxygen partial pressure field is obtained

PLOS ONE
Fast numerical method for tissue oxygen supply with a normal tissue consumption (M 0 = 1.7 × 10 −4 cm 3 O 2 /cm 3 /s) (see Fig 4(a)), whereas the partial pressure field is obtained with a reduced tissue consumption (M 0 = 2.89 × 10 −5 cm 3 O 2 / cm 3 /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].

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 M 0 = 2 × 10 −3 cm 3 O 2 /cm 3 /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.

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. 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 P b 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%).

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 E n is computed by the normalized L 2 -norm of the difference between the solutions obtained with an N k mesh and a 2N k 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. The iteration convergence is shown in Fig 8(b), where the relative difference is shown by the normalized L 2 -norm of the difference of the corresponding functions, q and P O , 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

PLOS ONE
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.
In our simulations, the step size along the blood vessels is set to be h k ¼ h 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  [47], in macaques by [41], in rats by [48] and in a human retina model [31]. 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 N v ¼ L v h k , where L v 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), R n is relatively invariant with the mesh size n and the total vessel length L v . 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 10 6 vessel segments, the estimated total time is about 70390 seconds for h = 10μm.

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.

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 P O is assumed to be equal to P b 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 P O on four nearest mesh points to x + outside the vessel are used to fit a linear function P O (x, y)�P b + 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. ϕ(x) satisfies the following constraints �ðxÞ ¼ 0 for jxj > 2; P j even �ðx À jÞ ¼ P j odd �ðx À jÞ ¼ 1 2 for all real x; P j ðx À jÞ�ðx À jÞ ¼ 0 for all real x; P j ð�ðx À jÞÞ 2 ¼ const for all real x: 8 > > > > > > > > > < > > > > > > > > > :

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

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

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.