Computational feasibility of calculating the steady-state blood flow rate through the vasculature of the entire human body

The vascular geometry used in this study was constructed with a symmetric, fractal approach. Vessels were represented as rigid, cylindrical tubes that were connected at junctions to form a closed network. The term junction denotes any location where vessels meet, or the start and end points of the network, where boundary conditions are applied. The fractal model was conveniently scalable to billions of vessels and had only two boundary conditions, namely the blood pressure at the inlet and the blood flow rate at the outlet. The number of vessels (N_v) in a network was related to the number of bifurcations or generations, denoted by N_g, by

$$\begin{eqnarray}&&{N}_{{\rm{v}}}=2\left({2}^{{N}_{{\rm{g}}}}-1\right)={2}^{{N}_{{\rm{g}}}+1}-2.\end{eqnarray} \tag{ 1 }$$

In the model, each parent vessel branches into two child vessels. The smallest vessels in the network were 2.5 μm in radius and 55 μm in length (Lauwers et al 2008). Between each generation, the vessel radius was reduced by 2^−1/3 following Murray's Law for bifurcation (Adam 2011) and the length was reduced by 0.8 to reduce vessel overlap for visualization (Donahue and Newhauser 2020). Figure 1 plots an illustrative example of a 2-dimensional network created using this algorithm. The algorithm used to construct this network geometry was described previously (Donahue and Newhauser 2020). Each vessel was modeled with a start location, end location, and inner radius, as well as junction indexes, and a list of immediately connecting vessels.

Zoom In Zoom Out Reset image size

At a point ${\boldsymbol{x}}=\left({x}_{1},{x}_{2},{x}_{3}\right)$ in an arbitrary 3-D coordinate system, the time-dependent form of the Naiver-Stokes equations are

$$\begin{eqnarray}&&\rho \displaystyle \frac{{\rm{D}}{\boldsymbol{v}}\left({\boldsymbol{x}},t\right)}{{\rm{D}}t}-{\rm{\nabla }}\cdot {\boldsymbol{T}}\left({\boldsymbol{x}},t\right)={\boldsymbol{f}}\left({\boldsymbol{x}},t\right)\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot v\left({\boldsymbol{x}},t\right)=0,\end{eqnarray} \tag{ 3 }$$

where ρ is the fluid density, ${\boldsymbol{v}}\left({\boldsymbol{x}},t\right)$ is the fluid velocity, ${\boldsymbol{T}}$ is the fluid stress tensor, and ${\boldsymbol{f}}\left({\boldsymbol{x}},t\right)$ represents volumetric forces acting on the fluid (Quarteroni et al 2016). Equation (3) describes the continuity of the blood flow through the vascular network, primarily that the blood is incompressible.

In this work, we simplified the Naiver-Stokes equations by assuming the flow to be steady, laminar, and fully developed. The solution gives the volumetric blood flow rate (Q) in a vessel of radius r and length L as

$$\begin{eqnarray}\begin{array}{rcl}Q & = & \displaystyle \frac{\pi {r}^{4}}{8\eta L}\left({P}_{\mathrm{in}}-{P}_{\mathrm{out}}\right)\\ & = & C\cdot {\rm{\Delta }}P\end{array}\end{eqnarray} \tag{ 4 }$$

where η is the blood viscosity, C is the conductance of the vessel, and P_in and P_out are the blood pressures at the inlet and outlet of the vessel, respectively (Linninger et al 2013). Equation (4) is also known as Poiseuille's equation (Sutera and Skalak 1993). To specify the viscosity, η, of the blood, we used an empirical model for in-vitro blood viscosity adapted from Pries et al (1996) for a constant hematocrit of 0.45, or

$$\begin{eqnarray}&&\eta ={\eta }_{0}\left(220\cdot {e}^{-2.6r}+3.2-2.44\cdot {e}^{-0.06{\left(2r\right)}^{0.645}}\right),\end{eqnarray} \tag{ 5 }$$

where r is the radius of the vessel in millimeters and the value of η₀ was 35 mPa s (Shibeshi and Collins 2005).

Because blood is an incompressible fluid, the net blood flow rate at each junction must be zero, or, mathematically,

$$\begin{eqnarray}&&\displaystyle \sum _{m=0}^{n}{Q}_{m}=0\end{eqnarray} \tag{ 6 }$$

where n is the total number of vessels at the junction and Q_m is the blood flow rate of the mth vessel at that junction (by convention the sign of Q_m is positive if blood is entering the junction and negative if leaving). Using this constraint and equation (4), we constructed a system of linear equations that described the net flow rates at each junction to determine the corresponding fluid pressures. For example, in a network with 6 vessels (figure 2) the corresponding system of equations is

