Voxelized simulation of cerebral oxygen perfusion elucidates hypoxia in aged mouse cortex

Departures of normal blood flow and metabolite distribution from the cerebral microvasculature into neuronal tissue have been implicated with age-related neurodegeneration. Mathematical models informed by spatially and temporally distributed neuroimage data are becoming instrumental for reconstructing a coherent picture of normal and pathological oxygen delivery throughout the brain. Unfortunately, current mathematical models of cerebral blood flow and oxygen exchange become excessively large in size. They further suffer from boundary effects due to incomplete or physiologically inaccurate computational domains, numerical instabilities due to enormous length scale differences, and convergence problems associated with condition number deterioration at fine mesh resolutions. Our proposed simple finite volume discretization scheme for blood and oxygen microperfusion simulations does not require expensive mesh generation leading to the critical benefit that it drastically reduces matrix size and bandwidth of the coupled oxygen transfer problem. The compact problem formulation yields rapid and stable convergence. Moreover, boundary effects can effectively be suppressed by generating very large replica of the cortical microcirculation in silico using an image-based cerebrovascular network synthesis algorithm, so that boundaries of the perfusion simulations are far removed from the regions of interest. Massive simulations over sizeable portions of the cortex with feature resolution down to the micron scale become tractable with even modest computer resources. The feasibility and accuracy of the novel method is demonstrated and validated with in vivo oxygen perfusion data in cohorts of young and aged mice. Our oxygen exchange simulations quantify steep gradients near penetrating blood vessels and point towards pathological changes that might cause neurodegeneration in aged brains. This research aims to explain mechanistic interactions between anatomical structures and how they might change in diseases or with age. Rigorous quantification of age-related changes is of significant interest because it might aide in the search for imaging biomarkers for dementia and Alzheimer’s disease.

