Parallel generation of extensive vascular networks with application to an archetypal human kidney model

Given the relevance of the inextricable coupling between microcirculation and physiology, and the relation to organ function and disease progression, the construction of synthetic vascular networks for mathematical modelling and computer simulation is becoming an increasingly broad field of research. Building vascular networks that mimic in vivo morphometry is feasible through algorithms such as constrained constructive optimization (CCO) and variations. Nevertheless, these methods are limited by the maximum number of vessels to be generated due to the whole network update required at each vessel addition. In this work, we propose a CCO-based approach endowed with a domain decomposition strategy to concurrently create vascular networks. The performance of this approach is evaluated by analysing the agreement with the sequentially generated networks and studying the scalability when building vascular networks up to 200 000 vascular segments. Finally, we apply our method to vascularize a highly complex geometry corresponding to the cortex of a prototypical human kidney. The technique presented in this work enables the automatic generation of extensive vascular networks, removing the limitation from previous works. Thus, we can extend vascular networks (e.g. obtained from medical images) to pre-arteriolar level, yielding patient-specific whole-organ vascular models with an unprecedented level of detail.


Introduction
Nowadays, current multiscale modelling approaches and computational resources have stretched the boundaries of cardiovascular models and hemodynamic simulations. These in silico tools are starting to become part of medical research and play a major role in understanding the intricate physiological phenomena, developing novel therapeutic strategies, and improving diagnostic procedures [1][2][3][4].
Hemodynamics in the cardiovascular system (CVS) is a challenging problem involving blood flow dynamics and physiological phenomena occurring at different scales. For instance, the microcirculation determines the peripheral hemodynamic environment, which affects the macro-circulation through the local and global regulation of peripheral resistance and compliance. In turn, the macrocirculation dictates the hemodynamic forces that promote the physiological interactions that unfold at the micro-circulation scale.
On the one side, complex three-dimensional models have gained momentum for the study of refined blood flow dynamics in localized regions of the CVS [5][6][7]. On the other side, one-dimensional models have been preferred to assess the hemodynamics at the scale of the entire CVS, ranging from large vessels to small arteries [8][9][10]. Most modelling approaches available in the literature employ simplified mathematical representations for the peripheral vascular beds [11][12][13]. In this sense, the oversimplification of the micro-circulation does not allow the study of sophisticated physiology at the level of smaller vascular structures, such as the size-dependent arterioles vasodilatory/remodelling responses. In turn, several works have studied the circulation at the scale of small arterioles and capillaries [14][15][16][17] employing lumped models (0D models: ordinary differential and algebraic equations), focusing on small pieces of tissue to analyse large vascular networks and investigating the relationship between blood flow, pressure, shear forces and mass exchange.
Medical images serve as input data to generate vascular domains. According to the modality, imaging techniques enable the visualization of a wide range of vascular vessels. Technological limitations do not allow us to analyse blood flow in the micro-circulation in vivo. Hemodynamics at the scale of arterioles and capillaries was investigated through mathematical models operating on top of vascular models constructed through high-resolution micro-tomographic images obtained from animal models [18,19]. In this context, automatic vascularization algorithms emerged more than two decades ago as a systematic approach to generate networks of interconnected vessels in regions of interest. These algorithms follow the hypothesis that network topology and geometry achieve an energetically efficient (or semi-efficient) perfusion of the vascularized tissues.
Several strategies have been devised to achieve that goal, from fractal approaches [20,21] to spacefilling methods based on Voronoid tesselations [22,23] and cost function optimization methods [24][25][26][27]. The latter methods are termed the constrained constructive optimization (CCO) approach. Based on these CCO methods, several ramifications have been proposed to evolve in terms of model flexibility [28][29][30], enabling morphometric studies and blood flow simulations at the smaller scales in the CVS. Moreover, building in vitro vascular networks from scratch is being pursued by current technological advances [31][32][33]. Prior knowledge about network topologies suited to improve mass exchange and promote growth signalling is desirable for the optimized construction of these networks. In this experimental context, CCO-based methods can also provide insight into the advantage of certain topological configurations over others.
However, CCO-based approaches rely on a sequential optimization procedure in which points are sorted out, and a connectivity decision is made regarding the minimum of a cost function. These methods have been successfully employed to generate networks of up to 20 000 vascular segments. The main drawback in reaching larger numbers lies in the sequential approach and the ever-growing computational burden related to memory copying, neighbour searching, and network scaling when adding a new vessel. This has confined the possibilities of studying networks with a large number of vessels. Hence, modelling organ function accounting for the specific cellular physiology and its coupling with the systemic hemodynamic environment has remained an open challenge.
In light of the previous context, the main novelty of the present work is to develop a parallel strategy to efficiently generate extensive vascular networks using a CCO-based approach. Given a spatial domain to be vascularized, the first stage consists of creating a baseline network using a recently proposed aDaptive Constrained Constructive Optimization method (DCCO) [30]. At a second stage, we split the domain into non-overlapping subdomains, each of which is independently vascularized using the baseline network as the underlying vascular substrate. At a third stage, we consistently merge the networks. For further reference, we name the proposed parallel approach as PDCCO. The proposed PDCCO strategy exploits the multiscale nature of the CVS, for which the macro-scale circulation is royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 approached using a conventional CCO-based strategy, while the refined vascular networks at the microscale circulation are generated in a decoupled fashion, understanding the vascular generation problem as a local phenomenon.
Several numerical experiments are reported to investigate the capabilities and limitations of the proposed method. By comparing with sequentially generated networks (DCCO-networks), we study: (i) the impact of subdomain shapes in geometric and functional statistical indexes; (ii) the impact of the size of the baseline network in the resulting vascularization and (iii) the impact of the cost functional in the vascularization of a domain featuring multiple inlets. To illustrate the potential of the proposed method, we generate the arteriolar network in a piece of tissue resembling a portion of cortical renal tissue featuring 100 000 vascular segments. Finally, as an application towards the simulation of whole-organ function, we propose to create the arterial network in an archetypal model of the human kidney. Specifically, we construct the renal vasculature from scratch, starting at the renal artery until reaching the interlobular vessels, where afferent vessels are connected to glomeruli where filtration occurs. Specifically, we report the generation of an arterial tree in a highly non-convex geometry representing the renal cortex, featuring 100 000 terminal segments.

