Modeling of Kidney Hemodynamics: Probability-Based Topology of an Arterial Network

Through regulation of the extracellular fluid volume, the kidneys provide important long-term regulation of blood pressure. At the level of the individual functional unit (the nephron), pressure and flow control involves two different mechanisms that both produce oscillations. The nephrons are arranged in a complex branching structure that delivers blood to each nephron and, at the same time, provides a basis for an interaction between adjacent nephrons. The functional consequences of this interaction are not understood, and at present it is not possible to address this question experimentally. We provide experimental data and a new modeling approach to clarify this problem. To resolve details of microvascular structure, we collected 3D data from more than 150 afferent arterioles in an optically cleared rat kidney. Using these results together with published micro-computed tomography (μCT) data we develop an algorithm for generating the renal arterial network. We then introduce a mathematical model describing blood flow dynamics and nephron to nephron interaction in the network. The model includes an implementation of electrical signal propagation along a vascular wall. Simulation results show that the renal arterial architecture plays an important role in maintaining adequate pressure levels and the self-sustained dynamics of nephrons.

By maintaining the volume and composition of the body fluids within narrow ranges, and by producing a set of hormones that affect the blood vessels, the kidneys provide important long-term regulation of blood pressure. Disturbances of kidney function can cause hypertension, a prevalent disease in modern societies. The kidneys protect their own function against short-term variations in blood pressure at the level of the individual unit (the nephron). In recent years, it has become clear that there is an interaction between nephrons, and that this interaction is mediated through the arterial network of the kidney. The renal vacular network has a complex topology, and at present there are no computational models of this topology, precluding a computational assessment of the consequences of Several recent modeling studies have investigated features of nephron-nephron interactions [20][21][22]. The complex renal arterial network in these studies was approximated by a simple bifurcating tree-like structure and included a group of 10-20 nephrons. Postnov et al. [20] showed that such a topology induces desynchronization between originally identical nephron models (coupling-induced inhomogeneity). Marsh et al. [21] demonstrated the presence of steady state, quasiperiodic and chaotic dynamics in an ensemble of cortical and medullary nephrons depending on the interaction strengths and the arterial blood pressure. Bayram et al. [22] suggested that several nephrons originating from the same small artery are more likely to be in an oscillatory state than a single nephron will. We recently demonstrated [23] that (i) neighboring nephrons at the same branch level were synchronized in-phase due to the vascular propagated electrical coupling, (ii) nephrons separated by one branch point tended to display a phase-shifted pattern due to hemodynamic coupling, and (iii) distantly located groups of nephrons showed asynchronous behavior. Marsh et al. [24] introduced heterogeneity into the nephro-arterial network, and observed similar dynamical patterns.
However, all of the studies cited used an overly simplified representation of the renal arterial network, and it is not clear to what extent the observed dynamics might be a consequence of the simplification. Further progress requires the use of more realistic models of the renal arterial network. In this study we focus on understanding how renal arterial structure affects blood flow dynamics. To address this question: (i) we perform experiments using optical tissue clearance to resolve details of microvascular structure not previously available; (ii) we develop an algorithm for generating a realistic model of the renal arterial network using the data obtained from our experiments and from micro-computed tomography study by Nordsletten et al [15]; (iii) we present a mathematical model that describes blood flow dynamics and nephron to nephron interaction in the renal nephro-arterial network including a new implementation of electrical signal propagation along a vascular wall; and (iv) we simulate how renal specific vascular structure can affect renal blood flow patterns and nephron-to-nephron interactions.

Algorithm for the vascular structure
We suggest a general algorithm to create an asymmetrical bifurcating tree based on experimental data, with modifications that are specific for the renal microcirculation: • An asymmetric bifurcating tree (ABT) with a gradual decrease of diameter, and with afferent arterioles present only at the end of the tree; • A kidney specific asymmetric bifurcating tree (KSABT) with a fraction of afferent arterioles branching directly from larger vessels.
In both cases the structure of the vascular networks is simulated using a probability-based bifurcating algorithm together with Murray's law [1]. Main features and assumptions of the suggested algorithm are described below: 1. Simulations start from a vessel with the initial diameter D initial .