Introduction Late-onset neurodegenerative diseases are believed to interrupt or degrade the coordination between vascular blood flow and neural tissue oxygenation. For example, hypoxic events have been implicated in beta-amyloid production in Alzheimer's Disease (AD) [1,2]. Patients with various forms of dementia also experience diminished oxygen tension in addition to significant morphological changes in cerebrovascular angioarchitecture [3][4][5][6][7][8]. Unfortunately, normal aging (aging without dementia) also entails similar morphological changes [3], so that normal age related changes are hard to distinguish from diseases. In order to better diagnose and eventually treat age-related dementia, it is imperative to be able to precisely quantify changes in vascular topology that potentially impair adequate oxygen supply to the brain. While the link between vascular restructuring and reduced tissue oxygenation have been established [4,[9][10][11][12], the quantitative contribution of hemodynamic factors to abnormal oxygen distribution have not been characterized due to the difficulty in posing and solving problems on a massive scale.
Multimodal medical images techniques enable the acquisition of detailed measurements about the microcirculation such as hematocrit, blood flow velocities, and oxygen concentration in the vasculature or the extravascular space [1,2,4,9,[11][12][13]. However, these invaluable data are typically collected from different specimen at different locations or time points with widely varying resolutions pertaining to a host of imaging modalities (μCT, 2PLSM, MRI). To direct the systematic search for imaging biomarkers of neurodegenerative diseases, it would be ideal to be able to focus data and observations across different length and time scales from diverse imaging modalities into a single coherent picture of hemodynamics and microperfusion in normal and pathological states. With this objective in mind, several groups have pursued computer simulations to quantify mechanistic interactions between cerebral blood flow and oxygen exchange supported by multimodal imaging data at the microlevel [14][15][16][17][18][19][20][21][22][23][24][25]. Especially, anatomically accurate models which recreate the detailed cerebral micro-circuitry bring the advantage that measurements at all length scales or acquired with separate imaging modalities can be combined in silico to make quantitative predictions about the relationship between blood flow and metabolism. Unfortunately, the solution of anatomically realistic mathematical models for cerebral circulation and oxygen exchange is a massive computational task.
Analytical and semi-analytical solutions of oxygen exchange. Analytical approaches such as the Green function method [31,32] approximate tissue oxygen tension by evaluation of an infinite series for oxygen supply from point or line sources representing blood vessels into a threedimensional tissue domain. Most analytic and semi-analytic singularity removal methods require summations of analytical expressions for concentrations and fluxes from single point or line sources into the domain, while satisfying homogenous boundary conditions. While such formulations are available for simple geometries (e.g. sphere, cube, cylinder), it is not obvious how to generate analytical formulations delineating complex biological spaces such as the highly gyrated cerebral cortex.
Methods using unstructured finite element meshes. Numerous groups have combined finite element discretization with 1D network graphs to create so called 1D-3D coupling approaches [26][27][28]30,36]. These elegant formulations will be reviewed in more detail in the discussion but a few selections are introduced here.
In a recent example, an FEM method was used to simulate extracellular species transport was in high resolution mesh measuring~4x4x4μm 3 with~82-84 million tetrahedrons [37]. As mentioned, FEM body-fitted meshes require contiguous representation of the space covering vascular segments and extravascular space. Even though excellent meshing tools are available such as ANSYS ICEM (Canonsburg PA), HyperMesh (Troy MI) and GMSH [38], segmentation of the cerebral microdomain for generating body fitted meshes precisely delineating capillaries, endothelium and parenchyma is a limiting factor in our experience. This bottleneck exists because automatic meshing often fails in the tortuous and highly bifurcated capillary bed, thus requiring manual corrections which is impractical due to the large number of segments. More importantly, tetrahedral segmentation of mixed vascular and extravascular domains entails a prohibitively large number of mesh elements needed to resolve the interface. Apart from the mesh size, it should be noted that the convergence of extremely large meshes suffers from a resolution-dependent deterioration of the matrix condition number, as shown by Briggs [39] even for the 1D diffusion. In effect, finely meshed domains necessarily result in an escalating iteration count, thus leading in practice to cases that do not converge at all or stagnate at poor approximations that does not satisfy all equations (= insufficiently small residual errors).
Another approach is to use dual mesh techniques where a 1D vascular network is coupled to a 3D extravascular mesh. These approaches benefit from problem reduction which improves solvability and avoids solving the 3D nonlinear Navier Stokes equations to predict vascular blood flow. Gagnon and Boas solved oxygen exchange between the microvasculature embedded in brain tissue using a tetrahedral mesh with sharp interface between the blood vessels and tissue [14,28,30]. Blood flow was computed separately using a simplified Hagen-Poiseuille network model; then oxygen extraction to tissue was solved by projecting the segment oxygen tension to the refined triangular surface mesh of each vessel segment. Other groups avoided the body-fitted meshing by registering a vascular network to a tetrahedral mesh and distributing the transmembrane flux for each segment across many tetrahedral mesh elements [33,34,[40][41][42][43]. Unfortunately, reported problem sizes merely encompass up to a few thousand vascular segments, because the required mesh sizes for fitting a contiguous extravascular mesh with the vascular network graph increase dramatically especially at the microscale.
A similar dual mesh technique has been proposed that registered the vascular network graph to a surrounding tissue domain represented by a coarse tetrahedral mesh [26,27,36]. Mass exchange was approximated as the total flux from the entire cylindrical segment to the mesh cell containing its center. This method already enjoyed the benefit that it did not require a body fitted mesh, thus enabling oxygen simulations for sizable portions of the cerebral cortex; 3x3x3mm 3 in humans [36]; and 1x1x1mm 3 in mouse [27]. Furthermore, the overall CPU time could be reduced by solving oxygen tension over the two meshes simultaneously, instead of having to iteratively solve oxygen fields in the blood and tissue domains separately. Unfortunately, point coupling between the centerline cylinder with surrounding tissue cell does not precisely orient and distribute the transmembrane flux across the endothelial interface, although the issue could be partially resolved using anisotropic and/or adaptive meshes.
Dual-meshes were also used to solve dynamic flow problems solving the Navier Stokes equations for fluid flowing within deformable domain boundaries [44,45]. A technique floating a moving (= Lagrangian) mesh inside a fixed (= Eulerian) mesh frame avoided remeshing in every time step, which would otherwise be necessary when simulating distensible blood vessels or a suspension of deformable red blood cells.
It is well known that parametric, structured meshing can drastically accelerate the convergence of transport problems. Several authors have developed homogenization methods to overcome the need for massive tetrahedral meshes by replacing body-fitted meshing with Cartesian meshes and anisotropic diffusion tensors [46,47]. Unfortunately, these methods do not capture discrete paths the blood takes through the capillary bed, so it is not yet obvious how homogenized model results can be interpreted or validated. These methods are discussed in more detail in the discussion section. Other groups performed dynamic 3D hemodynamic simulations effectively solving blood flow in the visible portion of the main cerebral arteries using parametric structured body fitted meshes [48][49][50]. The Sarntinoranont group introduced an elegant voxel-based method for predicting drug dispersion in spinal tissue; her method had the advantage that it seamlessly integrated the information flow between image data and simulations by aligning simulation results with image data at the same Cartesian "voxel" grid [51]. Here, we introduce a novel voxel based meshing-less discretization method that uses multiscale, parametric, structured grid techniques that drastically reduces problem size and accelerate convergence of oxygen transport simulations in the cerebral microcirculation.
Here we propose a method that rests on a structured Cartesian grid encompassing both the vascular and extravascular domain with a crisp separation interface obtained by a simple masking procedure. Masking of the interface only requires centerline a diameter information of the vascular network graph thus avoiding the need for body-fitted mesh generation altogether. This approach offers several critical advantages: (i) Masking operations can be executed using simple Euclidean geometry operators, eliminating the need for cell connectivity indexing (adjacency information). (ii) Memory storage for the Cartesian grid is minimal, since it requires only a few scalar parameters (domain bounding box and element counters for the main dimensions in x, y, and z-directions) which substantially reduces RAM requirements. (iii) Cortical oxygen perfusion problems discretized with our voxelized, mesh generation-less method enjoys excellent mathematical solvability. (iv) The cell masking paradigm is readily expandable to more specific anatomical labeling (smooth muscle domain, nonuniform endothelial thickness, or perivascular spaces). Finally, (v) the discrete labeling algorithm enforces that each vascular segment exchanges mass fluxes throughout its entire endothelial surface layer with the surrounding neighborhood of extravascular cells, which improves the computational fidelity of the predictions compared to point source 1D-3D coupling. This method is capable of computing oxygen gradients around vessels in the microvascular network at all relevant length scales. This new methodology converges simulations of significantly larger cortical microcirculatory networks than previously attempted. Moreover, the system is stable up tõ 100M mesh elements and~1M vascular segments on a desktop workstation. The substantial performance improvement puts us a step closer to whole-brain simulations.

Structural Cartesian mesh properties
First we show that the proposed structured, cuboid Cartesian mesh discretization significantly improves the matrix structure of the transport equations compared to traditional tetrahedral meshing as exemplified by Fig 1. The matrix representation of finite volume or finite diffusion fluxes using tetrahedral mesh has a very wide bandwidth, because many elements (= rows) The tetrahedral mesh matrix has a large bandwidth with no apparent block structure. In FEV discretization of unstructured domains, unknown nodal values typically occur in dozens of balance equations which augments matrix bandwidth thus substantially hampering the solvability of the underlying transport problem. (B) The structured Cartesian (cuboid) mesh covering the same domain with similar characteristic edge length has a much narrower bandwidth and has block structure which is ideally suited for block diagonalization preconditioning. Tetrahedral meshes also require almost ten times (in our example 9.2 times) more grid cells (nz) than the analogue cuboid mesh. have more than several dozen entries, each one indicating exchange fluxes with other cells. On the other hand, the Cartesian mesh has a highly organized block diagonal structure with a constant element-to-element connectivity of seven. Thus, the bandwidth of the mass exchange/ transport matrix stemming from the Cartesian mesh is much smaller than in the tetrahedral counterpart. Moreover, the transport coefficients in the matrix have typically the same value. Due to its superior structure, connectivity and uniformity of entries, the Cartesian mesh is more amenable to block preconditioning which enjoys more rapid and robust convergence.
Moreover, the number of elements required to represent a tissue domain with a tetrahedral mesh is significantly larger than for the Cartesian mesh. A hexahedral mesh encompassing the same space with the same characteristic edge length is, by default only~1/8 the size of a tetrahedral mesh. Specifically, expressing the cortical section into tetrahedral elements using ANSYS ICEM [Canonsburg, PA] resulted in a ratio of tetrahedral volumes to hexahedral volumes of 9.2:1 with similar characteristic edge lengths. Accordingly, parametric structured Cartesian meshes yield a problem reduction of almost one order of magnitude.
Voxelized vascular network representation-masking. As described above, we also use the Cartesian mesh domain to precisely delineate the cortical microcirculation. This choice gives the critical advantage that it eliminate the need for automatic body-fitted mesh generation, which often fails on microcirculatory meshes or requires time consuming manual adjustments to ensure adequate mesh quality at the interface between vascular luman and extravascular space especially at the microvessel length scale of~1μ m. the simple masking logic described in section 3.4 was applied to partition the "voxels" in the three dimensional domain into one of three groups; (i) interior vascular cells (red label), (ii) endothelial boundary cell (grey label), and (iii) tissue element (blue label for extravascular space). The structure of the resulting computational domain resembles a 3D image obtained from an imaging modality where the "grid resolution" determines the smallest features in the simulation. The edge length of the finest cuboid cell also delineates the minimum feature scale without the need for creating a body-fitted tetrahedral mesh. In effect, all advantages of the structural meshing described in the previous section carry over for the integrated vascular and tissue domains.
Results of the beneficial mesh generation of the voxelized vascular network are described next. The implementation of masking logic followed section 3.3. with more details in Fig 1 and S1 Text. The outcomes of the masking logic is illustrated using a microcirculatory network (~1x1x1mm 3 ) of the somatosensory cortex in one of the mouse Kleinfeld dataset [27,52,53]. The labeling results for the vascular network, endothelial interface and extravascular space are depicted in Fig 2. The interior lumen of the pial arteries, penetrating arterioles, capillaries, venules and veins is labeled in red. Large caliber vessels show blocks of red elements delineated by a distinct boundary layer (white elements) as shown in Fig 2C. Thinner vessels depicted in Fig 2D and 2E may have an endothelial layer measuring only one element (= isolated grey cells). Outside the endothelial layer, blue cells constitute the extravascular tissue space (= extravascular elements). The large number of tissue cells also underscores that the bulk of the computational effort is spent on solving oxygen tension in the tissue, while the number of equations to compute flow and oxygen in blood is typically about two orders of magnitude smaller. It should also be noted that the proposed Cartesian network masking preserves the block-diagonal matrix structure discussed above as can be seen in Fig 1, Right.

Divergence free and mesh independence simulations
We tested the ability of the Cartesian "meshing-free" method to consistently resolve microgradients around microvascular blood vessels in two large specimen of the somatosensory cortex of mouse. To this end, we tracked 25 rays through the cortical sample domain as a function of characteristic edge length. In these simulations, the convective flow field was solved a-priori using the matrix Eqs in (1) and (2). The oxygen fields in the vascular network and in the tissue mesh were converged simultaneously. For the oxygen dispersion in the vasculature, the convective flux uses the bulk blood flow and hematocrit; the mass transfer from the segment to the tissue is computed by summing the contributions of all endothelial boundary elements surrounding this segment. In effect, the oxygen field computation for the vasculature is executed on a tubular network but coupled to the surrounding tissue using the fuzzy representation of the boundary endothelial layer. With the convection equation for the vasculature and the mass transfer coupling equations for the interface and the diffusive transport and reaction equation for the tissue, the resulting vascular-tissue system in (3) is solved simultaneously using a preconditioned GMRES [54,55].
Divergence-free simulation. We calculated the overall mass balance error for blood flow simulation as well as for oxygen tension. Table 1 lists the results showing almost perfect closure A) The cortical network (E1.1) is labeled by pressure visualizing arteries (red), veins (blue) and connecting capillaries. B-C) This cortical network has been labeled a mesh visualized at many resolutions. Coloration of mesh is tri-color; red = blood, grey = endothelial, and blue = extravascular cells. Note, most of the cells belong to the extravascular space (blue). The blood lumen is marked by intravascular elements (red). Mass transfer occurs at the interface elements (grey). D-F) Visualization of mesh labeling at three resolutions reflects the capability of the labeling to resolve the connectivity between vasculature and tissue at all length scales. F) The distribution of endothelial layer, intravascular, and extravascular elements can be seen at the finest resolution.
https://doi.org/10.1371/journal.pcbi.1008584.g002 of both bulk blood flow and oxygen balances for the test cases. We also validated divergence free simulation of the oxygen field using spot checks in all locations of the vascular and tissue domains with a maximum flow balance error of |e| 1 <1.8�10 −10 %. Because the mass balance fluxes and reactive source terms were enforced using the finite volume approach, divergencefree computations were achieved both from the theoretical as well as the numerical perspective. More details in the finite volume discretization are given in S3 Text. Raytracing. Concentration profiles along rays in Fig 3 typically show sharp peaks near penetrating arterioles and valleys in locations further away and around venules. The computed oxygen concentration profiles along each ray for the first specimen (E4.1) further converged to a stable function even in the presence of steep gradients around blood vessels with high oxygen tension (= penetrating arterioles). Typically, mesh independence is achieved for meshes with 205 mesh elements per side, or a grid resolution of~5 microns. We also show successful convergence for more refined mesh sizes down to~3 microns, which was the limit in problem size imposed by available memory in our chosen inexpensive PC hardware (Intel Xeon W-2025, 1 core, 256GB RAM). A second large scale specimen (E1.1) exhibits similar satisfactory mesh independence with results shown in S2 Text. The oxygen simulation typically converged within 0.001-2.58 CPU hours (depending on mesh resolution). This investigation using two large scale cortical microdomains demonstrate excellent stability and mesh independence.