Material and methods
In this section, we briefly revisit the DCCO algorithm, and describe the proposed partitioned DCCO, termed PDCCO. Also, we describe the statistical analysis employed to compare vascular networks produced by PDCCO and DCCO.

Adaptive constrained constructive optimization
In turn, V and T are defined as the sets of admissible vessels and trees, respectively. Given a vessel cost functional F : V ! R, we define the total tree cost F : T ! R as Most commonly F measures total intravascular volume, i.e. F(v) = π|x d − x p |r 2 . Nonetheless, we can change the optimality criterion by modifying F. For instance, to account for the energetic cost associated with angiogenesis [30], F is defined as follows: where V ref and l ref are characteristic volume and length of the perfusion domain, r ref is a characteristic radius for the arterial tree, c v , c p and c d are, respectively, the volumetric, proteolytic and diffusion nonnegative coefficients such that c v + c p + c d = 1. Parameter γ characterizes the branching power law which governs the relation between parent and daughter vessels. Note that for γ = 3, we recover the classical volumetric cost functional, while for other values of gamma, the functional is consistent with the power law employed in the vascularization process.
In the CCO method, we search for T such that it minimizes a cost functional related to the genesis and maintenance of an arterial tree, namely As initial condition, the method requires a initial tree T base or the position of an inflow point, radius and supplied flow rate, x p 0 , r 0 , Q 0 .
If we cannot find any valid connections, we choose another random point and proceed again. Otherwise, we add to the enduring tree the vessels that resulted in the smallest value of dF .
Along the optimization process that underlies the vascularization, it is necessary to solve compartmental fluid mechanics equations in the network. In doing this, conventional hypotheses are steady-state regime and Poiseuille blood flow. Since we are targeting the generation of massive networks, we consider the model of blood rheology that accounts for the F _ ahraeus-Lindqvist effect [34]. This implies viscosity depends on the vessel diameter, which introduces a nonlinearity in the vascularization process. We use fixed-point iterations as linearization procedure. This approach proved to be convergent and stable for all the experiments reported in this study. However, it is important to highlight that lack of convergence can be obtained for unrealistically small vessels (smaller than capillaries), for which the viscosity model is also questionable. For a detailed explanation of the DCCO algorithm, the reader is referred to [30].
The sequential optimization required by CCO and DCCO algorithms becomes increasingly computationally intensive as more vessels are added to the vascular tree. This aspect poses a limit to the utilization of these algorithms for the generation of extensive vascular networks, where a few recent works have studied how to alleviate such restrictions [18,35]. In the next section, we propose a divide and conquer approach for the optimization procedure to widen the applicability of this type of vascular tree generation algorithm.