A vessel always bifurcates into two daughter vessels;
3. For a daughter vessel the distribution of its diameter is a function of the parent vessel diameter; 4. The diameter of the second branching daughter vessel is calculated according to Murray's law [1], i.e. D d2 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 6. Vessels bifurcate until they reach a given diameter (parameter D stop ). Vessels with a diameter less than D stop are counted as afferent arterioles; 7. (Only for KSABT). Additional afferent arterioles are distributed over the whole vascular tree and appear with a probability according to the distribution of distances between neighboring afferent arterioles. We use a segmentation algorithm for this step. (Fig 1,  A more detailed description can be found in Materials and Methods. Statistical data on network topology The available data were collected by means of μCT with 20-μm and 4-μm resolution and published by Nordsletten et al. [15]. For the distribution of daughter vessel diameters as a function of the parent vessel diameter we use a cubic approximation for mean values and a linear approximation for standard deviations (Fig 2, left panel). For simplicity, it is assumed that branch level corresponds to a Strahler order. For a given parent vessel diameter, we use the cubic approximation to find the mean and standard deviation for the distribution of daughter vessel diameters. The diameter for the first daughter vessel is then drawn at random from this distribution. The diameter of the second daughter vessel is calculated so that it complies with Murray's law. For the length distribution as a function of the vessel diameter, a linear approximation is used for both mean and standard deviations for the best fit in the range of [20;100] μm (Fig 2, right panel). Based on these approximations we build a vascular tree with gradually decreasing diameter and with afferent arterioles located only at the end of the tree (ABT). Nordsletten et al. found only a small deviation (%1%) between a regression analysis performed on all daughter-parent diameter relationships and Murray's law [15], a finding that justifies our use of it. Accordingly, we designed our algorithm to ensure that the vascular tree it generates also follows Murray's law. Because vessel diameters at a bifurcation fulfill Murray's law, they are not statistically independent. This reflects the fact that knowing the diameters of two of the three vessels in a bifurcation, the diameter of the last vessel can be calculated from Murray's law. It should also be noted that the assumption that Strahler order corresponds to branching order may introduce some errors. Two daughter vessels need not have the same Strahler order, as this depends on the number of bifurcations that occur downstream from the respective daughter vessels [25][26][27]. Thus, the Strahler ordering scheme categorizes the vessels starting from the smallest vessels going upstream to the larger vessels. In principle it is therefore not possible to assign a Strahler order to a given vessel before the downstream branching process has come to an end. Since our algorithm "works" from the larger vessels towards the smaller ones, the Strahler ordering scheme is incompatible with the bifurcating algorithm, since it, in a sense, only can be applied retrospectively. However, the effect of this is expected only to be of minor significance. The data obtained by μCT provide important information about the dimensions of the larger renal vessels. However, due to the limited spatial resolution, important information on afferent arterioles is missing. Different studies show that (i) afferent arterioles (AA) can branch not only from small vessels, but also from larger arteries [28,29] and (ii) two, three, or four afferent arterioles can originate together from a single site on an artery [30]. These features make the renal vascular structure significantly different from the vascular structure of other organs [31,32].
We assume that the distances between the origins of AAs on the vascular tree can be described by an exponential distribution. This assumption is based on several observations: 1. AAs often originate in doublets or triplets even from large vessels. This implies a high probability of zero or close to zero distances between neighboring AAs; 2. Vessel segments without AAs branching from them are rarely seen. This implies a low probability for distances between AAs of several hundreds of microns.
The simplest possible distribution that satisfies observations 1 and 2 is an exponential distribution. Finally, we assumed that the distribution of the distance between neighboring AAs is independent of the diameter of the feeding vessel.
To provide experimental support for using an exponential distribution for the distances between the origins of afferent arterioles, we performed experimental studies to visualize renal vascular architecture using confocal microscopic imaging of chemically cleared renal tissue (see details in Materials and Methods) [16,17]. We obtained two 3D stacks of renal vascular structure, from which we measured the distances between the origins of neighboring AAs. An example of the observed structures is shown in Fig 3, right panel. The data set contains information on distances between origins of more than 150 afferent arterioles branching from the feeding vessels with diameters of 25-100 μm. An exponential function (red) provides a good fit to the data (blue) (Fig 3, left panel). The exponential distribution is used to generate a vascular tree specific to the kidney (KSABT) where afferent arterioles may branch off any vessel on the arterial structure.
Using the algorithm described above and the statistical data on renal vascular morphology, a renal vascular structure is generated in the form of a connectivity matrix that maps the connections between vessels as well as the vessel length and diameter. The matrix is used as the input for full-scale computer simulations.