Precise quantification of oxygen microgradients in the cerebral cortex
We next wished to further analyze the pattern of oxygen dispersion in the extravascular space surrounding the cortical microcirculation. We performed simulations of blood flow and oxygen tension on large sections of the murine somatosensory cortex. To better highlight oxygen profiles around the microcirculation, we created contour maps and pO2 reliefs as shown in Fig 4 for the first specimen. Sharp peaks of high pO2 tension can be seen near oxygen rich penetrating arterioles. Note, that the high-resolution oxygen predictions inside the 3D computational domain also delineate a narrow plateau corresponding to the diameter of the vascular lumen.
Venules are typically located in flat basins of low oxygen tension. Accordingly, areas further away from arterioles but close to venules mark regions with relatively low oxygen tension depicted in blue. Note that the pO2 in venules is typically still above the tissue oxygenation level, which is lowest because of oxygen metabolism takings place inside the cell mitrochondria. We also could confirm that in some locations, venules may reabsorb oxygen from the surrounding tissue. However, our simulations seem to suggest that oxygen drainage is not a typical role for venules. The results in Fig 4 demonstrate successful simulations for a second cortical dataset and with trends confirming the consistent predictions seen in the first case study. Both case studies also quantified the evolution of the oxygen tension as a function of the cortical depth. For example, in Fig 4D a narrow band of low oxygen tension can be seen at a depth of 880μm below the pial surface. The concentration profile in the vicinity of major vessels exhibit steep oxygen gradients as expected. The precise quantification and detection of low oxygen pockets may be significant for the analysis of pathologies associated with aging.
These simulations show fine-grained resolution of the microgradients surrounding larger vessels while maintaining numerical stability at nPoints>100K vascular nodes, and nVols>8M tissue elements within only two CPU hours on a personal computer (Intel Xeon W-2025, 1 core, 256GB RAM). High resolution of the oxygen field on the micron scale constitutes the desired simulation objective.