Domain decomposition strategy
Let V be the perfusion domain with a region V 0 # V concurrently vascularized and D ¼ fV 1 , . . . , V Npart g a partition of V 0 . The domain decomposition strategy is fully described in algorithms 1 and 2, while figure 1 conceptually illustrates the major steps involved in the process. The proposed method follows the steps described below: (i) Create a baseline tree, T base , with N base terminal segments in V. This baseline tree accounts for the global macro-scale features of the circulation. (ii) For each V i [ D, using T base as a starting point, add N i terminal segments in V i to determine T i through the DCCO algorithm presented in §2.1. It is important that we consider only vessels v ∈ T base whose midpoint resides in V i . We call the set of such vessels V part , as potential parents for the new segments in T i . Thus, v cannot bifurcate in more than one subdomain, guaranteeing the consistency of the merged tree.   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 blood flow in such regions is homogeneously redistributed, thus, terminal blood flow distribution outside V i remains unchanged. This strategy yields a more consistent common tree across the parallel generated trees T i . That is, we maintain q out fixed for terminal vessels not in V part , and distribute the remaining flow Q part homogeneously for terminal vessels inside the domain. Such boundary conditions come into force on lines 22 and 30 of algorithm 1. This parallel stage accounts for the refined generation of micro-scale circulation networks, understood as a local phenomenon that does not significantly affect the macro-scale features of the circulation (see §3). (iii) We initialize T merged as T base and, for each L i , we sequentially add the vessels (x new , x bif ) whose parent is ðx p p , x d p Þ in the listed order. The radius of the added vessels are only scaled after the insertion of all tuples in the corresponding lists L i . So as not to search the whole tree for each insertion, we use a hash table whose keys are x p and x d to store the values x new and x bif . (iv) Once the tree T merged has vessel radii are scaled according to the hemodynamic and geometric constrains as detailed in [30]. This is described in line 13 of algorithm 2. The merging process consists in coupling the macro-and micro-scales by providing consistency to the concurrent networks as a whole unit.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 The simultaneous determination of the optimal points for each subdomain V i , as explained in Step 2 (see also figure 1), significantly reduces the computational burden, which amounts to the most computationally intensive stage in the sequential algorithm.

Data analysis
Vascular networks constructed using the PDCCO approach are compared to a vascular network built using the sequential DCCO approach. Given the random strategy to sample potential points, several instances are created for each case to properly characterize the stochasticity of the process and perform a fairer comparison between PDCCO and DCCO. From the set of instances, either in the sequentially constructed network (DCCO network) or in the parallel approach (PDCCO network), we compare morphometric and functional measures, specifically vessel radii, length, aspect ratio and pressure, as a function of the vessel generation. Pressure is computed through Poiseuille Law, consistently with the underlying physics in the model. We place the reference (null) pressure at the inlet of the network. In doing this, we gather, across the different instances all vessels of a certain generation and build the corresponding box-plot. We report profiles of these quantities as a function of vessel generation. Also, we compare the distribution of vessels in the length-radii space.
As further verification of the vascular tree optimality resulting from the PDCCO algorithm, we report the discrepancy between the intravascular volume in PDCCO-networks and the mean volume of all DCCO networks. This discrepancy is normalized by the mean DCCO network intravascular volume and is termed relative volume error.
For the prototypical model of the kidney vascular network, we compute the Strahler order of each vascular segment to provide an alternative characterization of the radius, pressure and flow rate in the network. From vessel radii, we compute the total cross-sectional lumen area of all vessels combined, within each Strahler order bin, and the intravascular volume of the subtended tree associated with each vessel. Lastly, we construct a connectivity matrix, as a function of the Strahler order [36].

Results
In this section, we first evaluate the performance of the parallel approach (PDCCO) compared with the sequential approach (DCCO) in three idealized experiments, and later we analyse the performance of PDCCO to generate complex trees with a large number of vessels. Unless stated otherwise, we make use of the volumetric cost functional to generate the PDCCO and DCCO vascular networks.
The first experiment investigates the sensitivity in the resulting tree to the number of segments N base present in the baseline tree, from which the parallel vascularization in PDCCO is initiated. The second experiment refers to the effect of the domain shape, particularly the fact that partitions can become royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 geometrically degenerated, in the resulting partitioned ( parallel) vascularization. For these two experiments, the baseline network is shared by both approaches (DCCO and PDCCO). Therefore, the sequential approach continues the construction of the baseline network until reaching the total number of terminals, while the parallel approach uses the baseline network to initialize the concurrent vascularization of the different subdomains. In the third experiment, we study the impact of the cost functional definition in the resulting vascular network. Specifically, we propose to vascularize a domain with multiple inlets, and investigate the lack of sub-tree balance, in terms of flow rate carried by each inlet, as a function of the parameters that define the cost functional.
Regarding the PDCCO performance generating complex trees, we choose to study the scalability of the proposed approach in a prototypical portion of tissue with a parallelepipedal shape. The goal is to generate a vascular network in which over 100 K vascular segments are placed. Finally, an example of application in which we vascularize a renal cortex structure demonstrates the potential utilization of the proposed approach in a challenging scenario.

