Abstract
Tubuloglomerular feedback (TGF) and the myogenic mechanism combine in each nephron to regulate blood flow and glomerular filtration rate. Both mechanisms are nonlinear, generate selfsustained oscillations, and interact as their signals converge on arteriolar smooth muscle, forming a regulatory ensemble. Ensembles may synchronize. Smooth muscle cells in the ensemble depolarize periodically, generating electrical signals that propagate along the vascular network. We developed a mathematical model of a nephronvascular network, with 16 versions of a single nephron model containing representations of both mechanisms in the regulatory ensemble, to examine the effects of network structure on nephron synchronization. Symmetry, as a property of a network, facilitates synchronization. Nephrons received blood from a symmetric electrically conductive vascular tree. Symmetry was created by using identical nephron models at each of the 16 sites and symmetry breaking by varying nephron length. The symmetric model achieved synchronization of all elements in the network. As little as 1% variation in nephron length caused extensive desynchronization, although synchronization was maintained in small nephron clusters. Inphase synchronization predominated among nephrons separated by one or three vascular nodes and antiphase synchronization for five or seven nodes of separation. Nephron dynamics were irregular and contained lowfrequency fluctuations. Results are consistent with simultaneous blood flow measurements in multiple nephrons. An interaction between electrical signals propagated through the network to cause synchronization; variation in vascular pressure at vessel bifurcations was a principal cause of desynchronization. The results suggest that the vasculature supplies blood to nephrons but also engages in robust information transfer.
 renal autoregulation
 tubuloglomerular feedback
 myogenic mechanism
 mathematical models
 network theory