Validation of tissue oxygen predictions against maps from 2 photon imaging
We further tested whether predicted tissue oxygenation fields would match actual in-vivo 2-photon phosphorescence lifetime imaging oxygen measurements. Accordingly, we acquired oxygen tension maps in three mouse specimen at various cortical depths (up to 250μm below the cortical surface) with the acquisition protocol given in the Materials and Methods section. We reconstructed the vascular microarchitecture in these specimen and simulated tissue oxygen tension in several plains below the cortical surface. In all three cases, predictions of oxygen tension utilizing this mesh masking technique shows qualitative and quantitative agreement with the experimental distributions as can be seen in Fig 5. For comparison, we first generated The experimental results were obtained with 2-photon phosphorescence lifetime oxygen imaging. The simulations were performed in planes with two different resolutions; a coarse 20x20 and a dense 200x200. There is reasonable agreement between experimental and simulated oxygen fields (2.9±1.4% error). Manual matching between the concentration profiles was used to determine tissue diffusivity, reaction rate, and mass transfer coefficient. The 3D perspective illustrates the steep concentration gradients along two penetrating arterioles.
We also produced high-resolution simulations (200x200 resolution = 400x400μm 2 ) which resolve the oxygen gradient at a much higher resolution than the data in the experiments. The high-resolution simulations also allow us to determine physiological parameters directly from measured oxygen data. Specifically, it is possible to vary diffusivity, D, mass transfer coefficient, U, and even oxygen metabolic rate in the simulation until there is optimal alignment with the experimental oxygen field. This procedure can be interpreted as a multi-dimensional least-square adjustment of unknown physiological parameters directly from experimental data. For the purpose of this demonstration, we performed the adjustment manually with resulting values listed in Table 2. Here, the measured oxygen concentration derived from the original image (= at edges of the field of view and in the center of the blood vessel) served as boundary conditions for the predicted oxygen field. In principle, one could formulate a multiparametric multi-dimension parameter estimation problem. A mathematical programming solution for the parameter estimation problem is unfortunately outside the scope of this paper.

Applications
Two larger scale applications are discussed in this subsection to demonstrate the significance of the proposed technique for assessing age-related metabolic decline and neurodegeneration.
Hypoxia in aged brains. The first application concerns a preliminary exploration of agerelated changes in the oxygen perfusion of the cerebral cortex. The aged brain suffers from a variety of vascular changes including an increased incidence (~30%) of micro-occlusions [11], systemic hematocrit decrease by~30% [4], and a~20% decrease microvascular vessel density [3][4][5][6][7][8]. It is further believed that these changes may lead to hypoxic pockets, which are regions below pO2<10 mmHg that may suffer damage to neural and glia cell lines. We wanted to test the hypothesis of the effect of aging by simulating oxygen tension in vascular networks, where numerous known age-related morphological or hemodynamic changes were imposed artificially in silico. The simulated reduction in oxygen tension was assumed to characterize possible deterioration due to different age-related pathologies.
Quantitative results of the preliminary study on the effects of aging are summarized in Fig 6  which shows comparisons of oxygen perfusion patterns in the normal cortex (young brain, row 1) with different aged brain scenarios (rows 2, 3, and 4). The second row of Fig 6 characterizes the effect of age-related occlusions, which were induced in the simulation by setting the diameter of 30% of the capillaries to a value of numerical zero. Despite the substantial reduction in capillary perfusion, the oxygen tension stayed above the hypoxic threshold for the  The observed lower tissue oxygenation in row 3 suggests that systemic hematocrit reduction has a stronger effect on lowering oxygen tension than micro-occlusions.
Finally, the combined effect of reduced hematocrit with increased occurrence of microocclusions is quantified in the oxygen maps of row 4 of Fig 6. In this combined effect scenario, hypoxic pockets with oxygen tension under a critical level of pO 2 <10mmHg begin to form especially in the deeper cortical layers (880 microns). The simulations results support the possibility of hypoxic pocket formation due to the combined effects of age-related metabolic decline in the angioarchitecture and hemodynamic functions.
Our simulations demonstrate that high-fidelity computations with the Cartesian meshing are promising for characterizing age-related changes to tissue oxygenation. It should be noted that these results are preliminary and a detailed investigation of age-related changes to cerebral oxygen tension is beyond the scope of this study.
Large-scale cortical network simulations to eliminate boundary effects. We observed that the depletion of vascular segments in the boundary zone is a main contributor to boundary effects in oxygen simulations. To overcome this undesired boundary sensitivity, one could artificially seed border regions with artificially added blood vessels so that the vascular density is not disturbed at the domain edge. However, it is not clear which metrics should inform such a "correction" for irregular edge vascularization. A more concise strategy entails generating large samples of the cortical microcirculation so that the domain boundaries are far away from the area of analysis. This idea is demonstrated next. We created a large artificial vascular network by using an image-based circulatory network synthesis (iCNS) algorithm which is described in detail elsewhere [52,59]. Specifically, we synthesized a topologically equivalent murine microcirculatory network sample covering 9mm 2 of its cortical surface to a depth of 1.2mm. The resulting microcirculatory network model depicted in Fig 7 is roughly one order of magnitude larger than the largest 2-photon Kleinfeld data set, and much larger than samples of previous studies. We limited our analysis to the inner region of interest (ROI) which fits into a box of 1x1x0.8mm 3 and is perfused by an uninterrupted, physiologically-consistent microcirculatory network with domain borders so far removed that they do not affect the oxygen simulation in the ROI. We conducted blood flow and oxygen simulations over this massive network with nPts~1M network segments and nVols~100M mesh elements. The Cartesian masking method converges stably for this massive domain sizes in~28 CPU hours on a personal PC (Intel Xeon W-2125 processor at 4.01GHz and 256 GB of Micron MTA36ASF4G72PZ-2G3B1 memory). Fig 7E shows the mesh domain, but also the unavoidable boundary effects at corner, bottom, and side edges of the domain. The calculated oxygen fields in the ROI showed no observable boundary effects.
For this massive cortical sample, we also simulated the effect of aging by the combined effect of reducing the inlet hematocrit and imposing micro-occlusions in 30% of the capillaries as described above. Inside the ROI, we found that young normal brains exhibited no areas with critically low oxygen tension (= young brain, Fig 7C). In contrast, multiple hypoxic pockets with oxygen concentration below the critical pO2<10mmHg threshold formed in the aged brain. An example of the spatial extent of a hypoxic pocket is depicted in Fig 7D. We also computed the volume fraction occupied by mild hypoxic regions (pO2<10mmHg) and tabulated the increase in occurrence for cohorts of normal (N = 7) to aged (N = 7) murine brains in Fig 7B. The simulated values were compared to detailed experimental results obtained in a different study using two-photon phosphorescence lifetime imaging [13], which studies oxygen tension in cohorts of young, middle aged and old mice. The preliminary comparative data show that the simulations matches the incidence of low oxygen regions in the normal brain. Moreover, the incidence of severe hypoxic pockets (= regions with oxygen tension less than 5mmHg) observed in the different age cohorts is closely matched in the massive network simulations presented here. The occurrence of regions in aged brains where oxygen tension was less than 10 and 15mmHg respectively were drastically increased in the aged population, both in the experiment as well as in our large-scale computer simulations.