Sensitivity to the baseline tree
To analyse how the number of terminal segments in the baseline tree (N base ) affects the final network, we executed nine sets of simulation scenarios, each set with a different N base value. All executions in a set use the same parameters except for the random generation seed used for the point generation (see Algorithm 1, lines 8 and 10). Each set has 10 instances of vascular networks with 5000 terminal vessels generated by the PDCCO approach, aiming to characterize the stochasticity in the generation process due to the point generation. The domain geometry is defined as the parallelepipedal region from figure 1 divided into N part = 4 subdomains. The number of terminal segments in the baseline tree for each set is given by N base ∈ {500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500}. Also, we built ten instances of trees using the sequential DCCO approach to use as baseline for comparison.
In figure 2, the relative volume error is reported (a) as a function of the number of segments in the baseline network (N base ). Recall that this error measures the normalized discrepancy between the intravascular volume in PDCCO-generated networks and the mean of DCCO-generated networks. Observe that the median error approximates a plateau for N base ∈ {500, 1000, 1500, 2000, 2500}, and then consistently declines as the number of elements in the baseline network is incremented. As expected, PDCCO renders a suboptimal solution for the optimization problem due to the decoupling royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 of the problem. Nonetheless, the median of the disagreement with respect to DCCO is below 8% in the luminal volume of the generated tree. A complemental comparison between PDCCO and DCCO generated networks is given in the right panel of the same figure, where the distribution in the feature space defined by the vessel length and radius is shown. Observe that both PDCCO networks and DCCO networks share the same structural features in terms of vessel radius and length. Moreover, this common behaviour is independent from N base . This piece of evidence suggests that the morphometric features of the vascular networks are invariant. Figure 3 displays the statistical behaviour of the vessel radius and blood pressure depending upon the vessel generation (i.e. the number of bifurcations from the root vessel, v 0 ) in the vascular network. For each vessel generation, the box-plots given by the PDCCO-networks and by the DCCO-networks are compared head-to-head for three selected numbers of segments in the baseline tree N base ∈ {500, 2500, 4500}. It can be observed that the radius interquartile ranges are in agreement between both PDCCO and DCCO algorithms along with the network depth. Similarly, whiskers are also in agreement, with subtle shifts in the middle part for N base = 2500. In turn, the pressure has a very compact distribution in the first generations, and then the distribution widens. This feature is shared by both DCCO and PDCCO algorithms regardless of N base . In more detail, we can observe that the pressure drifts upwards after the 6th generation in the PDCCO networks, compared to the DCCO networks. In other words, the pressure drop as predicted by PDCCO networks is slightly smaller than the pressure predicted by DCCO networks.  6 8 10 12 14 16 18 20 22 24 26   0 2 4 6 8 10 12 14 16 18 20 22 24 26  0 2 4 6 8 10 12 14 16 18 20 22 24 26   28  0 2 4 6 8 10 12 14 16 18 20 22 24 26 28   0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32  0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 10 -1 royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973

Sensitivity to the domain shape
We have designed an experiment to check whether the subdomain's shape functionally affects the generated vascular network. As with the first test, we ran 10 instances for each setup, all with the same parameters except for the random number generator seed. In total, we have six different setups where PDCCO and equivalent DCCO vascular networks are generated. For the PDCCO trees, we had N part = 4 and N base = 2500 with the final tree having 5000 terminal segments. To control the effects of changing the domain volume, we scaled the influx and root radius according to the allometric laws Q ∝ V and r 0 ∝ V 3/8 [37]. The domain is defined by a cuboid of dimensions 4L × 1 × 1 with four cuboid subdomains of dimensions L × 1 × 1 adjacently disposed on the vertical axis, as shown in figure 4. The influx position is taken as (L/10). We choose L ∈ {0.25, 0.5, 0.75, 1.0, 1.5, 2.0} so as to have the subdomain aspect ratio A r ∈ {0.25, 0.50, 0.75, 1.00, 1.50, 2.00} for our six different setups. Figure 5a displays the relative volume error (discrepancy between the intravascular volume in PDCCOgenerated networks and the mean of DCCO-generated networks). Observe that the median error for the different values of A r (aspect ratio) remains bounded, below 8%. It is worthwhile to note that the error is smaller for more degenerated subdomains (A r ∈ {0.25, 0.50}), grows and then stagnates for A r ∈ {1.00, 1.50, 2.00}. In figure 5b, we illustrate an example of the baseline network, the DCCO network and the PDCCO network for different values of A r . In figure 5c, the radius-length distribution is shown for the DCCO and PDCCO networks. Again, in this case, the PDCCO algorithm manages to deliver networks structurally similar to the DCCO algorithm. Also, note that by modifying A r , the change in the distribution shape is minimal. That is, the PDCCO algorithm emulates the overall morphometric behaviour of the DCCO network.
In figure 6, vessel radius and blood pressure distributions are reported as a function of the vessel generation. The box-plots for each generation from the generated PDCCO networks and DCCO networks are compared for three selected aspect ratios A r ∈ {0.25, 1.00, 2.00}. Radius interquartile range in PDCCO-generated networks closely follows the behaviour of the networks generated by the DCCO algorithm. Only subtle discrepancies are noticed at most distal locations, where the number of vessels per box drops significantly. Whiskers, in turn, are similar in the first generations and start to deviate as we go over the 10th generation. That behaviour is independent of A r . Due to the highly nonlinear relationship between pressure and vessel radius, the subtle differences observed in the left column plots become more evident in the right column plots. Remarkably, the   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 overall distribution shape of the pressure profile is in agreement in PDCCO networks and in DCCO networks for the different values of A r . A close-up analysis shows us that for A r = 0.25, the pressure range is greater in PDCCO networks for a given vessel generation, and the median remains slightly above that in DCCO networks. For larger values of A r , the pressure range drifts upwards, including the whole interquartile range, in PDCCO networks. Hence, the pressure drop is slightly smaller in PDCCO networks than in DCCO networks.