$$\begin{eqnarray}\begin{array}{rcl}{P}_{0} & = & {P}_{\mathrm{in}}\\ {C}_{A}\left({P}_{0}-{P}_{1}\right)-{C}_{B}\left({P}_{1}-{P}_{2}\right)-{C}_{C}\left({P}_{1}-{P}_{3}\right) & = & 0\\ {C}_{B}\left({P}_{1}-{P}_{2}\right)-{C}_{D}\left({P}_{2}-{P}_{4}\right) & = & 0\\ {C}_{C}\left({P}_{1}-{P}_{3}\right)-{C}_{E}\left({P}_{3}-{P}_{4}\right) & = & 0\\ {C}_{D}\left({P}_{2}-{P}_{4}\right)+{C}_{E}\left({P}_{3}-{P}_{4}\right)-{C}_{F}\left({P}_{4}-{P}_{5}\right) & = & 0\\ {C}_{F}\left({P}_{4}-{P}_{5}\right) & = & {Q}_{\mathrm{out}},\end{array}\end{eqnarray} \tag{ 7 }$$

where P_in in the inlet pressure and Q_out is the outlet blood flow rate. We specified the boundary conditions P_in and Q_out to ensure the problem is well-defined. These six equations are used to determine the six unknown pressures at the junctions. This system was cast in matrix form, or

$$\begin{eqnarray}\begin{array}{l}\left[\begin{array}{cccccc}1 & 0 & 0 & 0 & 0 & 0\\ {C}_{A} & -{C}_{A}-{C}_{B}-{C}_{C} & {C}_{B} & {C}_{C} & 0 & 0\\ 0 & {C}_{B} & -{C}_{B}-{C}_{D} & 0 & {C}_{D} & 0\\ 0 & {C}_{C} & 0 & -{C}_{C}-{C}_{E} & {C}_{E} & 0\\ 0 & 0 & {C}_{D} & {C}_{E} & -{C}_{D}-{C}_{E}-{C}_{F} & {C}_{F}\\ 0 & 0 & 0 & 0 & {C}_{F} & -{C}_{F}\end{array}\right]\\ \times \,\left[\begin{array}{c}{P}_{0}\\ {P}_{1}\\ {P}_{2}\\ {P}_{3}\\ {P}_{4}\\ {P}_{5}\end{array}\right]\ =\ {\left[\begin{array}{c}{P}_{{\rm{in}}}\\ 0\\ 0\\ 0\\ 0\\ {Q}_{{\rm{Out}}}\end{array}\right]}_{.}\end{array}\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