tubuloglomerular feedback (TGF) and the myogenic mechanism combine to regulate blood flow and glomerular filtration rate (GFR) in each nephron. Signals from the two mechanisms converge on the contractile mechanism in arteriolar smooth muscle cells, creating an obligatory interaction. Each mechanism is nonlinear and generates selfsustained oscillations of arteriolar vascular resistance and nephron blood flow in normotensive animals (18, 26, 50); synchronization can arise from the interaction and result in frequency and amplitude modulation (29, 31, 49). Frequency modulation is an effective signaling mechanism, and amplitude modulation permits the forcing from blood pressure fluctuations to reach the macula densa, providing a coordinated autoregulatory response. The myogenic:TGF frequency ratio takes on discrete integer values in experimental data, an adjustment that most likely arises from synchronization of these two oscillators (44).
TGF activation at the macula densa is transmitted to cells of the afferent arteriole, inducing depolarization and Ca^{2+} influx from extracellular fluid (32). The depolarization, and the vasoconstriction it causes, are propagated electrotonically to other afferent arterioles branching from a common feeder artery. TGF and myogenicactivated electrotonic signaling from several nephrons synchronizes their rhythms in normotensive animals (12, 44).
The frequency of TGF oscillation is in the range of [30,50] mHz. There is a slower change in proximal tubular pressure and renal blood flow, but no separate mechanism has been identified as its cause (36, 42).
Using laser speckle contrast imaging to provide simultaneous measurement of relative blood flow changes in many nephrons on the surface of normal rat kidneys, we found similar frequencies within clusters of two or three nephrons (19). Larger clusters form and dissolve, lasting only some fraction of the 30 to 40min observation periods. Injections of vasoactive agents alter blood flow in the whole kidney and in all individual nephrons, always in the same direction with a particular agent. The spontaneous fluctuations in blood flow reflect both vasoconstriction and vasodilatation occurring simultaneously in different clusters. We conclude that the fluctuations in the synchronization arise from an intrarenal source. Our hypothesis in the current work is that the network of interacting nephrons generates the variability in synchronization.
Symmetric networks can produce stable synchronized solutions (10). A system is symmetric if exchange of elements within the system leaves the system unchanged. A group of identical interacting nephrons, connected in a branching vascular network with identical vascular segments of a given order, is symmetric by this definition. Symmetry is broken if one or more of the nephrons, or one or more of the vascular segments, or both, has properties different from the others.
Interactions among oscillating nonlinear systems provide a basis for synchronization, and synchronization a basis for the formation of a selforganizing system. A question that arises in theoretical work on nonlinear networks is the extent to which synchronization and selforganization persist despite broken symmetry (10). Such a possibility exists, but its realization depends on the details of each system. From the results with laser speckle imaging, it would appear that fields of synchronized nephrons are limited in extent and variable with time (19). Our hypothesis is that these properties appear because of broken symmetry within the renal vascular network. We further propose that some synchrony persists and provides the basis for the formation of a selforganizing system of nephrons.
To test these hypotheses we have developed a model of a network of nephrons, interacting electrotonically and also through local changes in blood pressure induced by the dynamics of the nephrons. In animal experiments we found that proximal tubular pressure in animals with chronic hypertension has irregular fluctuations with broadband spectral power, rather than the regular oscillations found in normotensive rats (12, 14, 48, 49). Two potential sources of this change have been identified: an increase in TGF gain (8, 27) and an increase in the electrotonic conductance between nephrons (6, 47). Simulations with nephron models show that each independently supports the change in dynamics observed in hypertension (21, 30). The network model will permit us to evaluate the relative importance of each and the interaction between them.
METHODS
We report on the development and exercise of a mathematical model of nephrons in a vascular network. The nephron model has been presented and analyzed (21, 29, 31, 32); the vascular network is new.
Model Components
Vascular model.
The structure is based on an analysis of images from microcomputed tomography of a single rat kidney. Nordsletten et al. (34) analyzed data collected by GarciaSanz et al. (9). The analysis used a series of algorithms that produced length, radius, and branching data on arterial and venous vessels, ordered according to dimensions, and a connection matrix detailing the statistics of branching from each vessel order. The overall pattern was of a bifurcating network, each vessel giving rise to two branches on average. This description applied at all levels of both arterial and venous vessels.
Figure 1 shows our realization of the results of Nordsletten et al. (34). Blood enters the network through a single feeder artery and then flows through a series of bifurcations at nodes, ending at paired afferent arterioles. The structure is intended to be symmetric; broken symmetry is known to be important in determining the dynamic behavior of networks (10) and the symmetric structure provides a reference state from which to evaluate the effects of specific forms of symmetry breaking. The structure in Fig. 1 is intended to represent a threedimensional network, with no particular spatial orientation. Any combination of the model's nephrons can represent those on the surface of the kidney or below it.
Because the results reported by Nordsletten et al. (34) are from a single rat one cannot judge how representative their numerical values are, and so we used their results only as guidelines for selecting dimensions of the arterial network. The fundamental experimental observation of the phenomena we study is the set of measurements of proximal tubular hydrostatic pressure. Values of segmental arterial length and radius were adjusted until the individual nephron models predicted tubular pressure oscillations comparable in amplitude to experimental values, with the radii of the vascular branches following Murray's law (33). Table 1 lists the values selected.
Tubular pressure measurements are available from several laboratories. Mean proximal tubular pressure typically averages 8–12 mmHg, and the oscillations average ±1–2 mmHg. The model predicts a mean value of 9.8 mmHg and oscillations of ±1.3 mmHg when all nephrons in the network are identical. Both the average pressure and the magnitude of the tubular pressure oscillations may change when nephrons are not identical, depending on how they are made to differ, but often they are not different from the values cited for simulations with identical nephrons.
The electrical circuit of the vascular structure was modeled as a network of resistors, the resistance increasing with segmental length, with electrical nodes at vascular branch points. Applying Kirchhoff's Law, ∑I_{i} = 0, where the I are electrical currents flowing into a node. Referring to Fig. 1 and substituting the appropriate conductances and electrical potential differences yields the following equations for the series of nodes between nephrons 1 and 2, and the origin of the vascular tree, as follows:
Hydrostatic pressure in the vascular tree is calculated using the dimensions of the vessels and assuming Poiseuille flow. Blood viscosity varies as a function of vessel radius, an effect not accounted for in the HagenPoiseuille formula, but renal vascular dimensions are not known with sufficient precision to calculate these viscosity changes accurately. We selected a set of lengths and radii that satisfies Murray's law (33) and that provides hydraulic resistance to enable nephron models to generate measured nephron pressures, flows, and concentrations. We intend this approximation scheme to subsume the viscosity effect.
Tubular model.
For each nephron, pressure, volumetric flow rate, and NaCl concentration were simulated as functions of time and distance in a mechanically compliant tubule that reabsorbs water and NaCl. The model contains segments corresponding to the proximal tubule, descending limb of Henle's loop, ascending thick limb of Henle's loop, and early distal tubule, each segment with its own parameter values. The boundary conditions were glomerular filtration rate as the initial proximal tubular flow; hydrostatic pressure at the end of the early distal tubule as a nonlinear function of the flow rate in that segment, reflecting the compliance of nephron segments beyond the early distal tubule; and initial NaCl concentration in tubular fluid at the beginning of the descending limb of Henle's loop. NaCl is the only solute. The concentration of NaCl remains invariant along the length of the proximal tubule (17, 29).
GFR and TGF models.
Glomerular filtration rate was predicted using an equation derived by Deen et al. (7), based on conservation of protein mass in glomerular capillaries; local determinants of filtration rate are glomerular capillary hydrostatic pressure, early proximal tubule hydrostatic pressure, colloid osmotic pressure, and glomerular capillary filtration coefficient. The action of TGF was predicted as a function of NaCl concentration in tubular fluid at the level of the macula densa, with an algebraic equation used to fit data from tubular perfusion experiments. TGF action was used to modify Ca^{2+} conductance in vascular smooth muscle cells of the afferent arteriole.
Myogenic mechanism model.
We modeled afferent arteriole cells with an adaption of the GonzalezFernandez and Ermentrout (11) model of vascular smooth muscle cells, originally derived for simulating cerebral arteriolar vasomotion. The six dependent variables in each segment were as follows: fraction of open K^{+} channels, intracellular Ca^{2+} concentration, electrical potential difference across the plasma membrane, myosin light chain phosphorylation, length of the contractile mechanism, and length of a set of parallel elastic components, all as functions of time, in a single cell. The open probabilities of both the K^{+} and the Ca^{2+} channel were membrane potential dependent, the open K^{+} channel probability was also dependent on the intracellular Ca^{2+} concentration, and the interaction of the membrane currents of these two ions produced an instability in membrane potential that gave rise to a limit cycle oscillation of the electrical potential and the ionic currents. The oscillation of intracellular Ca^{2+} concentration serves as a forcing that causes oscillations in myosin light chain phosphorylation, in the contractile mechanism, and in the length of the elastic elements. The strength of contraction in the model is dependent on the length of the contractile mechanism, enabling the model to respond to changes in wall tension caused by variation in arteriolar pressure.
TGF activation initiates contraction of the afferent arteriole (32) at its point of contact with the juxtaglomerular apparatus; the vasoconstriction declines with distance from the point of initial contraction (47). We divided each of the 16 afferent arterioles in Fig. 1 into 2 segments and used the vascular smooth muscle cell model to simulate each of the resulting 32 segments. The two segments in each arteriole were intended as a twopoint approximation to a continuous length dependence of the strength of the TGF signal. As a part of the twopoint approximation, the effect of the TGF signal in the segment farther from the glomerulus was 0.33 of the effect in the nearer segment; the two segments were made to interact through an electrical conductance between them. The details of these model components, including the parameter values, have been published in a number of studies, most recently in Marsh et al. (32); a causal model diagram is in Laugesen et al. (21) and is reproduced here as Fig. 2.
This single nephron model, including the TGF, glomerular filtration, and myogenic components, was used without modification from the models that have been published, with a single exception. In the earlier work, a constant vascular resistance term accounted for the pressure drop along upstream renal arterial segments. That constant resistance term has been deleted in the present work because the pressure drops across the upstream vessels are now calculated explicitly.
Numerical Methods
Twelve nonlinear ordinary differential equations (ODEs) describe the myogenic mechanism in the afferent arterioles of each of the 16 nephrons, forming a set of 192 equations. We solved this set numerically with a routine that provides a variable step size, generates numerical values for the Jacobian, uses backward differentiation, and iterates at each time step until error limits are satisfied. The absolute error limits were set at 1 × 10^{−8} for all of the ODEs with the exception of those representing the lengths of the contractile mechanism and the elastic elements. These length variables were at least an order of magnitude smaller than the other variables, and were assigned an absolute error limit of 1 × 10^{−9}. The solutions for the electrical potentials of the arteriolar cells cross 0 mV (11, 29). The absolute error limit requirements for these variables were therefore eliminated when their values were in the interval [−0.1 mV, 0.1 mV]. Relative error limits for each variable were set at 5 × 10^{−8}. The ODE solver routine, known as DVODE, was developed at the Lawrence Livermore National Laboratories and is available at https:https://computation.llnl.gov/casc/software.html.
A set of nodes like those of the network shown in Fig. 1, with conductances among them, is a conductor bridge. Electrical potentials imposed at the periphery of the set uniquely determine the potentials on each of the internal nodes. These internal node voltages are given by v_{unknown} = −D^{−1}Cv_{known}, where C is a 17 × 15 matrix and D is a symmetric 15 × 15 matrix; i designates columns and j designates rows in both.
The C_{i,j} are the conductances between nodes on the periphery, i, and those in the interior, j; the offdiagonal terms, D_{i,j}, are the conductances among the interior nodes, while the diagonal terms, D_{i,i} are the negative of the sum over all of the conductances in row j = i of C and D; v_{known} is the vector of known electrical potentials; and v_{unknown} is the vector of unknown electrical potentials. The 17 elements of v_{known} are the electrical potentials in each of the 16 afferent arterioles and electrical ground. The elements of v_{unknown} are the electrical potentials of the 15 vascular nodes as defined in Eqs. 1–4, where nodes 1–8 are defined by equations similar to Eq. 1, nodes 9–12 by equations similar to Eq. 2, and nodes 13 and 14 by equations similar to Eq. 3 and node 15 by Eq. 4. The quantity D^{−1}C is only dependent on network topology and conductances, which are constant throughout the simulation and therefore require the matrix multiplication to be performed only once. The v_{known} are calculated with each call to the ODE solver routine, and the v_{unknown} are updated by calculating D^{−1}Cv_{known}. The ODE solver iterates the sequence until the error limits are satisfied.
The partial differential equations describing tubule pressure, flow, and NaCl concentration were approximated as a set of finite difference equations using a centered difference scheme (46). Nephrons were discretized into 195 segments, each 100μm long. Solutions were iterated to convergence for the vectors of the dependent variables independently in each nephron. The convergence criteria were 1 × 10^{−3} for the combined vector of pressures and flows and 1 × 10^{−5} for the vector of NaCl concentrations. The solution routine for each nephron was assigned to its own processor, and the processors ran concurrently.
The main program controls all initialization procedures, advances the solution for all components at constant time intervals of 1 ms, calculates blood pressures and flows in the vascular network, and manages message passing with the nephron routines. Operator splitting was used, as follows: using a time step for the whole system of Δt s, at time t the solver for the ordinary differential algebraic equations was run over the time interval [t, t + Δt/2], producing updated values for the dependent variables of the myogenic mechanism, the network electrical potentials, and the vascular hydraulic resistances. These new values were supplied to the nephron routines, which were then executed over the time interval [t, t + Δt]. New values for nephron blood flow, TGF signal to the afferent arteriole, and glomerular capillary pressure were supplied to the ODE solver, and that routine was run over the time interval [t + Δt/2, t + Δt], to complete the work in the interval [ t, t + Δt]. Unless otherwise noted, runs simulated a 1,500s interval. The first 600 s were allowed for the relaxation from initial conditions and were discarded; results of the last 900 s were used for data analysis and in the figures. The real time required to complete a 1,500s simulation was typically ∼1,500 s. The time required for the simulation varied primarily because different conditions affected the speed with which the ODE solver could satisfy the error limits.
A machine with 48 microprocessors, maintained in the Center for Computation and Visualization at Brown University, was used to generate all results we report. Preliminary calculations and program development were carried out at the High Performance Computing Center of the University of Southern California.
Data analysis.
Synchronization among elements of a network is an important outcome of interactions within the network (1, 38). In the renal vascular network, these interactions arise among nephrons because TGF and the myogenic mechanism develop selfsustained oscillations. Over time periods of 30 min or longer the oscillations are made irregular by blood pressure fluctuations (13, 28) and by intrinsic factors, like deterministic chaos (21, 30, 48, 49, 51). We evaluated phase locking between pairs of nephrons to estimate the extent of synchronization. The phase angle of an oscillating signal is a linear function of its frequency. Two oscillators are said to be phase locked over some time interval if the difference between their phase angles remains constant over that interval. The presence of phase locking signals means that the two frequencies are similar, which is the standard definition of synchronization (1, 38). The dispersion of the phase angle difference over time is therefore a measure of the extent of synchronization, and its mean value characterizes the type of synchronization, as discussed below. We used techniques of directional statistics (4) to provide these estimates from simulated nephron blood flows.
Simulations provided time series of nephron blood flow over 900 s, sampled at 4/s, generating 3,600 values. Mean values were subtracted from each individual value of a given pair of time series, yielding two realvalued processes (45), x(k) and y(k), 1 ≤ k ≤ 3,600. The instantaneous phase of a realvalued discrete process, x(k), can be defined as the argument of the corresponding analytic signal,
We organized nephronpair interactions according to the number of vascular nodes separating them, and used one way analysis of variance and the NewmanKeuls post hoc test to assess significance of differences in