Domain with multiple inlets
In this section, we build a vascular network in a domain supplied by multiple inlets. We consider a cubic domain centred at the origin (side length of 2) and an external network supplying the domain through four opposed planes ( figure 7). The tuple of parameters (c v , c p , c d ) that characterize the cost functional (see expression (2.2)) is defined by taking steps of 0.1 for all the parameters, from 0 to 1 for each parameter, with the constraint that c v + c p + c d = 1. Hence, we build a total of 55 baseline tree vascular networks, and investigate the flow rate distribution among the four major supply branches. Out of these 55 baseline trees, we analyse two exemplar cases: (i) a close to an 'equal' flow distribution; and (ii) an 'unequal'    The scenario termed 'equal' is the one for which the flow rate standard deviation in the four branches achieves the smallest value. This leads to inlet branches with nearly equal size. The scenario called 'unequal' is the one for which the flow rate standard deviation in the four branches reaches the highest value, resulting in inlet branches with unequal size, and networks whose extension varies substantially. Next to each inlet branch, we inform the corresponding flow rate fraction. We also report in the same figure the radius and pressure distributions for these two extreme cases.
Macro-scale features of the vascular network are determined in the generation of the baseline network, which is accomplished using the sequential DCCO algorithm. This structure depends on the cost functional parameters as expected. The PDCCO approach initiates the process with an already balanced/unbalanced tree, and it continues to grow within the domain in a parallel fashion, regardless of the extension of each vascular subnetwork and flow rate balance. This is confirmed by the subtle changes in the flow distribution when moving from the baseline tree T base to the merged one T merged . Hence, the partitioned stage is responsible for the definition of refined vascular features encountered at smaller scales. Indeed, as also shown in figure 7, the distribution of vessel radius and pressure as a function of the vessel level changes substantially in the first levels, although higher levels present similar values.

Scalability
In this section, we design an experiment to test the scalability of the parallel approach to generate a vast network. This allows us to study the computational gain obtained with the proposed strategy. Let us consider the parallelepipedal region from figure 8. The domain is partitioned into 16, 25 and 36 equally sized subdomains, and the number of terminal vessels in the baseline network is N base ≈ 5000-this is to keep the remaining 100 000 − N base vessels as a number divisible by N part . First, we royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 generate the baseline tree in the entire domain, and then the parallel generation is executed, after which the N part networks are merged to deliver the final network. The flow rate at the inlet of the network is set to be Q in = 2 mm 3 s −1 and the inlet radius is r 1 = 0.075 cm. Figure 9 displays the resulting vascular networks for the different partitions employed. Over the top of the figure, by thresholding the network with the vessel radius, it is possible to see that the dominant vascular structures remain unchanged when altering the partitioning and the parallel generation of the vascular tree. This outcome is fundamental to confirm that the parallel version of the DCCO algorithm does not introduce artefacts in the resulting network. The equivalence between the sequential and the parallel algorithms can also be verified in the radius and pressure profiles reported at the bottom of the figure. The box-plot for the three different partitions are compared alongside as a function of vessel generation. Both geometrically and functionally, the three networks behave equivalently. Over the left panel in figure 9, the intravascular volumes of the different networks are compared, where no  For this experiment, we measure the cost in terms of time spent to address the different phases in the algorithm. Table 1 reports the time spent (in hours) to build the baseline network (N base = 5000), as well as the time spent in the construction of the vascular networks for each subdomain (mean and s.d.). Note that, when removing the shared baseline network, the PDCCO algorithm scales linearly concerning N part as the subdomain vascularization problems are completely decoupled.
To compare the estimated cost in building such a network using the sequential DCCO algorithm, starting from the baseline network, we executed the sequential algorithm for different numbers of terminal segments, specifically {1, 2,3,4,5,6,7,8,9,10,15,20,25,30, 35} · 10 3 terminals. Regression analysis yielded a cost of the order of N 3 seq , wirh N seq being the number of terminal segments added to the baseline network sequentially. Then, we estimated the cost of sequentially adding N seq = 95 000 terminals to the baseline network. This estimated cost is also reported in table 1.