In this form, the 2-dimensional matrix mostly describes the conductance values of the vessels that connect junctions while the pressure vector describes the unknown pressures at each junction (e.g., P_i 's) and the vector at the right mostly specifies the net flow rate at each junction. For brevity, the 2-dimensional matrix will be referred to as the conductance matrix for the remainder of the manuscript. We solved this system of equations using iterative methods (van der Vorst 2003) to determine the pressures. With the pressures calculated at all the junctions in the network, we calculated the blood flow rate, Q, through each vessel using equation (4).

We designed a two-step algorithm (figure 3) to calculate blood flow rate in a vascular network of arbitrary size. The first step partitioned the vascular network's geometry and created a computationally-efficient distribution of data using graph partitioning algorithms (see section 2.2.2). The second step calculated blood flow rates through the vessels using matrix inversion techniques (see section 2.2.3).

Zoom In Zoom Out Reset image size

We called the first step Network Preprocessing and the second step Blood Flow Rate Calculations. We developed separate modules for each step to facilitate validation and testing. Additionally, this modular design simplifies future extensions by allowing additional modules to manipulate the vascular data prior to calculating the blood flow.

Limitations on the memory available in a single compute node, as well as matrix solver requirements, motivated the first step of our algorithm. For networks with greater than 67 million vessels, the required amount of memory to describe the vessel network geometry and the conductance matrix exceeded 64 GB available on each of our compute nodes. Our solution was to distribute the rows of the conductance matrix across multiple compute nodes. The matrix solver used in this study stipulated that rows of the conductance matrix be distributed in contiguous groups among the compute nodes (Falgout and Yang 2002), e.g., rows 1 to n on compute node 1, rows n + 1 to 2n on compute node 2, etc. To limit the volume of inter-node communication, the network geometry was partitioned to reduce the number of compute nodes where a given junction appeared.

The algorithm that constructed the network (Donahue and Newhauser 2020) stored the vascular data sequentially in the order of generation number, N_G. If the data was read and distributed among compute between nodes sequentially and contiguously, it resulted in many junctions referenced on multiple compute nodes, reducing computational efficiency (figure 4(a)). We used a graph partitioning algorithm (Karypis 1999) to group and sort the vessels to minimize the number of multi-partition junctions (figure 4(b)). A multi-partition junction is a junction where vessels belonging to multiple partitions meet (figure 4). In practice, each partition is assigned to a different compute node during the Blood Flow Rate Calculations. Thus, to construct a given row in the matrix representing a multi-partition junction, the conductances of the vessels that meet at that junction must be communicated to the compute node that contains that row of the matrix. To facilitate this communication, one of the partitions was selected as the master partition. The master partition was used to decide which compute node was responsible creating the multi-partition row in the conductance matrix, including collecting the conductance values. The lists of immediately connected vessel indexes and the junction indexes were renumbered to correspond with the new partitions . The newly partitioned vessel data was written to the data file for subsequent processing by the Blood Flow Rate Calculation algorithm.

Zoom In Zoom Out Reset image size

Calculating the blood flow rates in the preprocessed network comprised four sub-tasks: importing the vessel data, constructing the conductance matrix, solving for the blood pressure at each junction, and calculating the blood flow rate through each individual vessel. To accomplish the first sub-task, Each compute node used in the simulation was assigned one of the partitions created by the Network Preprocessing application. The compute node imported only the relevant junction indexes, boundary conditions, and vessel geometry from the data file.

We calculated the blood pressures at each junction by constructing the conductance matrix and solving for the unknown pressures. The second sub-task was to construct the conductance matrix (section 2.1.2) in compressed sparse row format. This required the communication of conductance values for vessels that met at a multi-partition junction. The third sub-task used an iterative Krylov method (van der Vorst 2003) to solve the system of linear equations. We applied a boundary condition at the network inlet of P_in = 133.3 Pa, which is the average blood pressure of the cardiac cycle (Aribisala et al 2014), and a boundary condition was assigned to the blood flow rate at the outlet using

$$\begin{eqnarray}&&{Q}_{{\rm{out}}}({N}_{G})={2}^{{N}_{{\rm{G}}}-1}\cdot {v}_{{\rm{c}}}\cdot \pi {r}_{{\rm{c}}}^{2}\end{eqnarray} \tag{ 9 }$$

where ${Q}_{{\rm{out}}}$ is the outlet boundary condition, N_G is the number of generations in the network, v_c is the blood velocity in the smallest capillary (of generation G), and r_c is the corresponding radius of the capillary. In our geometric model, the radius of the smallest capillary (r_c) was 2.5 μm (Duvernoy et al 1981, Lauwers et al 2008, Cassot et al 2010). Ivanov et al (1981) measured the blood flow velocities in capillaries at approximately 1 mm s⁻¹, which we assigned to v_c.

To accomplish the final sub-task, the blood flow rate through each vessel was calculated according to the following procedure. First, the pressure for each of the multi-partition junctions was communicated to all the compute nodes sharing a given junction. Then, the blood flow rate of each vessel was calculated using equation (4). To optimize this procedure for computational speed and efficiency, we converted the problem into a matrix-vector multiplication operation.

In order to verify the calculations were free of blunders, each compute node calculated the blood flow rate using the following secondary method. For the specific-case of symmetric vascular networks used in this study on the assumption of steady-state flow, we calculated the blood flow rate (Q) of each vessel of network according to

$$\begin{eqnarray}&&Q=\displaystyle \frac{{Q}_{{\rm{out}}}({N}_{G})}{{2}^{l-1}}\end{eqnarray} \tag{ 10 }$$

where l is the generation number of the vessel. This simple procedure allowed us to verify the calculations for networks with up to 17 billion vessels. Following this verification, the blood flow rates were written to persistent storage media.

To assess the accuracy of the calculated blood flow rates (equation (4)), we compared the blood flow rates calculated by the steady-state model to those predicted by a commercial computational-fluid-dynamics (CFD) package (Autodesk CFD 2018, San Rafael, CA, USA). The commercial package used a finite-element model to calculate the blood flow using the Naiver-Stokes equations. This package was previously validated (Alexa et al 2014, Autodesk 2015).

We used two different types of vessel network configurations to validate the algorithm. The first types was a healthy vascular network, where all vessels were open (figure 1). For the healthy networks, we calculated the flow rates through networks with between 6 and 126 vessels (2 ≤ N_G ≤ 6). The second type was a damaged network, where a single vessel was completely closed, i.e., the vessel radius was set to zero. For the damaged networks, we used a 126 vessel network (N_G = 6) calculated the flow when three different vessels were closed (figure 5). The networks sizes were limited by the computational expense of CFD calculations.

Zoom In Zoom Out Reset image size

To perform the calculations in the CFD package, it was necessary to convert the geometry created by the fractal method (Donahue and Newhauser 2020) into a fully specified 3-D model. This was done using a commercial computer-aided design (CAD) package (Inventor 2018, Autodesk, San Rafael, CA, USA). First, we imported the bifurcation locations and vessel radii from our model geometry. Next, linearly-tapered vessels were constructed between junctions. Tapering vessels avoided nonphysical sharp edges at the bifurcations in the CFD model. The initial radius of each vessel matched the radius of the corresponding generation in our model while the final radius was the radius of the next generation. We chose this approach to simplify the implementation of the tapers. The vessels in the final generation had a constant radius of 2.5 μm. We converted the vessel network geometries manually. Closed vessels were completely removed from the geometric representation in the CAD model. Using the CFD package, we specified the material properties of the blood and vessel walls using the built-in materials library. Automated tools, provided by the package, were used to compute the simulation mesh for each geometry. We used the default solver of the CFD package to perform the blood-flow calculations. The solver was set up to perform laminar calculations of the blood flow, in accordance with the assumptions of the steady-state model. Additionally, it was found that the Reynolds number was less than 100 for all networks simulated, implying that turbulent flow was not present. Table 1 lists the mesh generation and solver settings for the commercial package. We used the boundary conditions specified in the steady-state model (section 2.2.3) for the CFD package calculations.

We recorded the calculated flow rates from the steady-state model (equation (4)) and the CFD package. Our model calculated a single value of the blood flow rate through each vessel. For the CFD package, which calculated the distribution of blood flow rates throughout each vessel, we recorded the blood flow rate at the middle of each vessel. Theoretically, blood vessels directly downstream of a closed vessel have a zero flow rate. The iterative solvers used in both our algorithm and the CFD package reported negligible but non-zero values in these vessels that were caused by limitations in machine precision (e.g., rounding and truncation errors). We removed these artifacts by setting blood flow rates less than 1 μm³/s to zero. This value, which is 20,000 times smaller than the theoretical minimum flow rate in the network (1.9 × 10⁴ μm³ s⁻¹), is negligible. We computed the unsigned absolute ($\left|{{\rm{\Delta }}}_{\mathrm{abs}}\right|$) and relative ($\left|{{\rm{\Delta }}}_{\mathrm{rel}}\right|$) differences in blood flow reported by the two calculation methods.

Validation measurements were performed using a personal computer. The system had a 4-core 4.1 GHz CPU (Core i5-3570K, Intel Corporation, Intel Corporation, Santa Clara, CA, United States) and 32 GB of memory. The operating system was Windows 10 (Microsoft, Redmond, Washington, United States).

We measured the performance of the algorithm on the QB2 cluster hosted by the Louisiana Optical Network Initiative. This shared resource has 504 compute nodes. Each compute node contains 64 GB of memory and two 10-core processors (2.8 GHz Ivy Bridge-EP E5-2680 Intel Xeon 64-bit). An Infiniband Interconnect (Pfister 2001) provides internode communication and access to persistent storage, comprising a 2.8-PB parallel file system. The file system uses a Lustre architecture (Schwan 2003), enabling up to 16 compute nodes to operate on a single file simultaneously.

We compiled all software developed in this study with a commercial compiler (Intel C++ Compiler Cluster Edition v18.0, Intel Corporation, Santa Clara, CA, United States). The applications utilized shared-memory parallelism, distributed-memory parallelism, and linear algebra routines. We used OpenMP (Dagum and Menon 1998) for shared-memory parallelism, an implementation of the message passing interface (MPI) for distributed memory communication (Intel MPI v18.0), and the Math Kernel Library (MKL v18.0) of optimized linear algebra routines. These were included in the compiler.

All input and output (I/O) operations to persistent storage for our applications utilized the Hierarchical Data Format version 5 (HDF5) library (Folk et al 1999). Our applications used the parallel I/O capabilities of the HDF5 library and the Lustre file system to increase speed. We tuned the I/O operations following the suggestions of Howison et al (2010).

Graph partitioning in the Network Preprocessing application utilized a k-way partitioning algorithm implemented in the parallel METIS library (ParMETIS version 4.0.3) (Karypis and Kumar 1998, Karypis 1999). We used the k-way partitioning routine on an edge graph of our vascular network. The ParMETIS library was compiled with support for large indexes (64-bit) and double-precision real (64-bit) values to accommodate networks with up to 2.3 × 10¹⁸ vessels.

Matrix solving was completed with the Loose Generalized Minimum Residual algorithm (LGMRES) (Baker et al 2005), a Krylov-based iterative solver. We used the BoomerAMG solver (Henson and Yang 2002) as a matrix preconditioner to accelerate the convergence of the iterative solution. Implementations of these algorithms were provided by the Hypre Solver Library (Falgout and Yang 2002, Stüben 2001). We compiled the library with support for OpenMP and large (64-bit) indexes. Table 2 lists the user-specified parameters and options for the matrix solver and preconditioner. These parameters were determined through an iterative tuning process.

One of the primary metrics for computational feasibility and utility is execution time. We measured the execution times to calculate the blood flow rate in networks containing between 14 vessels and 17 billion vessels (3 and 33 generations, respectively). Each calculation was repeated 5 times to quantify the variation in execution time caused by resource contention in the computational cluster. For large vessel networks (i.e. N_G > 25 and N_V > 67,000,000), the volume of data required to describe the network and perform the calculations exceeded the memory of a single compute node (64 GB). The number of compute nodes, N_c, used to perform the calculations was

$$\begin{eqnarray}{N}_{{\rm{c}}}=\left\{\begin{array}{ll}1, & \mathrm{for}\ {N}_{{\rm{G}}}\leqslant 25\\ {2}^{({N}_{{\rm{G}}}-25)}, & \mathrm{for}\ {N}_{{\rm{G}}}\gt 25.\end{array}\right.\end{eqnarray} \tag{ 11 }$$

The Network Preprocessing application was only required for large networks (N_G > 25). Computational execution times were recorded for the Blood Flow Rate Calculations and the Network Preprocessing independently. The total time, vessel-network-partitioning time, vessel-index-renumbering time, junction-renumbering time, and I/O time were recorded for the Network Preprocessing application. For the Blood Flow Rate Calculations, the time to complete the I/O operations, conductance matrix construction, conductance matrix solution, blood flow rate calculation, and total times were recorded.

Weak and strong scaling are two metrics used to evaluate the scalability of an algorithm in high-performance computing (Gustafson 1988). Weak scaling characterizes the execution time as the number of processors increases when the workload per processor is fixed, e.g., as the number of processors doubles so does the workload. Ideally, as the number of processors increases the wall-clock execution time remains constant. We measured the weak scaling for both Network Preprocessing and Blood Flow Rate Calculations. For the Blood Flow Rate Calculations, the smallest vessel network used was N_G = 25 on a single compute node and doubled the number of vessels and compute nodes up to N_G = 33 on 256 compute nodes. Because the Network Preprocessing was only used for more than one compute node, we tested the scalability in the range from N_G = 26 on 2 compute nodes and to N_G = 33 on 256 compute nodes.

Strong scaling characterizes how adding more processors reduces the computation time of a fixed size problem. This is typically measured using the speedup factor. Ideally, the speedup factor is equal to the ratio of the new number of processors to the original number used for the problem (Gustafson 1988). The strong scaling was determined for the Blood Flow Rate Calculations by constructing a fixed-size vessel network (N_G = 25) and varying the numbers of compute nodes (1 ≤ N_c ≤ 256). Strong scaling of the Network Processing considered a 26-generation network, and 1 ≤ N_c ≤ 256. Strong scaling was quantified by calculating the speedup factor (S) according to

$$\begin{eqnarray}&&S=\displaystyle \frac{{T}_{q}}{q\cdot {T}_{n}}\end{eqnarray} \tag{ 12 }$$

where T_q is the time it takes a reference number of processors, q, to complete the calculation and T_n is the time it takes n processors to perform the same calculation. For the Blood Flow Rate Calculations, q = 1. The Network Preprocessing had a reference number of processors of q = 2.

Memory requirements for Network Preprocessing and Blood Flow Rate Calculations were collected separately using a heap memory profiling tool (Valgrind Massif version 3.13 (Nethercote and Seward 2007)). We measured the peak and average memory usage for vessel networks with 3, 13, 19, 25, 26, 28, 30, and 33 generations, which spanned the operational range of our two-step algorithm and compute cluster (1 ≤ N_c ≤ 256). For network sizes requiring two or more compute nodes, the memory used per node was recorded.