Discussion
Voxelized simulation domain. The proposed "voxelized" Cartesian meshing-free methodology offers a robust approach for stable convergence of blood flow and oxygen exchange simulations for very large, anatomically detailed sections of the cerebral cortex. This method succeeded in predicting blood and oxygen micro-perfusion in massive cortical microcirculatory networks with stable convergence up to relatively dense mesh resolutions in the micron range.
The use of a voxelized grid with cubic elements (= cuboids) delivered two critical benefits: First, a mask based logic introduces a structured Cartesian grid representation which resolves the interface between the endothelium and the extravascular space without the need of manual body-fitted segmentation required for tetrahedral grids. Therefore, we termed the method mesh generation-free, which means that it does not require the creation of a body fitted mesh. Second, the stencil of the resulting Cartesian mesh always contains exactly seven cells, so that the bandwidth of the coupled mass transfer problem is always much smaller than in tetrahedral grids. Moreover, since the Cartesian meshes are structured and parametric in the main coordinate directions (x, y, z-directions) leading to a mesh comprised entirely of uniform Comparison between oxygen tension in young mouse and aged mouse suffering from reduced hematocrit in a large size synthetic vascular network spanning 3x3x1.2mm 3 of the mouse cortex. A) The vascular structure colored by diameter (red = thick, blue = thin) shows a dense set of tortuous vessels connected through a complete capillary bed. B) The tabulated occurrence rate of low-oxygen mesh elements in a portion (= a 1x1x0.8mm 3 tissue subset from the interior of the domain) of the predicted oxygen field shows the same trends as in the experimental data. C) Predicted pressure distribution throughout the network matches the topology of arteries, veins, and capillary closure. D) A cut through the dense structure reveals regions of hypoxia and pockets of hypoxia occur more frequently in the lower regions of the cortex. D) A zoomed view of oxygen tension in the interior of the vascular network. There are some segments with low oxygen content (but above the hypoxic threshold) especially in deeper regions). E) The oxygen levels in the tissue reveal the lower regions are more susceptible to lower oxygen tension then regions near the cortical surface. A detailed view compares oxygen tension in young to aged specimen. In the young specimen, no hypoxic areas form. In the simulated age mouse, hypoxic regions and their spatial extent were identified; hypoxic pockets are indicated in red and are larger and more prevalent in the aged brain.
https://doi.org/10.1371/journal.pcbi.1008584.g007 PLOS COMPUTATIONAL BIOLOGY cuboid elements, there is no need to store the tissue domain grid, in effect drastically lowering the memory requirements of the simulation code. That property also justifies the naming of a technique that is mesh-free, because it does not require the cumbersome generation of memory expensive fitted meshes. The compact structured cubic mesh with narrower matrix bandwidth reduces problem size, stabilizes the perfusion problem against inaccurate initial guesses and drastically accelerates convergence. The drastic improvements in the structural matrix properties of the transport problem were illustrated in Fig 1. The proposed binary masking technique also enables automatic assignment of mechanistic expressions for oxygen transfer fluxes between vascular segments and the adjacent tissue elements through a fuzzy endothelial cell layer. The sharply defined interface obtained with this technique does not introduce singularities in the mass transfer flux computations. In other numerical techniques, singularities may occur when a massive flux discharges into a small volume cell that cannot accommodate the flux without incurring an extreme concentration drop. This does not occur in the current discretization technique, because large mass exchange fluxes emanating from thicker vessels are partitioned and distributed across the endothelial interface in proportion to the tissue cell edge length. In other words, the massive oxygen efflux is divided among the set of adjacent cuboids in agreement with the characteristic edge length thus avoiding singular concentration gradients.
Moreover, proposed "voxelized" simulation domain facilitates the direct comparison to stacked slices of neuroimage data. Thus, microcirculation simulations can be conveniently set up directly from image stacks without the need of first segmenting the vasculature-tissue interface with a body-fitted tetrahedral grids, a task that is currently infeasible for automatic mesh generation software.
Boundary effects. Lorthois [15,60] pointed out that boundary conditions severely affect microcirculatory hemodynamic predictions. In their microcircular blood flow simulations, Schmid et. al [61] assembled an array of cloned samples to limit border effects. Gould et al [27] cyclically extended the tissue space with vascular segment connections across opposing domain borders to diminish boundary effects. In the current study, we observed that boundary effects in oxygen tension computations mainly stem from the fact that neuroimaging data typically do not capture vascular topology accurately at the edge of the imaging domain. This irregularity especially affects branches of large penetrating arteriole trees (causing dangling segments that often have to be removed). The removal causes the boundary region in simulations to appear undersupplied with oxygen carrying blood vessels, thus creating a marked disturbance of oxygen perfusion in the border zone. To overcome this undesirable disruption, one could populate the boundary region with artificial blood vessels so as to compensate for the severed segments at the domain edge. However, it is not clear which metrics should inform the "correction" for irregular edge vascularization. Our more concise solution strategy entails generating large samples of the cortical microcirculation so that the domain borders are distant from the area of interest.
We applied an image-based cerebrovascular network synthesis (iCNS) algorithm [59] to create a replica of a sizable portion of the murine somatosensory cortex with dimensions of 3x3x1.2mm 3 . Note that these dimensions are much larger than previously reported microcirculatory oxygen simulation domains. Then we limited our analysis to a region of interest (ROI) at the interior core. The ROI measuring only 1x1x0.8mm 3 has all its sides far removed from domain boundaries. Edge effects in oxygen supply due to the disruption of the integrity of the vascular bed near the boundary were thus kept away from the predicted oxygen fields. However, expanding the computational domain to eliminate boundary effects necessarily requires the large-scale solution of oxygen perfusion problems which is fortunately enabled by the proposed voxelized simulation technique.
Multiscale acceleration. Linear algebraic system equations for ultra-dense grids (= massive cell numbers) do not always converge from poor initial guesses. Moreover, the numbers of iterations to achieve convergence may be prohibitive due to the deterioration of the system condition number that comes about with fine mesh resolutions as was mentioned in the introduction and is discussed in detail elsewhere [39].
To ensure robust convergence for arbitrary initial values, we took advantage of multiscaling effects by obtaining approximate oxygen fields initially on a coarse mesh, then computing intermediate updates on Cartesian grids with finer resolution, before converging the final solutions at the finest scale. The masking procedure needed for meshes of different resolution can be performed easily by simply setting the edge length to a smaller value. Coarse-grain simulations rapidly converged due to their reduced problem size, but lacked details such as local oxygen micro-gradients. Computed coarse-grain oxygen fields were then linearly interpolated down to the finer partitions (= finite element interpolation) and used as starting conditions for fine grid simulations. In all cases, we succeeded in converging from uninformed, of "poor", initial guesses where frontal fine-grained attempts tend to diverge.
Our method differs from classical deflation as it does not solve dual-stage problems using the residual vectors to compute reduced-order updates, but merely uses consecutively more refine-grid solutions as better initial guesses. In its current form, our method does not benefit from such accelerated convergence methods, but does seem to have a stabilizing effect by suppressing high frequency updates and preventing them from overtaking the solution.
An interesting, unexplored approach using our mesh generation-free Cartesian technique would be to combine it with deflation as proposed by Nicolaides [62]. It can be expected that the geometric (= positions in 3D space) and the logical regularity of the Cartesian mesh (= block diagonal arrangements of matrix entries) tremendously aides in the implementation of geometric multigrid techniques [39,63]. The numerical accuracy of proposed computational technique might be further improved with the use of logarithmic profile assumptions [35]. We did not yet attempt to deploy geometric multigrid algorithms [64] with the structured meshingless perfusion problem formulations described in this manuscript, although these can be used in the future. Reported convergence could be further accelerated by algebraic multigrid solver algorithms [65] to further performance gains, but these ideas are outside the scope of this study.
Fuzzy 3D blood flow computation. We also tried to solve blood flow and oxygen convection over the "voxelized" 3D vascular grid at differing resolutions but found this to be unstable in our implementation [66]. Thus, we avoided numerical problems with frontal 3D blood flow computations by instead, merely solving the blood flow and pressure fields once over the original network (before segmentation). The fluxes and pressure fields were interpolated to assign the hemodynamic states for segment subdivisions.
Comparison to prior methods. To better assess the significance of the proposed technique, a selection of representative prior work of microcirculatory simulations are compared in Table 3. The list is not complete and many valuable contributions could not be incorporated due to space limitations; specifically we did not list homogenization techniques [46,47] because they do not currently address oxygen simulations, but approximate capillary blood perfusion using a continuum surrogate for the microvessels. The proposed "voxelized" method compared favorably in problem size as measured by the number of vascular segments and mesh elements in relation to previous techniques. The simulation time for the E4.1 data set (with a resolution of 305x305x305 elements) was~1 CPU hour and the simulation time for the 3x3x1 (with 455x455x455 mesh elements) was 19.4 CPU hours. Note, the current simulations of the microcirculatory blood flow and oxygen fields used an unprecedented fine grid resolution (= micron scale).
Analytic and semi-analytic approaches that were introduced in the Introduction section resolve the oxygen exchange fluxes, nested infinite series approximating the oxygen fields for the vascular network and the extravascular mesh domain have to be repeatedly evaluated until convergence is achieved. Accordingly, the computational complexity grows considerably with the number of source points (= vascular segments). So far, oxygen exchange for networks with a few hundred segments have been demonstrated [31,32,35]. The reported problem sizes should be seen in relation to the required vascular density of~11,600 vascular segments/mm 3 of tissue. The second semi-analytical method in this group is splitting. Gjerde et al [35] substantially reduced the computational effort by separating the solution into an analytic base component and a corrective field, which needs to be determined numerically. Although this offers considerable performance gains, these methods are mathematically intractable for networks at the scale we use here.
Other groups used a dual-mesh technique for modeling water exchange on bifurcating vascular models up to a few dozen segments with variable permeability to simulate the difference between healthy and diseased capillaries [67]. To avoid singularities from the 1D-3D coupling, the authors distributed extravascular discharge flux using a nonuniform shape function, which was able to stabilize solvability even at coarse grid resolutions [68]. Shape functions are also used in algorithms to discharge the vascular flux into the homogenized domain [46]. Other groups developed sophisticated methods that coregistered the circumference of the tangential disc representing each cylinder with the surrounding tetrahedral mesh [33,34,[40][41][42][43] to avoid singularity issues. High quality dense meshing is required to effectively resolve gradients surrounding the smallest capillaries, so that adaptive meshing becomes necessary [33]. Despite their mathematical elegance, simulations for sizeable sections of the cortical microcirculation have so far eluded these techniques due to their fine micro-resolution requirement.
It is critical to underscore the benefits of the computational approach with simulations on a domain sizes comparable to in vivo image data. The results here demonstrate several key innovations: i. The ability to solve flow and oxygen exchange over substantial domain sizes enables validation of predicted oxygen gradients and fields against in vivo multiphoton image data. While several prior approaches laid out elegant mathematical foundations for highly accurate 1D- 3D coupling computations, these approaches engender intractable problem sizes when bridging the micro to macroscales characteristic of perfusion in the cerebral circulation.
ii. Successful convergence of massive problems over domains that are large enough that their boundaries are far removed from the region of interest. Our simulations further show that oxygen perfusion simulations over domains smaller than 1x1x1mm 3 are severely affected by boundary choices, so that predicted results may be excessively affect the simulated metabolic conditions instead of physiological oxygen exchange. Severe boundary effects are depicted in the oxygen field simulations shown in S4 Text.
iii. The simulations were highly sensitive allowing the prediction of subtle perfusion changes related to age-related decline. Specifically, we are able to confirm the experimental discovery of hypoxic pockets acquired with in vivo multiphoton observations, and predict with fidelity size and location of hypoxic tissue locations.
Limitations. The preliminary results reported in this manuscript on the effects of aging should not be interpreted as in-depth physiologically validated conclusions about the aging brain. The study merely aimed at showing that the proposed method is applicable to problem size with fine-grid resolution needed to quantify even subtle effects of age-related hemodynamic changes on tissue oxygen tension. A separate study to investigate the mechanistic sources of age-related oxygen tension changes is in preparation.
Furthermore, we would like to point out that the change in oxygen extraction also impacts the tissue pH which further feeds back on the kinetics of the oxygen dissociation from hemoglobin to blood plasma [69,70]. Without accounting for pH effects, the predicted oxygen tension may be inconsistent with the complete physiological response seen in vivo. However, we are currently not aware of any methodology that would simultaneously solve for blood flow, hematocrit, oxygenation, dissociation kinetics, oxygen and CO 2 metabolism and its effect on the pH Accordingly, we defer the inclusion of pH changes to future work.
We also note that the assumed pressure drop across the microcirculatory network of Δp = 115mmHg is somewhat large. We have pointed out previously [27,52,58] that our simulations of the Kleinfeld dataset use the original diameter specification which may report some diameters in the smallest capillaries too narrow so that larger pressure driving forces are needed to obtain physiological flow rates.
Conclusions. In conclusion, the proposed voxelized technique which requires no mesh generation offers a scalable method for performing computer simulations of solute exchange and metabolism for sizable sections of the cerebral cortex. The numerical techniques presented here compliment recent advancements in microcirculatory network synthesis and simulation [52,59]. A combined approach that integrates neuroimage data from multiple sources into mechanistic perfusion simulations of the mouse cortex was illustrated. Several anatomically detailed applications showed stable mesh convergence up to~100 million equations. It should be easy to further expand this limit, which in our case studies was limited only due to the memory capacity of our local personal computers. Predictions derived from anatomically detailed mechanistic models are expected to be useful for testing medical hypotheses for normal or pathological states.

Ethics statement
All protocols executed at the University of California, San Diego were approved by the Institutional Animal Care and Use Committee at University of California, San Diego. Animal handling and surgical procedures performed at the Montreal Heart Institute were approved by the ethics committee of the research center of the Montreal Heart Institute. All experiments at the Montreal Heart Institute were performed in accordance with the Canadian Council on Animal Care recommendations.

Vascular structure acquisition
Four sections of the murine vibrissa primary sensory cortex [53] were reconstructed from two-photon laser scanning microscopy (2PLSM) images. This dataset is referred to in this paper as the Kleinfeld data set. These images represented the length and orientation of the blood vessels [53,71,72]. These structures were then labeled with an automated algorithm as described in the original manuscript. The categorization used size and branching level information (Strahler order) to differentiate between pial vessels, penetrators and capillaries. Capillaries were identified by a diameter cutoff of 6 μm and penetrating vessels differentiated from pial vessels with depth and diameter thresholds. The final reconstructed network topology and diameter information was stored using sparse connectivity matrices. More details on image acquisition [53,71,72], image reconstruction [73], as well as the formulation of the network equations [58] can be found elsewhere. All protocols were approved by the Institutional Animal Care and Use Committee at University of California, San Diego.

Oxygen tension measurements in young and aged brain
Data acquisition and methods are described in the original experimental study reported in [13]. In short, the O2-sensitive molecule PtP-C343 [74] was used to measure tissue oxygen in awake animals using a cranial window and head fixation. Following fixation, two-photon imaging was performed using a custom-built laser-scanning microscope that used a consecutive sequence of 820 nm, 80 MHz, 150 fs pulses through an electro-optic modulator (ConOptics, USA) to adjust the gain and allow the generation of alternating on and off laser pulse periods for microsecond lifetime imaging required for oxygen quantification. Cortical tissue oxygen data was acquired in young, middle aged, and old mice, scanning 400um x 400um adjacent planes to estimate O2 in large areas. Twenty (20) mice were used for tissue oxygen measurements, 7 young (8.8±0.1 months), 6 middle-aged (15.3±0.1 months) and 7 old (27±0.1 months) at various depths up to 250-300μm below the cortical surface. Animal handling and surgical procedures were approved by the ethics committee of the research center of the Montreal Heart Institute. All experiments were performed in accordance with the Canadian Council on Animal Care recommendations.

Mathematical blood flow and oxygen model
Blood flow and hematocrit were computed first, then the oxygen fields in the vasculature and the tissue are solved with the blood and hematocrit fields as inputs.
Blood flow. Blood flow is modeled as a biphasic suspension with details in [26,52] and summarized in system (1)-(2), which has nonlinear coupling between the flow, f, pressure, p, and hematocrit fields, h. To solve the coupled flow problems efficiently, we perform fixed point iteration which computes the flow rate, f, using the diameter dependent resistance matrix, R, for an assumed hematocrit field, h. Approximate blood flows can then be computed with (1) using linear algebraic methods. With the converged flow field, f, the linear subsystem for hematocrit in (2) is solved consecutively to update the prior hematocrit field. Hematocrit updates are in-turn used to recalculate the flow field in (1) until the two fields converge. In our implementation convergence is typically achieved in less than 30 iterations.
Here, C 1 2 R nArcs x nPts is the fundamental incident matrix (= arc connectivity matrix) as expressed by Tellegen's theorem [75] and R 2 R nArcs x nArcs is the matrix of vessel segment resistances. Eq (2) is used to compute the hematocrit field with plasma skimming [52,58].
More details can be found in [26,52,59]. The matrix � M c 2 R nPts x nPts encodes hematocrit convection with � h given inlet boundary conditions enforced by the decision matrix D 0 2 R nVol x nVol . The counter nArcs indicates the number of vascular segments, and nPts is the number of vascular nodes.
Oxygen field. The blood flow, f, and hematocrit fields, h, computed from (1)-(2) serve as inputs to a convection-advection oxygen transport model with metabolic reactions in brain tissue. The final system is summarized in simplified conceptual matrix form in (3), which we solve simultaneously for oxygen tension in blood vessels, c v 2 R nPts , and in tissue, c t 2 R nVol . Oxygenated blood enters the cortical microcircuitry through pial arteries; oxygen travels by convection through the microcirculatory network and exits through the pial venous outlets, M c 2 R nPts x nPts . Oxygen also transfers across the blood brain barrier (BBB) into the tissue (= extravascular space), where it metabolizes to satisfy the energy demand of neurons and glial cells. The connectivity matrix C 3 2 R nPts x nPtsþnVol enforces the mass transfer between the vascular nodes and endothelial mesh elements which is scaled by the mass transfer conductivity, The matrix M d 2 R nVol x nVol encodes the diffusion in endothelial and extravascular (tissue) mesh elements. The diagonal matrix R 2 is used to implement the cerebral oxygen demand (CMRO 2 ) and contains the reaction rate constants, k met , and the cuboid volume (V t ).
In system (3), � c t and � c v are the known Dirichlet boundary conditions in the tissue and vasculature respectively aided by the decision matrices D 1 2 R nPts x nPts and D 2 2 R nVol x nVol . In the simplified notation of (3), we show the reduced system in which the interior vessel lumen cells have been eliminated (by identity with the vascular graph). More implementation details are offered in S1 Text. The coupled oxygen problem is further presented in component form for better description of the system equations as follows: Oxygen transport in the microvasculature. According to mass conservation, the change in convective oxygen flux (= the divergence of oxygen convection) carried by the blood stream must be equal to oxygen transfer to the tissue as in (4). As shown above, this relation can be conveniently expressed with the connectivity matrices with the added mass transfer from the vasculature into the tissue. Here, nVol, stands for the number of Cartesian tissue volumes (= cuboids). Moreover, U is the transmembrane permeability of the endothelial layer, w is the endothelial layer thickness, A is the endothelial layer surface area.
Oxygen transport in tissue. Oxygen transported to the tissue domain is permitted to further disperse by diffusion or undergo cellular reactions according to the cerebral oxygen metabolism in (5). The tissue oxygenation model is a diffusion-reaction system with coupled mass transfer sources from the vasculature. The parameters here are; D, the isotropic molecular diffusivity and, k, the first order cerebral metabolic rate of oxygen (CMRO 2 ). In extravascular mesh elements (labeled as blue in Fig 1) there is diffusion and reactions. In the endothelial cells (grey in Fig 1), there is in addition mass exchange with the vascular network. Note, we chose to neglect intravascular diffusion (= diffusion inside the blood vessel lumen), because the magnitude of diffusion is significantly smaller than advective fluxes (Peclet number >> 1) in the axial direction. Inclusion of a radial diffusion term would be possible, but was not undertaken here, because the required continuum assumption of radial diffusion in single phase does not match physiological discrete RBC interaction with the glycoclyx, so we felt justified to neglect radial diffusion at the microscale. Thus, radial concentration profiles were not resolved in our model.
Note that system (3) above is a simplified concatenated version of the component form of the vascular in (4) and the tissue equations in (5). We also note that the intravascular mesh elements (labeled as red hexahedrons in Fig 1) do not participate in mass transfer, diffusion, or reactions, but can be eliminated from the system by the equality with the associated vascular network node as in (6).
Where c i t is the concentration of the tissue element and c j v is the concentration of the nearest vascular node.
All physiological parameters and boundary conditions are similar to previously reported values [27,36,58] for oxygen in the brain as listed by Table 2. These values were estimated from experimental data and were found within the range of reported values.

Automatic equation generation and solution
A masking procedure with edge detection logic was applied to identify connectivity between the two domains. We first produce similar characteristic lengths between the Cartesian cuboids and vascular network segments by re-segmenting the vasculature as described in S1 Text. By re-segmenting the vascular network, we ensure the vascular segments have a similar characteristic length as the edges of the cuboid mesh. This stabilizes the simulation, ensuring the vascular segments are able to display realistic axial gradients and that the mass transfer between the vasculature and the surrounding tissue is sensitive to these gradients. An example of partitioning vasculature into many sub-segments is offered in Fig 8B-8D. The process then undergoes two important steps: (i) masking and (ii) equation generation.
Step-1. Masking. The simple fuzzy algorithm label3DMesh traverses every vascular segment to identify a bounding box for the centerline as described in S1 Text. To ensure full enclosure of the entire cylindrical segment, the bounding box is then expanded by the width of the segment radius. Once the set of boundary indices is determined, the indices of this subdomain of cells are temporarily stored, denoted as vesselNeighborhood. For each cell in the vesselNeighborhood, the Euclidean distance, d, to the vessel centerline is calculated. A cell is labeled as an edge cell, if its center lies inside the endothelial layer defined by a thickness, w/2, from the vessel edge. It is labeled as an interior cell, if it is not an edge volume and d<r. Otherwise, a cell is labeled an extravascular cell, which is also the default assignment. The pseudocode for masking the vesselNeighborhood is listed in S1 Text. An example of mesh labeling is offered in Fig 8A-8C.
Step-2. Equation generation. Once all cells of the Cartesian domain are labeled and connectivity matrices defined as described in S1 Text, we begin with generating the mass balance equations for each element in the simulation domain. First, we assembled the convection matrix for the vascular network (M c ) using an upwinding scheme. For its assembly, it is beneficial to reuse incidence matrix of the vascular network (C 1 ). Furthermore, M c , is asymmetric because of upwinding of the flow field. We did not account for oxygen convection by interstitial water fluxes in the extravascular domain, because such convective transport is not known to play a role in tissue oxygenation.
We then assembled the mass transfer connectivity matrix, (C 3 ) which connects the graph network nodes to the endothelial mesh cells. The connectivity information is obtained by the masking procedure as described in S1 Text and shown in Fig 2. When multiplying C 3 with the diagonal matrix G 1 , we obtain the mass transfer flux balances. The mass transfer flux balance is A) A penetrating arterial tree with the overlay of a Cartesian mesh reflects a network whose straight segments span many cuboid cell volumes. B) The vasculature is re-segmented (partitioned) until all segments span no more than two adjacent volumes. C) A labeled vessel inside a Cartesian mesh shows the extravascular elements in blue, intravascular segments in red, and endothelial segments in grey. D) The original centerline before and after partitioning for a penetrating arteriole sub-tree. In this illustrative case, the algorithm requires 8 steps to fully partition the tree.
https://doi.org/10.1371/journal.pcbi.1008584.g008 then accordingly C T 3 G 1 C 3 , forming a symmetric submatrix system. This mass transfer submatrix also shows structurally very nicely the coupling between the network graph and the Cartesian tissue mesh in our dual-mesh technique.
The diffusion in the tissue is enforced using the diffusion matrix M d . More details on how to assemble the diffusion matrix can be found in S1 Text. Note, the diffusion submatrix for diffusion is symmetric.
The oxygen consumption in the tissue is implemented with the help of the diagonal reactivity matrix R 2 . Oxygen consumption in the fourth submatrix is also symmetric.
Finally, on the right hand side, we show conceptually the implementation of boundary conditions for the vascular nodes, � c v , and for the tissue domain, � c t . Here the sparse decision matrices D 1 and D 2 formally indicate those nodes where known boundary values (Dirichlet boundary condition) need to be enforced. Neumann flux boundary conditions can easily be implemented with the logic indicated in (1), but has been omitted here to avoid detraction from the main structure of the overall coupled mass transfer system. Implementation note. In the formulation of (3) we have chosen not to distinguish the tissue from the endothelial nodes so as to avoid clutter in the equation. However, the connectivity matrix of the mass transfer correctly associates network nodes c v with endothelial nodes. Moreover, we have also omitted the equality between interior vascular nodes on the Cartesian grid and the associated network graph nodes (c v ) as in (6). We did not include the interior vascular lumen nodes because they can be eliminated from the system altogether, at which point (3) describes the reduced system. More details given in S1 Text.
All mesh cells in the tissue domain undergo metabolic reactions, R 2 matrix. Very small segments that are entirely embedded in a single extravascular cell exchange oxygen with this one cell only. This happened in the smallest segments of the capillary bed, depending on the Cartesian mesh edge length. In general, numerous cells delineate the endothelial layer of a blood vessel segment. According to the fuzzy mesh resolution, the total segment surface exchange area is evenly divided amongst the mass transfer mesh elements.
Solution algorithms and Implementation. In summary, the system to be solved is summarized in systems (1)-(2) for blood flow and hematocrit; and (3) for oxygen concentration in vasculature and tissue. Equation generation was written with object-oriented codes. All linear algebraic computations were performed in PETSc using the GMRES solver and a block Jacobi preconditioner [54,55]. We did not employ an algebraic multigrid solution algorithm.
Supporting information S1 Text. Implementation details. An in-depth description of how to programmatically implement the anatomical labeling and formulate the equations.