Renal vascular network
In this last experiment, the goal is to illustrate the application of the PDCCO method in a complex scenario: the vascularization of the renal cortex structure. To this end, we have created a prototypical geometry of a left human kidney from reported anatomical data [38]. For a body surface area of 1.65 m 2 [8,39] and from the data reported in [38], we determined total, cortical and medullary volumes for our model. Then, we scaled the geometric model until achieving appropriate dimensions. Morphometric data for the kidney model is detailed in table 2. From the point of view of the PDCCO algorithm, this is the input data that characterizes the domains to be vascularized.
At the beginning of the renal vasculature, the renal artery bifurcates into many segmental vessels, that continue to branch into interlobar arteries that extend through the renal columns surrounding the renal pyramids, to finally give rise to the arcuate arteries. Arcuate arteries branch into interlobular vessels. The renal cortex is one of the most densely vascularized regions in the human body. Within 72 cm 3 , up to 1  million glomeruli can be found. Glomeruli are connected to tubules, forming the nephrons, which are the basic functional units responsible for blood filtration. Each glomerulus is a capillary bed contained in Bowman's capsule, where the blood crossing the glomerular membrane filtration turns into the socalled filtrate. An afferent arteriole supplies each glomerulus (diameter in the order of 25 μm), which in turn branches off from the aforementioned interlobular arteries (diameter in the order of 50− 100 μm). Such a compact vasculature is guested in a spatial domain containing many obstacles. The set formed by renal pelvis, calyxes and pyramids poses intrinsic spatial restrictions from early embryonic stages that constrain the development of the vascular network. The initial vascular model, the main components of renal anatomy, and the initial vascular network, from which the proposed PDCCO algorithm is applied, are displayed in figure 10. In this experiment, the goal was to create a vascular network of the entire kidney. The kidney vascularization was achieved through a multi-stage strategy as described in [30], and detailed in table 3.
All the main ingredients involved in the process are schematically presented in figure 10. First, an initial renal vascular network containing the three major arterial vessels, inferior, anterior and medial segmental arteries (see the three-vessel network in figure 10), was handcrafted from the renal artery [8,39]. Then, we triggered the first vascularization stage from these three major segmental vessels (network T 0 ) using a sprouting cost functional with angle constraint and spanning the medulla domain, called V core (see purple region in figure 10). The second and third stages are responsible for placing 5000 terminal segments in the V cortex region (see region which surrounds V core in figure 10). These stages, and the ones that follow, optimize a volumetric cost functional to drive the vascularization. The difference between these two stages lies in the size of the neighbourhood (determined by f n ), where potential connections between the new terminal point and each vessel are sought. The outcome of these two stages is the vascular network T 3 , which is used as the baseline network to launch the parallel partitioning proposed in this work. The kidney domain is divided into eight parallelepipedal subdomains placed at different transverse plane locations (see vascularization of V i , i = 1, …, 8 in figure 11). Four stages with progressive reduction in the neighbouring factor f n were used to vascularize the V cortex domain, until reaching a vascular network featuring 100 000 terminal vessels.
Concerning the setting of the PDCCO algorithm, the baseline network (called T 3 in table 3) contains N base = 5000 terminal vessels. As said above, eight partitions were created by slicing the kidney with royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 Table 3. PDCCO parameter setup for kidney vascularization. Coefficients and reference values for (1.0, 0.9, 1.0, 7) (1.0, 0.9, 0.5, 7) (1.0, 0.9, 0.25, 7) (1.0, 0.9, 0.125, 7)  Figure 11 illustrates the result of the PDCCO algorithm for the vascularization of a human kidney. On the top of the figure, the baseline network containing 5000 vessels is shown. This network was generated sequentially. After this first stage, the vascularization is parallelized into eight subdomains. The individual vascular networks are also shown next to the baseline network. Finally, the merging procedure is illustrated on the right, where the final network is assembled, reaching 100 000 terminal vessels. In the middle of the same figure, the radius and pressure profiles are displayed as a function of vessel generation. The pressure drop in the network is around 5 mmHg. On the right, the radius, pressure, and flow rate are displayed across the entire tree. Finally, at the bottom of the figure, we report a sequence of panels that allows us to visualize the resulting network structure by thresholding the merged tree with the vessel radius. Note the hierarchical arborization surrounding the pelvis and pyramids to finally reach the cortex, where the arterial vessels are closely packed.
In figure 12, we show several quantities classified in terms of the Strahler order. Specifically, that figure displays the total cross-sectional lumen area associated with all the vessels with a given Strahler order, the degree of connectivity of vessels, the pressure and flow rate as a function of the Strahler order, and, finally, the intravascular volume of the subtended tree to a vessel of a certain radius.  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973