Hemodynamic interaction
A hemodynamic interaction occurs when one nephron is stimulated by its TGF-mechanism to contract its afferent arteriole, causing the hydrostatic pressure to rise at the upstream branch node, increasing the blood flow to the second nephron. Half a period later when the increased blood flow activates the TGF-mechanism in the second nephron and causes its afferent arteriole to contract, the blood flow to this nephron is reduced, and the blood flow to the first nephron increases. Due to the time lag in the response, this type of interaction tends to induce antiphase synchronization between nephrons.
To describe the flow in each vessel of the vascular tree, one needs the hemodynamic resistance R hdr calculated from the Poiseuille relationship: Here L is the length of the vessel, η is the viscosity of blood. The viscosity depends on vessel diameter, and we adopt the expression for viscosity from Ref. [33], assuming an invariant hematocrit of 0.45: where η is in units of centipose and D is the inner diameter of the vessel in micrometers. The pressure variation in each node, derived from expressions for conservation of flow in each node (Fig 4), is given by: other : C hdr dP j dt Here C hdr is the vessel compliance. P pn j , P dn1 j , P dn2 j , and P j denote pressure in the parent node, the first daughter node, the second daughter node, and the current node, respectively. R hdr vj , R hdr dv1j , and R hdr dv2j are the hemodynamic resistance of the current vessel, the first daughter vessel and the second daughter vessel, respectively. Pressure in the root of the tree and blood flow to the nephrons F neph j serve as boundary conditions, where F neph j is calculated for each nephron as the difference between the pressure in the node giving rise to the afferent arteriole and the glomerular capillary pressure divided by the hemodynamic resistance of the afferent arteriole. The hemodynamic resistance is determined from the nephron model (Eq (16).

Electrical interaction
Activation of TGF in one nephron causes membrane depolarization and contraction of the vascular smooth muscle cells of the corresponding afferent arteriole. Because the cells of the renal vessels are coupled electrically by gap junctions [34], the membrane depolarization spreads, through electrotonic conduction, into the vascular smooth cells of the neighboring AAs causing them to contract [35]. Thus, the cells of the renal vasculature constitute a well-coupled syncytium in which the endothelial cells form the path of least resistance [36]. The deflection in the membrane potential generated by the macula densa decays exponentially with the distance between the site of contact between the macula densa and the AA [35,37]. The electrotonic spread between the vascular cells is fast, and it acts to adjust the TGF-mediated tubular pressure oscillations to attain a state of in-phase synchronization [18].  A fragment of the vascular structure. pn, dn1, dn2 and j denote the parent node, the first daughter node, the second daughter node, and the current node, respectively. Each vessel has its own hemodynamic resistance.
doi:10.1371/journal.pcbi.1004922.g004 Electrotonic propagation of the electric current originating from the macula densa through vascular segments is described by the following equations: For segment sites that are connected to nephrons: where I n stands for the current induced by the macula densa signal and transmitted to the vessel. V prev describes the voltage in the previous (downstream) adjacent unit of the same segment. V rest is the resting membrane potential, C u is the unit capacitance. G u and G g are the unit and gap junction conductance, respectively. For inner units of the segment: where V prev and V next denote potentials for downstream and upstream adjacent units, respectively. At the branch point (downstream of the vessel): where V 1,2 , and G g1,2 represent the voltages at the ends of connected vessel segments and gap junction conductances, respectively. If the branch point is upstream relative to the current segment then V prev should be used instead of V next . If the current segment has branch points from both sides then G g (V next − V) should be replaced with the corresponding conductances and voltages for the connected vessels. For the root segment G g1 ( should be replaced with specific values for G and resting values for V 1 and V 2 .

Nephron model
The model of the individual nephron consists of six coupled ordinary differential equations, each representing an essential physiological relation and a number of algebraic functions. A sketch of the main components of the nephron model is given in Fig 6. The six coupled differential equations are (see Material and Methods, [38] for details): The first equation of the model (eq 8) represents the pressure variations in the proximal tubule in terms of the in-and outgoing fluid flows. Here, F filt is the single-nephron glomerular filtration rate and C tub is the elastic compliance of the tubule. The flow into the loop of Henle is determined by the difference (P t − P d ) between the proximal and the distal tubular pressures and by the flow resistance R Henle . The reabsorption in the proximal tubule F reab is assumed to be constant. Eqs 9 and 10 describe the dynamics of the afferent arteriole. Here, r represents the inner radius of the vessel and v is its rate of change. ω is a measure of the mass relative to the elastic compliance of the arteriolar wall, P av denotes the average pressure in the arteriole, and P eq is the value of this pressure for which the arteriole is in equilibrium with its present radius and level of muscular activation. The expressions for F filt , P av and P eq involve a number of algebraic equations that must be solved along with the integration.
The remaining equations (eqs [11][12][13] in the single-nephron model describe the delay T in the TGF regulation. This delay arises both from the transit time through the loop of Henle and from the cascaded enzymatic processes between the macula densa cells and the smooth muscle cells that control the contraction of the afferent arteriole. Although the model is significantly simplified and does not contain a detailed description of all physiological mechanisms, its dynamical features include TGF oscillations, the response of afferent arteriole to increasing inlet pressure and, as a consequence, the autoregulation of efferent arteriolar blood flow (Fig 7).

Topological properties of the network
We used the algorithms described above to construct models of both the renal vascular tree and a simple asymmetric bifurcating tree.
Simple asymmetric bifurcating tree. Fig 8 shows an example of a simulation of the ABT structure based on the experimental data described by Nordsletten et al. [15]. The plots in the bottom panel show a high correlation of the simulated structure with the experimental data. Note, however, the divergence in the relation between daughter and parent vessel diameter in the simulated and the experimental data. As noted above, the experimental data are grouped according to Strahler order, and two daughter vessels will not necessarily belong to the same Strahler order. This is especially likely if the two daughter vessels have a large difference in diameter. In this case the larger daughter vessel will, on average, be expected to have a greater Strahler order than the one with the smaller diameter, since there is a greater likelihood for further downstream bifurcations in the branch coming from the larger daughter vessel. It is therefore not to be expected that the simulated and experimental distributions will be identical.
Bifurcating tree with exponentially distributed afferent arterioles. An example of a simulation of a KSABT structure is shown in Fig 9. Note that there is an even more pronounced difference in the dependence of daughter vessel diameter on parent vessel diameter compared to the experimental data [15]. The reason is that KSABT structure allows afferent arterioles to branch from any vessel (except the renal artery) and recalculates the vessel diameter after each bifurcation in accordance with Murray's law. This algorithm leads to segmentation of large vessels so that at the end of each segment the vessel divides into an afferent arteriole and a segment with a slightly reduced diameter. The distribution of the distances between afferent arterioles (Fig 9, bottom panel) is based on the experimental data obtained from optical clearing methods [16,17].
Notice that the number of afferent arterioles in the ABT and KSABT structures differ even though the parameters were similar in the two cases (see Table 1), and that in the KSABT structure the vessels feeding the afferent arterioles are, on average, wider than for the ABT structure where afferent arterioles appear only at the top of the tree. The larger diameter of the feeding vessels in the KSABT structure reduces the hemodymamic resistance of the vessels and, consequently, the pressure drop from the renal artery to afferent arterioles will be smaller when compared to that in the ABT structure. Dynamical properties of the network As discussed above, nephrons can interact with each other via hemodynamic and/or electrical coupling. These interactions can lead to different dynamical patterns, and can affect the physiological properties of the whole kidney. The vascular structure is important for the characteristics of both types of interactions. It affects the electrical properties of the vascular tree in the following ways: • The longer the distance between the sites where the electrical signal is generated and received, the weaker the signal. Hence, nephrons with long afferent arterioles and longer distances between each other show weaker interaction compared to nephrons connected via shorter distances (Fig 10) [39]; • At a branch point, the signal propagates better to the vessel with the larger diameter. This leads to weaker interaction between nephrons whose afferent arterioles originate from large vessels, while the strongest coupling is observed at the top of the tree, where the diameters of the daughter vessels are of similar scale.  To estimate the effect of vessel dimensions on the electrical coupling on a pair of nephrons, we performed simulations for a minimal branching structure consisting of one branch point and three vessels, a root vessel and two afferent arterioles. The lengths of the root vessel and the first afferent arteriole are fixed at 300 μm and 50 μm, respectively, and the diameter of the afferent arteriole is fixed at 20 μm. The length and diameter of the second afferent arteriole are varied from 50 to 250 μm and from 20 to 40 μm, respectively. The root vessel diameter is adjusted according to Murray's law. A characteristic feature of the renal vascular tree is that the pressure drop from the renal artery to the start of the afferent arteriole is small, whereas there is a relatively large pressure drop across the afferent arteriole [40]. This allows for efficient regulation of the glomerular filtration pressure through regulation of the resistance of the afferent arteriole.
Simulations on the two types of vascular structures show that the vascular structure with exponentially distributed AAs (KSABT) has much lower pressure drop compared to the ABT. Simulations for a small branch of 40 μm of initial diameter (Fig 11) with realistic hemodynamic resistances and root feeding pressure P root clearly show that in the case of KSABT (right panel) the tubular pressure P t of all nephrons is quite high and TGF oscillations can be observed. For the ABT, however, the root feeding pressure for the nephrons is low and outside the working range of the single nephron model that leads to negative values of P t (left panel). In the nephron model that we use the negative values of tubular pressure are observed when the root feeding pressure is low and cannot balance flows and pressures on the venous side. The higher pressures in the KSABT structure can be explained by two related factors: (i) the total number of branching levels is smaller in the KSABT than in the ABT network, and, (ii) the average diameter of feeding vessels is larger in KSABT, reducing the pressure drop between the renal artery and the afferent arterioles and providing higher pressure to the nephrons. The working pressure interval is depicted on the bottom panel of Fig 11. All nephrons in the KSABT work properly from 8 kPa in the root artery. This allows us to use a realistic pressure For the ABT structure, the nephrons are inactive for the whole range of the root feeding pressure, while for the KSABT structure all nephrons are active within a wide range of the root feeding pressure. Bottom panel: Dynamics of tubular pressure P t variations in all nephrons for pressure in the feeding vessel equal to 10 kPa (* 75 mmHg). This pressure level is insufficient for nephrons in the ABT structure to oscillate while for the KSABT structure the nephrons show oscillatory dynamics. Negative values of tubular pressure for the ABT structure are the consequence of insufficient pressure in the afferent arteriole to balance venous and filtration pressure and flow. Tubular pressure shown on the middle and bottom panels is color coded for each nephron according to the color markers on the top panel. (around 13-13.5 kPa) in the renal artery. In the case of ABT there is no oscillatory behavior even at "hypertensive" values of the pressure.
Autoregulation properties of the two types of structures differ. Nephrons in ABT are perfused at inadequate pressure to support autoregulation of total efferent flow or glomerular filtration rate (Fig 12, left). In contrast, the structure of KSABT provides nephrons with a higher pressure, supporting autoregulation at the level of individual nephrons and net flow (Fig 12,  right). A significant shift in nephron feed pressure in KSABT leads to prolonged and smoothed range of autoregulation for net flows.

Discussion
The systemic circulatory system plays a central role in supplying organs and tissues with oxygen and nutrients and removing metabolic waste products. On a large scale, the circulatory system can be viewed as a branching network where larger vessels on the arterial side branch into smaller and smaller vessels until the capillary level is reached. It is clear however that there are large differences in the structure of the vascular system among different tissues and organs, and that this difference relates to the specific functions of the tissue and organ. In other words, the vascular topology is optimized to subserve the specific needs of a given organ or tissue. In skeletal muscle, the main need is to allow for large variations in total blood flow so as to match the supply with demand during muscle work. In other organs, like the brain, there is no need for large variations in overall flow, but there is a need to redirect flow to local areas where the metabolic activity is high. In the kidney, the critical aspect of the circulation is not only to supply the tissues with oxygen and nutrients, but rather to support the filtration of plasma at the glomerulus of each nephron. For example, in the human kidneys an amount of fluid corresponding to approximately 3 times the body weight is filtered every 24 hours. This difference in function poses restrictions on the design of the vascular system in the specific organs and tissues. Many diseases, e.g. hypertension and diabetes, affect the vascular system, and the morbidity and mortality of these diseases are to a large extent associated with the diseases' effects on the vascular system. The effect of the disease in a given organ is not only the result of the disease process itself, but also depends on the specific topology of the vascular system in the organ [10].
In recent years several new imaging modalities, e.g. μCT, and confocal microscopy, have made it possible to obtain images of the entire microcirculation in a given organ [9,15]. This has paved the way for a better understanding of the complex interrelationship between organ function and vascular topology. To obtain a thorough understanding of the role of vascular topology in organ function it is necessary to combine experimental and modeling approaches, since many questions cannot be addressed experimentally. In this connection, a major task is to create realistic models of the vascular tree in a given organ as a basis for further model studies.
In this study we present new experimental findings of the 3D arterial and afferent arteriolar structure and we develop a probability-based algorithm for generating such structures from the data.To study the role of kidney specific vascular structure we also introduce a mathematical model that describes blood flow dynamics and nephron-to-nephron interactions in the arterionephron network. The simulation results indicate that the arterial structure in the kidney minimizes the pressure drop between the main renal artery and glomeruli, and it affects nephronnephron interactions.
One of the main results is that the distances between individual AAs are exponentially distributed. This conclusion is the result of the analysis of 3D data sets with more than 150 afferent arterioles obtained with the optical clearance method. This result differs from that in most other organs where the interbranch distance seems to follow a lognormal distribution [31,32]. Exponential distributions are typically seen in Poisson processes where the time or distance between events are independent of each other. The genetic/molecular mechanisms that underlies such a branching pattern are unknown, and an interesting subject for further theoretical work. Statistical distributions derived from our data and from the literature [15] form the basis for our structure generating algorithm.
The microcirculation of the kidney is unusual because it has 2 capillary networks, the glomerular capillaries and the peritubular capillaries, connected through the efferent arteriole. Compared to capillary networks in other organs, the glomerular capillaries operate with high intravascular pressure, a condition required to to drive the filtration of fluid across the capillary wall into the lumen of the proximal tubule. The high filtration rate is a prerequisite for the ability of the kidneys to regulate the volume and composition of the extracellular fluid. When we used a simple bifurcating tree (the ABT algorithm), it was apparent that with the vessel dimensions reported by Norsletten et al. [15], the pressure drop from the renal artery to the glomerular capillaries exceeds the value found experimentally [40]. This is not surprising, since in the ABT the AAs only appear at the terminal branch points of the tree, an assumption that maximizes the hemodynamic resistance between the renal artery and the glomerulus. In contrast, when afferent arterioles are distributed exponentially across the vascular tree and allowed to branch from any arterial segment, the resulting pressure in the glomerular capillaries is significantly greater, and in a range that is compatible with normal nephron function.
It is of interest to note that only the KSABT algorithm resulted in a realistic number of nephrons for the whole kidney (Table 1), i.e. around 33,000 [41,42]. Because the ABT only has glomeruli at the terminal branches, the number of nephrons becomes much lower than in the KSABT, where AAs, and thus nephrons, also originates from the larger vessels. For the ABT to give a realistic number of nephrons, it will be necessary to start it with a vessel of an unrealistically large diameter.
The topology of the arterial network influences not only the pressure drop along the network, but also electrically mediated nephron-nephron interactions. As shown in Fig 10 signal propagation depends on the dimensions of the branching vessels. A signal from the macula densa is conducted through the nephron's afferent arteriole and, at the next bifurcation, propagates preferentially to the largest diameter branch. A nephron whose AA originates from a larger vessel will therefore be less efficient in synchronizing with its neighbors than will a nephron at the distal end of the arterial tree, where the vessels tend to be of a similar scale. Previous studies on nephron-nephron interactions have used a simple bifurcating tree when investigating the network properties of the renal circulation [23,24]. Such simplifications, by ignoring significant structural asymmetries, could therefore overestimate synchronization. Another important feature is that asymmetrical connections of endothelial cells can lead to anisotropic propagation of electrical signal between two afferent arterioles and can affect nephron-tonephron interaction.
Although our algorithm results in architecturally realistic arterial networks for the renal circulation, the tubular pressure distribution between nephrons is exaggerated, cf. Fig 11. Experimental data has shown that the proximal tubular pressure is quite uniform in nephrons on the surface of the kidney [43]. This suggests the presence of additional processes that acts to homogenize the pressure distribution between nephrons. One possibility is that the length, and thereby the resistance, of the AA may vary according to the branch site, i.e. if the AA branches from a larger vessel with a high intravascular pressure it may be longer compared to a similar vessel that from a smaller vessel. At present, there is no data on renal vascular morphology to test this possibility, and it has therefore not been included in the present algorithm. Another possibility is that the individual vessels actively adjust their radius to compensate for differences in the feeding pressures. The AA has a myogenic response [44] that allows it to adjust its radius in response to changes in the intravascular pressure. An increase in pressure will provoke vasoconstriction, increasing the hemodynamic resistance of vessel, reducing glomerular pressure towards the control value. This mechanism plays an important role in renal autoregulation of blood flow [45] and is implemented in our model in a simplified way.
An anatomical feature not explicitly included in the present algorithm is the presence of triplets and quadruplets at the branch sites, as reported in [30]. Such triplets or quadruplets of AAs mostly occur at the top of the tree, since the small arteries will always split onto two AAs. We do not explicitly model the population of juxtamedullary nephrons in the present work. However, it is known that this subpopulation of nephrons predominantly derive their afferent arterioles from the larger vessels in the arterial network, e.g. the arcuate and subarcuate arteries [28]. However, the available data does not allow specific modeling of this subpopulation.
In the present work we have chosen to use a simple nephron model [38]. The model, detailed below, has primitive representations of TGF and the afferent arteriole that cause a single combined effect on afferent arteriolar diameter in response to arterial pressure change. When several copies are combined in a network configuration the model generates interactions [20,21,23] similar to those produced with more detailed representations of TGF and the afferent arteriole [24]. Because of the complexity of the network a minimal model such as the one we use here best serves the purposes of the study by providing the clearest opportunity for understanding the contributions of the vascular structure. To test stability of the model we ran simulations for different values of parameters C hdr (0.3-5nL/kPa), β(0.4-0.67), α (12)(13)(14)(15)(16)(17)(18)(19)(20). We found minor quantitative changes but the dynamics of the model with ABT and KSABT structures remained qualitatively unchanged.
In conclusion, we have demonstrated that the renal vasculature has specific characteristics that were not taken into account in previous modeling studies. We have developed a new algorithm for generating a renal arterial network using experimental data on the length and radius distribution of renal vessels. The resultant network topology closely resembles that found in normal rats. We have found that in the contrast to the simpler known algorithm, the kidney specific vascular tree (based on our data and from the literature) contained a realistic number of nephrons, displayed adequate pressure levels in the nephrons, and tended to weaken interactions between nephrons whose afferent arterioles originated from larger vessels compared to the nephrons on top of the tree.

Materials and Methods Experiment
All experimental protocols were approved by the Danish National Animal Experiments Inspectorate and were conducted in accordance with guidelines of the American Physiological Society. A 12 week old Sprague-Dawley rat was anesthetized (intraperitoneal injection of pentobarbital), a catheter was inserted into the aorta, and the vena cava was opened. The blood was removed by perfusing the animal with PBS with nifedipine and HEPES to prevent blood clotting. Thereafter the animal was perfused with biotin (2 mg/ml) that binds to endothelial cells, washed with PBS and perfused with Alexa-647 Streptavidin (20 mg/ml). The animal was then fixed with 4% paraformaldehyde in PBS using a pressure of 90 mmHg. The kidneys were removed and a section of one of the kidneys was clarified according to Erturk et al. [17], using dehydration steps in THF and final clearing in DBE [16]. Tissue dehydrated by this method shrinks approximately 20% in each dimension. The kidney slice was scanned with a 633 nm laser in a Zeiss LSM 710 equipped with a 10X/0.3 NA objective. Two 3D stacks were recorded, covering 850 x 850 x 497 μm and 850 x 850 x 1018 μm. The second image stack was obtained with spectral imaging to obtain both green autofluorescence from tubules and the Alexa-633 signal from both vessels and tubuli. The combined image showing tubuli in one channel and vessels and tubuli in the other was used as input for the segmentation protocol. The dual channel image was then automatically segmented in the open source image processing package "FiJi", using pixel-based segmentation plugin "Trainable Weka Segmentation". After segmentation the 150 images in the stack were manually corrected in FiJi for missing features and holes in vessels. The 3D representations (Figs 3 and 13) were performed with open platform for bioimage informatics "Icy".

Computational implementation
To build a vascular structure using the present algorithm four distributions are needed: G ddp (daughter-parent diameter distribution), G vlvd (vessel length-vessel diameter distribution), E aad (distribution of distances between two neighboring AAs) and G aad (distribution for afferent arterioles diameter), where G ddp , G vlvd and G aad are assumed to be Gaussian distributions and E aad is an exponential distribution. All distributions were obtained from experimental data. Two parameters are used to control the structure size: D initial describes the diameter of the first vessel and D stop is the vessel diameter where bifurcations stop. 5. If the next segment diameter becomes smaller than D stop , stop branching of this vessel.

Description of electrical coupling
We simplified the model described in Ref. [46] with some additional assumptions: 1. All related processes in the endothelium are much faster than the TGF feedback and, thus, all membrane potentials are instantly adjusted to the signal produced by the macula densa; 2. We use a linear representation of the whole-cell conductance. This allows us to calculate a matrix of connectivity coefficients between nephrons, instead of running an additional integration loop for each integration step of the whole model; 3. At the branch points, cells from each vessel split in numbers according to the diameters of the other two vessels. The conduction of such a connection is proportional to the number of connected cells from both sides; 4. We neglect latitudinal variability of endothelial membrane potentials, but still take into account longitudinal discreteness of the vascular bed. This allows us to build a computationally simple, but physiologically relevant, model for the propagation of the electrical signals generated by the nephrons. Fig 15 represents the main notations of the coupling. Each endothelial cell is described by the whole-cell conductance G c which can be linear or nonlinear, and instantaneous (very fast activated) or inertial (governed by the dynamics of gating variables). We use the simplest representation in the form of an instantaneous and linear whole-cell conductance G c (assumptions 1 and 2) being connected in parallel with the whole-cell capacitance C c and coupled with neighbors by means of gap junction conductances G gj .
There is considerable variability with regard to the reported morphological and electrical parameters of endothelial cells. For the present work we chose parameters to be in accordance with those in Ref. [46]. We assume that a vessel segment of diameter D contains N = πD/W c endothelial cells in its cross-section, where W c % 5μm is a typical width of an individual cell. Similarly, we assume that the length L of the vessel segment is composed of M = L/L c endothelial cells, where L c % 50μm is a typical effective length of an endothelial cell. The length of an endothelial cell is around 100 μm. However, due to overlap of the ECs in the vessel wall, the effective length is approximately half the the real length of the cell [47].
We ignore possible spatial inhomogeneity of voltage distribution in a vessel cross-section. Thus, we describe each unit (piece of vessel with length L c ) with its unit conductance G u = NG c , unit capacitance C u = NC c , and total gap junction conductance G g = NG gj . A vessel segment is represented by a chain of M such coupled units. At the branch (bifurcation) point three vessel segments are connected to each other. At this point we assume a triangle-shaped coupling geometry as shown in Fig 15. To calculate G g1 , G g2 , G g3 we use assumption 3 listed above. Cells from the i-th segment connected to the j-th and k-th segments are split in two parts according to circumference (or equivalently-diameters) of the jth and k-th segments. So, N i = N j + N k , where N j /N k = D j /D k and N is the number of cells which is always an integer. However, this rule gives different number of cells for different vessels at the same branch point. For example, the number of cells N j−k from j-th to k-th segment will differ from the number of cells N k−j from the k-th to j-th segment. The existence and value of this mismatch can not be predicted in advance. Thus, we assume that the number of gap junctions between two segments is equal to the arithmetic average between the number of cells in the connected vessels. In this way we have some intermediate, but symmetric value for inter-segment coupled cells: N jk = N kj = (N j−k + N k−j )/2. The conductance of such coupling is G g jk = G g kj = N jk G gj .
At the root of the tree the same rule for the conductance is applied. However, we assume that the diameter of vessels downstream from the root is large enough, so that each endothelial cell from the root segment will be connected to two endothelial cells from downstream vessels. This gives us G g root = 2N i G gj for each connected vessel, where N i is the number of cells in the circumference of the root segment.

Description of nephron model
Physiological justification for all equations and expressions are given in Ref [38]. Here we focus on computational implementation of the model. P t , r, v r , X 1 , X 2 , X 3 denote dynamical variables of the system. F Hen , F Filt , P g , P av , P eq , P el , P act , C, R a , R, C e are functions that have physiological meaning and are needed to calculate derivatives. A, B, C, D are coefficients of the 3rd order polynomial equation for the efferent plasma protein concentration, C e .
P v , P d , R a0 , R e , R Hen , C tub , H a , F reab , F Hen0 , C a , a, b, ω, K, β, C max , C min , C eq are constants. P a , T, α are control parameters. Full list of parameters can be found in Table 2.
The functions that are used in the model are described below. Flow in the loop of Henle: Preglomerular resistance We solve the third-order equation AC 3 e þ BC 2 e þ CC e þ D ¼ 0 to find C e , the efferent arteriolar plasma protein concentration. For appropriate parameter valse the equation has a single, positive solution. This solution is determined numerically for each integration step. A, B, C, D are found from the model equations: C ¼ P t À P v þ R Â a Â C a Â ð1 À H a Þ þ R Â ðP t À P a Þ Â H a ð19Þ D ¼ ðP t À P a Þ Â R Â C a Â ð1 À H a Þ ð 20Þ Functions of plasma protein concentrations: P av ¼ 1 2 Â P a À ðP a À P g Þ Â b Â R a0 R a þ P g ð22Þ The tubuloglomerular feedback function C: The depolarisation of the cells in the afferent arteriole due to the TGF signal from the macula densa is assumed to be directly proportional to C. Equilibrium pressure: P el ¼ 1:6 Â ðr À 1Þ þ 2:4 Â e ð10ÂðrÀ1:4ÞÞ ð25Þ P act ¼ 4:7 1 þ e ð13Âð0:4ÀrÞÞ þ 7:2 Â r þ 6:3 ð26Þ Note that the parameter P a in the single nephron-model becomes a variable in the nephro-arterial network.