RESULTS
All simulations were performed with time invariant arterial pressure of 100 mmHg at the entrance to the feeder artery; the dynamics we report arise from actions of network elements.
Symmetric Network
Figure 3 shows the results of using nephrons with identical properties at each of the 16 nodes in Fig. 1. The nephrons generate identical time series that have properties similar to those from one and twonephron versions of the model (21, 29, 32). With all nephrons identical and starting from identical initial conditions, any given variable in a single nephron will evolve at the same pace as the same variable in all other nephrons, gradients of electrical potential or of hydrostatic pressure cannot develop, and the simulation result is therefore expected.
The phase space diagram, also in Fig. 3, was constructed using three variables: intracellular Ca^{2+} concentration in a vascular smooth muscle cell, reflecting primarily dynamics of the myogenic mechanism; the NaCl concentration in tubular fluid at the level of the macula densa, with TGF dynamics; and nephron plasma flow, with dynamics of both oscillators. The data for the phase diagram were taken from a simulation over 900 s. The trajectory formed a closed loop that repeated itself, showing the two oscillations. The primary loop is due to the slower, larger oscillation of TGF and, embedded within it, are five myogenic oscillations. The behavior of the phase space trajectory is consistent with the interpretation that the two oscillations are limit cycle oscillators synchronized at a 5:1 frequency ratio.
As noted in methods, the version of the model presented here differs from earlier ones only by the addition of more nephrons and by a more detailed vascular architecture. This architecture permits us to calculate vascular hydraulic resistances and pressure drops. Figure 4 shows the blood pressures at each of four node levels. The pressures in the nodes oscillate at the same frequencies as the nephron blood flows shown in Fig. 3. The amplitudes of the pressure oscillations increase proceeding downstream toward the afferent arterioles. The pressure oscillations in the nodes are due to oscillations in the hydraulic resistances of afferent arterioles. These results suggest that blood pressure at the origin of the afferent arterioles varies in response to both systemic arterial blood pressure and the activity of the arterioles themselves.
Broken Symmetry
In these simulations we use changes in tubular properties to introduce asymmetry into the network.
Nephron length: one nephron.
The nephron model normally uses a spatial step size of 0.01 cm and a total of 195 spatial steps. Figure 5 shows examples of nephron time series with one additional spatial step inserted into the descending limb of nephron 3 and a matching spatial step into its ascending limb. All of the nephrons in the network produced irregular time series like those in Fig. 5, left, and all had phase space attractors like those in Fig. 5, right. The phase space attractors reflect the irregularities in the time series, are not those of limit cycle oscillators, and are similar to those generated by increasing the value of the coefficient coupling TGF with the myogenic mechanism in a twonephron model (21). A length increase of this or greater magnitude generates a limit cycle oscillation in one or twonephron versions of the model (32). When the twonephron model runs with two uncoupled nephrons of different lengths, the shorter nephron has a smaller faster TGF oscillation than the longer; when coupled the two nephrons become synchronized at a single frequency for TGF, the myogenic coupled at 5:1 to TGF in both, and the amplitudes of the TGF oscillation become more nearly equal than in the uncoupled state. The interactions occurring within the network affect the dynamics of its nephrons.
A single node separates nephrons 3 and 4; both nephrons are three nodes distant from nephrons 1 and 2, five nodes from nephrons 5 through 8, and seven nodes from nephrons 9 through 16. All nephrons in this last group had identical time series and phase space attractors; nephron 12 is presented as an example from this group. Nephron 3 generates the largest amplitude fluctuations and nephron 4 the next largest, nephrons 9 through 16 the smallest. Figure 6, top, superimposes the time series from nephrons 3 and 4 and shows them for the most part to be synchronized in phase. Figure 6, middle, superimposes the time series predicted for nephrons 3 and 12. An initial period of antiphase synchrony lasting ∼100 s is followed by desynchronization and then by an additional 150 s of antiphase synchronization.
Figure 6, bottom, presents the instantaneous values of absolute phase angle difference for the two pairings. The difference for the nephron 3 vs. 4 pairing is small throughout the record, consistent with continuous inphase synchrony. The nephron 3 vs. 12 pairing generates a more variable record: the value is ∼3 radians during the initial 100 s, drops toward 0, and returns to 3 over the 150 s. Phase angle differences lie in the interval [−π, π]. The use of the absolute value constrains the value to the interval [0, π]. The values of phase angle difference from this pairing are therefore equal approximately to π and are consistent with antiphase synchronization over some of the 900s duration of the simulation. Values for the mean phase angle and the synchronization metric applied to nephron pairs 3 and 4, 3 and 12, and 4 and 12, are in Table 2. Values listed in Table 2 are means for the 900s duration of the simulation. The small mean phase angle difference between nephrons 3 and 4 reflects inphase synchronization over the length of the simulation, and the synchronization metric near the limit of 1 shows minimal dispersion of the phase angle difference. Nephrons 3 vs. 12 and 4 vs. 12 show relatively large values of the mean phase angle difference but not as high as the limit, π, because the antiphase synchrony is interrupted during the 900 s of the observation period.
An additional feature of the results is a series of amplitude changes taking place over a longer time course than the length of a TGF period. Lowfrequency fluctuations have been described in experimental time series (36, 42), but no mechanism has been suggested. Our results support the hypothesis that network dynamics are responsible.
In summary, the effect of a 1% increase in the length of a single nephron, a change that in isolation merely slows the TGF oscillation and the myogenic oscillation it modulates (32), propagates throughout the network inducing complex dynamics in all nephrons. We ascribe this change in dynamics to interactions of nephron activities in the vascular network.
Nephron length: random variation in 16 nephrons.
Lengthening one nephron introduces the idea of using variation in nephron length to break network symmetry and illustrates the kind of changes one can expect from such a maneuver. A pattern of random differences in the lengths of all nephrons is likely a more realistic approximation of the length distribution in vivo. We therefore used a random number generator to vary the lengths of each of the 16 nephrons, adding or removing spatial segments in the descending limb of the loop of Henle in discrete steps over the range [−2, 2]. Matching changes were made in the ascending limbs. Using unique seeds for the random number generator, we simulated network dynamics with 12 different length patterns.
Examples of the dynamics with these randomly distributed length changes are in Fig. 7 and show some similarities to those presented in Fig. 5. All nephrons had fluctuations with irregular amplitudes and time dependent changes in frequency. Most nephrons underwent slower fluctuations than TGF. These slower changes were usually not synchronous among more than two or three nephrons. The phase space attractor patterns were like those in Fig. 5; none was typical of a limit cycle oscillator and all had complex attractors. Supplemental Videos S1 and S2 (Supplemental material for this article is available online at the Am J Physiol Renal Physiol website) show the electrical potentials predicted in each of the 15 nodes, one from the simulation with a symmetric network, and the other from a simulation with nephrons of different lengths. In the symmetric network, all electrical potentials vary in synchrony and are the same magnitude as the potentials generated in the afferent arterioles, the source of electrical currents in the network. Electrical potentials in the nodes of the asymmetric network differ from each other, particularly during depolarization. The depolarization is greater in the nodes closer to the nephrons than in those closer to the main feeder artery.
In the symmetric network, the node vascular pressures form time series that mirror the nephron blood flows, as seen in Fig. 4. Figure 8 shows those pressures in the asymmetric network simulation whose nephron plasma flows are in Fig. 7. In Figs. 7 and 8, node 1 supplies blood to nephrons 1 and 2, node 4 to nephrons 7 and 8, node 5 to nephrons 9 and 10, and node 8 to nephrons 15 and 16. In each case, the node vascular pressure has periods of irregularity, those periods do not occur at the same times in all the nodes, and the TGF and myogenic oscillations of any given pair of node pressures do not synchronize when these irregular periods appear. In addition, there are dissimilar patterns of node vascular pressure and of the plasma flow rate in the nephrons each node supplies, reflecting loss of synchronization throughout the network.
The results in all 12 simulations were quantitatively similar to those in Fig. 7 but differed from each other because of the different patterns of nephron lengths. We applied the measures of network synchrony described in methods to determine the dependence of nephron synchronization on the number of nodes of separation. The result is in Table 3. The network structure shown in Fig. 1 provides for interactions over odd numbered node pairings. There are 8 onenode interactions, 16 threenode, 32 fivenode, and 64 sevennode interactions in each simulation.
For the nephron pairings with one or three nodes of separation the mean phase angle differences, 
In summary, by the addition or subtraction of a small number of segments from randomly selected nephrons all elements of the network are forced to take on complex dynamics. For the most part the time series retain much of the periodic behavior shown in one or twonephron simulations, but additional features arise in the network context: variation in the amplitude of TGF oscillations and lowfrequency perturbations. Variation in amplitude is known to occur as periodic systems remain synchronized while undergoing a bifurcation to chaos (35, 38). Beat frequencies are also known to arise in dynamical networks (35), but our simulations do not extend long enough to justify an analysis of the periodic nature of this phenomenon.
The simulation results show localized synchronization of small clusters of nephrons and are therefore consistent with the findings of laser speckle imaging that highly synchronized pairings occur in small groups of nephrons (19).
Synchronization mechanisms: vascular network.
Electrical and hemodynamic coupling provide the two mechanisms available to the network model to achieve synchronization. We tested the two independently to evaluate their relative importance. In one simulation, the internephron conductance was set to 0 while in the second hydrostatic pressure in nodes 1–8 was clamped at 75 mmHg. The control was the simulation with random nephron lengths, using the standard internephron conductance, and vascular node pressures left free to change, as shown in Fig. 7. The same seed number for the random number generator was used for all three simulations described in this paragraph. The results are in Table 4.
With internephron conductance set to 0, 
The results of this set of simulations suggest that the two coupling mechanisms interact and modulate each other. Fluctuations in pressure in the vascular nodes tend to produce antiphase synchronization and to reduce synchronization, and internephron electrotonic signaling leads to inphase and increased synchronization. This interaction, we suggest, restricts the effective field size of synchronization to clusters comprising a subset of the 16 nephrons and is likely responsible for the emergence of lowfrequency phenomena in the network.
Synchronization mechanisms: intranephron.
Both TGF and the myogenic mechanism act on the contractile mechanism of vascular smooth cells of the afferent arteriole, converging to create an obligatory interaction. TGF gain has higher values in hypertensive than in normal rats (8, 27), and tubular pressure time series have more complex dynamics in hypertension. The coupling of TGF to the myogenic mechanism in the single nephron model is represented as a single parameter, ζ, that modulates the Ca^{2+} conductance of smooth muscle cells. Increasing ζ in both the one and twonephron models causes a bifurcation from a limit cycle oscillation to chaos (21). Reducing the value of ζ lowers the amplitude of TGF oscillations and induces quasiperiodicity with loss of TGFmyogenic synchronization. Further reduction of ζ extinguishes the TGF oscillation leaving the myogenic oscillation as the only dynamic process. Rats with normal arterial blood pressure typically show both oscillations in tubular hydrostatic pressure records, synchronized most often at a frequency ratio of 5:1; synchronization ratios of 4:1 and 6:1 are also found but are less common (44).
All nephrons in the minimally asymmetric network model develop irregular dynamics, even when using a value of ζ too low to cause such changes in simpler configurations. To discover how these changes depend on ζ when its value is altered we ran six simulations of the network model using values of ζ, which, in the one and twonephron simulations, produced no TGF oscillations, bimodal limit cycle oscillations, or chaotic dynamics (21). Nephron lengths were varied with the same seed value for the random number generator as used in the simulation of Fig. 7. All nephrons in a given simulation were assigned the same value of ζ.
Figure 9 shows time series of nephron 9 in the 6 simulations. The lowest value of ζ used for the network simulation produced low amplitude TGF fluctuations with a burst into higher amplitude TGF oscillations, myogenic oscillations, and lowfrequency variations of the two oscillations; in the one and twonephron models there was no TGF oscillation with this value of ζ. Increasing values of ζ increased the amplitude of TGF fluctuations, with the lower frequency changes persisting. The two largest values of ζ produced time series whose TGF fluctuations varied from one cycle to the next. Phase space attractors for all time series in this group of simulations were complex. The complexity of the attractors was greatest at the extremes of the range of values and least in the midrange.
Table 5 contains values of 
The fivenode and sevennode separations produced larger phase angle differences and less synchronization than the one and threenode separations, even with intermediate ζ values. The effects of changing ζ on the five and sevennode separations were smaller than for the one and threenode results.
Proximal tubule pressures in rats with at least two different forms of hypertension have more complex dynamics than rats with normal blood pressure (16, 48, 49). TGF gain is higher in the hypertensives (8, 27), a result we that have shown capable of causing the increased complexity in the one and twonephron models (21). Increasing ζ in the network context produces similar but not identical results. The addition of hemodynamic interactions with electronic signaling in the asymmetric vascular network produces complexity in nephron blood flow time series at all values of ζ, and increasing the coupling strength of TGFmyogenic interactions increases the complexity.
Lower frequency dynamics.
In addition to TGF and myogenic dynamics, all nephron time series in Fig. 7 showed slower changes, the patterns different from each other. To understand these changes we extended the simulation shown in Fig. 7 to 4,000 s. Figure 10 shows two normalized time series: blood flows in nephron 9 and for the entire network of 16 nephrons. To facilitate comparison of the amplitudes, the values of each time series were divided by the respective mean values. Both time series contain infrequent events, but show little evidence of synchronization between them. The synchronization metric, ρ̄, was 0.25, a low value indicative of little synchronization. The mean phase angle difference was 0.40. The relative amplitude of these events is greater in the nephron blood flow time series than in the whole network result. Application of a standard FFT routine showed peaks representing TGF and myogenic rhythms in both time series, similar to those reported previously (29, 31), and broadband power without peaks, at lower frequencies. These results suggest that the patterns of desynchronization found in the shorter simulations persist for more than 1 h.
DISCUSSION
Intranephron mechanisms regulating nephron blood flow generate selfsustained oscillations that interact and synchronize, creating a functional ensemble (29). Nephron ensembles interact by exchanging electrical signals over afferent arteriolar walls (32, 47) to become synchronized (12, 19, 32, 44). We developed and exercised a model of the renal vascular network to simulate the interaction of a number of nephrons. The network model includes 16 copies of a nephron blood flow regulation model used previously (21, 29, 32), each communicating through its afferent arteriole with a branching, electrically conductive vascular tree. Branches in the network distribute blood to nephrons and serve as nodes at which electrical signals converge. The periodic constriction and relaxation of afferent arterioles also affect the hydrostatic pressure in the nodes, providing a second mechanism for nephrons to interact (30, 39, 40).
We used the symmetry of the system as a framework for the formulation of specific hypotheses and the interpretation of results. Elements in a symmetric network may synchronize (10). Nephrons and arterial branches are the elements of our network model. The network structure we used is symmetric because each branch supplies a pair of vessels with identical dimensions and electrical conductivity, and all branches of the same order elsewhere in the network do the same. The entire system is symmetric when the same nephron model is used in each of 16 locations coupled to 8 nodes in the vascular network, and the same parameter set is used for all nephrons. Each of the 16 nephrons develops an autonomous bimodal oscillation, each with the same properties. Blood pressures at the vascular nodes also show a bimodal oscillation synchronized with the nephron oscillations, and electrical potentials at the vascular nodes are synchronized with electrical potentials in the afferent arterioles. In this symmetric configuration the entire system acts as one, showing that the renal vascular network operates as can be predicted from theoretical arguments for all symmetric networks and from simulations of one or two nephrons.
We introduced an asymmetry by lengthening one of the 16 nephrons by ∼1%. Increased nephron length reduces the frequency and increases the magnitude of the TGF oscillation, a change that lowers the frequency of the myogenic mechanism in that nephron because of the interaction of the two (32). The changes propagate through the network, interfering with synchronization and causing every nephron to develop complex dynamics, unlike the simpler twonephron model, which synchronized their bimodal limit cycle oscillations with length differences as great as 10% (32). In addition to the number of nephrons, the models differ because blood pressure at the origins of afferent arterioles is allowed to vary in the network model but not in the twonephron model.
The length change in a single nephron introduces a single focus of symmetry breaking into the network and desynchronization and complex nephron dynamics ensue. Nephrons separated from each other by a single node tend to be more completely synchronized with each other than with nephrons at a greater remove. Nephrons separated by seven nodes from the single longer nephron remain synchronized with each other but less so with the single longer one. The synchronization between nephrons separated by a single node was in phase and minimally dispersed. Synchronization is in phase when electrotonic interactions predominate. When a single arterial branch supplies blood to two nephrons, there is no difference in the entrance pressures to their afferent arterioles, effectively eliminating the hemodynamic source of desynchronization and leaving electrotonic signaling to act on them alone. Nephrons separated by five or seven nodes exhibit intermittent antiphase synchronization.
The asymmetry in the kidney is likely to involve more than 1 of a group of 16 nephrons. We therefore used a random number generator to provide different patterns of length variation among the 16. A similar response emerged in each of 12 simulations, each with a unique pattern of length variation. Each nephron, regardless of its own length, developed complex blood flow dynamics. Synchronization of TGF with myogenic dynamics in individual nephrons was not preserved and fluctuations at frequencies lower than that of TGF emerged, a result consistent with experimental observations (36, 41). Synchronization among nephrons varied inversely with the number of nodes of separation, and, as might be expected, synchronization of electrical events among nodes of the vascular network was weaker in the asymmetric configuration than in the symmetric one, as shown in Supplemental Videos S1 and S2. In any given simulation, some nephrons appeared to have experienced only a minor disturbance to their fundamental periodic behavior while others showed more substantial changes. All nephrons in every simulation had nonperiodic changes in the amplitude of the TGF oscillation accompanied in some cases by changes in frequency. These very lowfrequency disturbances occurred in all nephrons but usually not all at the same time.
Dynamical systems of order 3 or higher may undergo bifurcations from regular oscillations to chaos if the value of some parameter of the system changes. Despite the bifurcation, the system can retain significant periodic behavior and synchronization of chaotic systems is a wellknown phenomenon (1, 35, 37, 38). Whether the dynamics of the asymmetric configuration represents chaos is a question we have chosen not to address in this work because it would considerably expand the scope, and we have therefore described the dynamics as irregular or complex, without implying any difference in the application of one or the other adjective.
Clamping entrance hydrostatic pressures for all afferent arterioles allows electrotonic signaling to predominate throughout the network, synchronization is then in phase; eliminating electrotonic conduction and allowing vascular pressure to vary spontaneously favors antiphase synchronization among all nephron pairs. With both mechanisms intact, electrotonic signaling among nephrons of different lengths generates coupled nephron pairs with frequencies and amplitudes that differ among pairs, a change that disrupts synchronization of blood pressure among nodes.
Nephrons attached to a single node are supplied by blood at the pressure in that node. Electrotonic signaling between them is then free to establish inphase synchronization, as it does in the twonephron model (32). Antiphase synchronization occurs when nephron pairs are supplied from different nodes; synchronized vasoconstriction in one or two pairs of nephrons can raise the pressure in upstream nodes, a change to which other nephron pairs react. The reaction in these other nephrons occurs after a delay due to TGF, leading to antiphase synchronization. Local differences in blood pressure will affect glomerular filtration rate, causing variation in the rate of fluid delivery to the macula densa regions of the nephrons, leading to variation of signaling to afferent arteriolar smooth muscle cells, and causing variable changes in the amplitude and frequency of TGF oscillations. The patterns of nephron interactions with both vascular interaction mechanisms intact therefore reflect a dynamic balance between them, acting in a complex pattern throughout the network, the pattern changing with time.
The dependence of synchronization on the number of nodes of separation could have functional consequences. The effects of nephrons with different properties on dynamics within the network would be reduced because the loss of synchronization would provide functional isolation. Medullary nephrons, for example, are longer than the cortical nephrons we studied here; single vascular networks appear to supply both groups with blood. Because of their length, medullary nephrons might be expected to alter the dynamics of cortical nephrons in a shared vascular cluster, but the functional isolation provided by desynchronization could serve to reduce the impact of the interaction between the medullary and cortical nephrons.
The asymmetric network model predicts slow changes in blood flow of individual nephrons and of the whole network, as shown in Fig. 10. The individual nephrons in a network do not undergo these changes simultaneously, and neither in the network nor individual nephrons are these changes periodic. Siu et al. (42) reported 10mHz oscillations in rat whole kidney blood flow, and Pavlov et al. (36) detected oscillations in proximal tubule hydrostatic pressure at a similar frequency, also in rats. In both reports, the spontaneous lowfrequency activity was small in normotensive animals and enhanced by forcing renal artery pressure over a range of frequencies (42). The network model, by design, simulates a subset of the range of diversity present in nephrons and blood vessels in the kidney. It would be premature to conclude that the model will not simulate these experimental results.
The experimental observations that nephrons interact through their vascular connections (12, 20) and that proximal tubule pressure oscillations found in normal rats are replaced by more complex dynamics in hypertensives (16, 48) motivate the work reported in this paper. A number of modeling studies have addressed various aspects of the problem. Single nephron models without an active myogenic mechanism have produced selfsustained TGF oscillations (15, 17), and in some cases can bifurcate to chaos (2, 43). Layton et al. (22–25) developed a model of thick ascending limb dynamics to study the dynamics of TGF regulation; the model lacks representation of the myogenic mechanism, and blood pressure is not used as an input. They have used their model to show that coupling two or three nephrons can produce oscillations or irregular fluctuations, depending on parameter choices. The coupling was modeled on electrotonic signaling. Bayram and Pitman (3) used the model of Layton et al. (22) in a multinephron configuration and found enhanced stability but did not use a parameter set that produced complex dynamics. Postnov et al. (39) assembled a multinephron model to examine blood flow dynamics, viewing the vascular structure as a resource distribution network; they used the single nephron model of Barfred et al. (2) but did not include electrotonic messaging, Marsh et al. (30) studied a multinephron network, also using the Barfred model, and included both internephron signaling and vascular pressure dynamics. The model lacked a myogenic representation and used a simple vascular network structure; it developed some complex dynamics. None of the models cited in this paragraph had a complete set of the interactions known to be active in this system, and none developed the complexity we report here.
On close examination the experimental data from normotensive animals usually contain timedependent changes at lower frequencies than TGF, a result consistent with the dynamics our network model predicts (26, 48–50) but not previously simulated. We had interpreted our experimental findings to indicate the presence of limit cycle oscillations in tubular pressure, flow, and early distal tubule Cl^{−} concentration (18). Limit cycle oscillations are time invariant, a condition not satisfied by any of the results in this work when the network is not symmetric. The results of Fig. 7, for example, show different levels of complexity in different nephrons, and synchronization among small clusters. Those showing less complexity provide time series that are consistent with experimental findings.
The output of a complex system has greater bandwidth and information content than a limit cycle oscillator (37). In its symmetric configuration, the network was fully synchronized, all nephrons functioning as limit cycle oscillators forcing the elements of the vasculature at the same frequencies. Because the nephron oscillations were strictly periodic and the phase angle differences zero, the message content is highly redundant, and a full characterization of the network at any time requires a single specification of the frequency and amplitudes of the TGF and myogenic oscillations. In the asymmetric configuration, all nephrons have different amplitudes, frequencies, and phase angles, and the electrical potentials and pressures in the vascular nodes are dissimilar, all time dependent, and therefore with reduced redundancy of information in the message. The collection of nephrons and blood vessels, which is understood as a resource distribution network, is therefore also an information network. The specification of the information content of a signal does not provide semantic meaning, and none of the simulations reported here were designed to provide insight on this matter. It seems reasonable to speculate that the signals in this network are involved in nephron blood flow regulation, but whether this information influences other nephron functions and at what time scales is not known.
Nordsletten et al (34), making use of a scaling system to rank vessel diameter and length, found that each artery in the renal vascular network could connect with vessels of two or three different orders, and their figures show variable branching patterns, confirming the earlier report of Casellas et al. (5). The vascular structures are therefore likely to be asymmetric. The effects of variation in vascular anatomy remain to be studied.
Significance
The focus in this work is on the interactions that occur when nephrons are aggregated in a network that supplies them with blood and conducts electrical signals among them. TGF and the myogenic mechanism interact and form a regulatory ensemble in each nephron, and nephron ensembles signal each other electrically. Both sets of interactions affect afferent arteriolar hydraulic resistance, and the changes in nephron blood flow modulate vascular pressure in the network, creating an additional mechanism for the interaction of nephron ensembles. Complex dynamics result from these interactions so long as the network is asymmetric, as normal variation in nephron structure assures it must be. The results suggest that asymmetries due to the presence of medullary nephrons and to variations in vascular branching patterns will have additional effects on synchronization within renal vascular networks.
Analysis of experimental data has suggested that nephron pressures and flows oscillate in animals with normal blood pressure and exhibit chaotic dynamics in those with chronic hypertension. Simpler nephron models or those with many fewer nephrons have supported this interpretation. The results presented in this study suggest an alternative: nephron dynamics are normally complex, and parameter changes that occur in hypertension increase their complexity. Messages passed through the symmetric network are highly redundant and therefore have relatively low information content. The complex dynamics we found implies that electrical message passing through the network, sent and received by afferent arterioles supplying individual nephrons, is much less redundant and therefore has greater information content, and that the information content varies as nephrons adjust their dynamics.
GRANTS
The work of D. J. Marsh was supported by grants from the Lundbeck Foundation of Denmark, the University Kidney Research Organization of Los Angeles, California, and the Office of the Dean of Medicine and Biological Sciences at Brown University. The work of N.H. HolsteinRathlou was supported by the Danish Medical Research Council and the Novo Nordisk Foundation.
DISCLOSURES
No conflicts of interest, financial or otherwise, are declared by the author(s).
 Copyright © 2013 the American Physiological Society