Discussion
The relevance of generating networks of vessels to vascularize entire organs or even portions of tissues has turned into a valuable tool for research in drug delivery [42], and functional assessment of liver [43], heart [44], brain [23,45,46], among other applications. Thus, the ability of automatic algorithms to cover vascular territories with ensembles of vessels started to be tested in increasingly realistic applications and challenging problems in the field of biomedical engineering. CCO strategies emerged as a promising approach to address this problem [24], leading to variants, such as the DCCO approach [30] to enhance the capabilities of the methodology to deal with the construction of intricate vascular networks. The intrinsic limitation of the sequential CCO-based approaches is the limited number of vessels that the algorithm can generate in a reasonable timeframe. In this work, the main contribution was the proposition of a parallel algorithm consisting of three steps: (i) the generation of a baseline network using the DCCO approach; (ii) the partitioning of the vascular domain and the parallel application of the DCCO method to separately grow these networks and (iii) the merging step in which all networks are joined, and the complete network is adequately scaled to comply with the geometric constraints and the branching power law. . In lexicographic order, we report: combined cross-sectional lumen area, network connectivity, pressure and flow rate as a function of the Strahler order. In the bottom row, we report the subtended tree volume as a function of the vessel radius.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 In this work, we first proposed experiments to demonstrate that the PDCCO networks were in close agreement with DCCO networks in terms of geometry and function. We have seen that the vessel radius and length were in full agreement for both parallel and sequential algorithms. That was demonstrated for different numbers of segments in the baseline network and different subdomain shapes. These experiments were relevant to understanding the limitations in the setting of the parallel approach. Hence, we saw that even with baseline networks featuring a small number of segments, the error in the intravascular volume in PDCCO networks remained below 7% on average, while the radius and vessel length featured a similar distribution to DCCO networks. Regarding the aspect ratio of the partitions, the differences were also small.
Despite such a difference in vessel radius as a function of the depth in the vascular network (vessel generations), the computation of the pressure profile from proximal to distal locations does not change significantly between PDCCO and DCCO generated networks. The differences observed in the pressure drop were always below 1 mmHg on average. Pressure drops faster for PDCCO for deeper domains (i.e. A r > 1) but it is quite consistent for broad and cubic-like domains (i.e. A r ≤ 1). Thus, the partitioning of the domain should favour the creation of subdomains of regular dimensions or transversally larger sections (concerning the perfusion inlets) to avoid a bias of the pressure drop due to the parallelization.
At this point, we should emphasize that the embarrasingly parallel nature of the proposed algorithm does not support global optimality in the sense of conventional CCO-based algorithms. As a matter of fact, the original CCO approach does not yield a true optimal solution either. Indeed, in the sequential approach, the current state of the network is fixed, and only local bifurcation perturbation and whole-network diameter scaling are considered. In [47], the authors discuss what they call postoptimization, in which the state of the network, after connecting a new segment, is modified by altering all the bifurcation points. This enables the algorithm to reach a better degree of optimality. However, as noted by these authors, even if a more optimal network is achieved, the functional significance of these modifications ( pressure profiles) remains unaffected. The optimality of a certain network depends mainly on the placement of the first vessel generations [47]. Moreover, postoptimization adds significant computational burden, and is not practical for large problems. As in the original CCO method, in our approach, the optimality is mostly dominated through the first steps of the algorithm, when generating the baseline tree, which is accomplished through the conventional sequential algorithm (see §3.1).
The proposed PDCCO strategy can also be understood as a multiscale method for the generation of vascular networks. In fact, the baseline tree accounts for the macro-scale features of the vasculature, while the remainder of the network architecture responds to local micro-scale phenomena. In turn, the point at which the partitioning procedure is executed, after the baseline tree is generated, can be seen as the point at which we assume that the scale separation hypothesis holds. The multiscale nature of the approach is also evidenced in the experiments reported in §3.3 where the baseline tree topology was determined by the configuration of the cost functional, while the partitioned algorithm dealt with the micro-scale vascular features in a parallel fashion regardless of the macro-scale connectivity. Although there is no definite answer about which is the stage at which scale separation can be safely assumed, numerical experiments reported in this study evidenced a negligible departure from the optimal solution delivered by the sequential algorithm, both in terms of geometry and function. In turn, the merging procedure can be considered as a consistency step, in which the entire tree geometrically accommodates to the small-scale perturbations introduced in the micro-scale vascularization processes.
Concerning the scalability of the PDCCO method, specifically, the ability to efficiently generate extensive vascular networks, we have proposed a test to vascularize a domain with 100 000 vessels, using partitions with a different number of subdomains. In such a test, the construction of the baseline network is common to all the experiments, and the real potential of the parallel approach is seen in the phase where the baseline network is grown to the full merged network. The test allowed us to exploit the embarrassingly parallel nature of the proposed algorithm to reach a vascular network of 100 000 terminal vessels (that is 200 000 segments in the network). Removing the cost in the baseline network generation, the time spent in building the network scaled linearly with the number of partitions, and the geometric (intravascular volume and vessel radius profile) and functional (blood pressure profile) features of the resulting networks remained on an equal footing.
Finally, we have employed the proposed algorithm to build the entire vascular network in a prototypical human kidney. While previous efforts in this enterprise were reported in the literature [48], this has been the first attempt to build the entire renal vasculature from scratch to the level of interlobular arteries. Only by using the parallel approach was it possible to reach the level of royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 interlobular arteries. Although the morphometric validation of the so-constructed vascular model is not the focus of the present paper, it is possible to note the physiological morphometric (vessel length and radius) and functional ( pressure and flow rate) features in the model. In fact, note that it was not until only recently that reliable three-dimensional data about the rat renal vasculature was reported [49]. By using a multi-stage approach, the algorithmic generation was initiated with three segmental vessels, to vascularize the renal pelvis region in the first stage. The subsequent vascularization of the cortical tissue was developed until reaching a baseline network with 5000 terminal segments. This was done using the sequential DCCO algorithm. Then, the domain was partitioned, and the PDCCO algorithm was applied until reaching the goal of 100 000 terminal vessels (i.e. 200 000 segments overall). In this process, the first generations of segmental vessels were automatically placed in the renal core, and from there, the interlobar arteries, the arcuate vessels, and finally, the interlobular arteries. The perfusion pressure in the vascular segments generated with the proposed approach are in accordance with the measurements of intra-renal pressures assessed in vivo [50]. Even with an arborization involving tens of thousands of vascular ramifications from the renal artery, the perfusion pressure slightly changes. Substantial pressure drops occur only in the afferent arteriole segment, immediately before the entry into the glomerulus [50]. The construction of such a complex architectural arrangement of vessels was only possible due to the combination of the parallel and multi-stage approaches, paving the way towards an integrative analysis of the macrocirculation, microcirculation and cellular physiology for the description of whole-organ kidney function.
In [51], the authors characterized the morphology of the rat renal vasculature using micro-computer tomography images. They provided a comprehensive analysis of both arterial and venous networks, reporting different geometric network characteristics as a function of the Strahler order. In the present study, for the prototypical model of the kidney vasculature created using the proposed PDCCO, we also classified the geometric and functional entities associated with the model in terms of the Strahler order. The network behaviour yields an asymptotic growth of the cross-sectional lumen area as the Stralher order approaches zero. Flow rate drops more than linearly (log-scale), and the connectivity matrix reveals a coupling between small vessels and high Strahler order vessels. As a consequence, even if the pressure drops almost linearly with the Stralher order, there are small vessels which are exposed to high blood pressure. This may have an important role in the damage suffered by glomeruli (connected right distal to these smaller vessels) in hypertensive scenarios, even if the flow rate is small. In turn, the relation between the subtended intravascular volume tree and the vessel radius is consistent with the data reported in [51] (cf. Figure 11). In terms of model validation, there is scarcity of data in the field. One remarkable exception is [51], which studied the vasculature in the rat kidney. This speaks about the highly challenging enterprise that carrying such a study implies. In view of the lack of access to anatomical data, we focused on the description and validation of an efficient technique that enables the generation of massive networks. To do that, we rely on the already demonstrated capabilities of one of the most well-established approaches available in the literature, the CCO approach. Based on the CCO methodology, whose physiological and anatomical consistency was already demonstrated in numerous studies [24,25,[52][53][54][55][56], we proposed an algorithmic modification that allowed us to scale the number of vascular segments. Here, we proposed a strategy to remove an algorithmic barrier from CCO, and communicated evidence in favour of an efficient technique, which can perform equivalently to a standard and well-established one.
As a limitation of the proposed strategy, we can mention that the algorithm may suffer from a lack of convergence (rare but possible) when merging the networks in the case of nonlinear rheological blood behaviour. Another limitation we can mention is the degeneration of vascular networks in extremely degenerated subdomains. The latter can be mitigated by performing regular partitions of the domain to be vascularized.

Final remarks
This work presented a novel strategy to exploit the parallel generation of extensive vascular networks by using an existing DCCO algorithm. The parallel approach removed limitations inherent to the sequential algorithm, enabling the automatic generation of vascular networks containing hundreds of thousands of vascular segments.
One of the relevant contributions of the present work is the ability to create complex networks that integrate the macro-and micro-circulation realms, thus bridging the different scales featured by the CVS.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 8: 210973 Furthermore, we reported the in silico vascularization of a prototypical geometry of the human kidney for the first time. This result was achieved by integrating a multi-stage strategy and the parallel generation of vascular networks, starting at the renal artery and reaching the scale of interlobular vessels.