## Abstract

Tubular pressure and nephron blood flow time series display two interacting oscillations in rats with normal blood pressure. Tubuloglomerular feedback (TGF) senses NaCl concentration in tubular fluid at the macula densa, adjusts vascular resistance of the nephron's afferent arteriole, and generates the slower, larger-amplitude oscillations (0.02–0.04 Hz). The faster smaller oscillations (0.1–0.2 Hz) result from spontaneous contractions of vascular smooth muscle triggered by cyclic variations in membrane electrical potential. The two mechanisms interact in each nephron and combine to act as a high-pass filter, adjusting diameter of the afferent arteriole to limit changes of glomerular pressure caused by fluctuations of blood pressure. The oscillations become irregular in animals with chronic high blood pressure. TGF feedback gain is increased in hypertensive rats, leading to a stronger interaction between the two mechanisms. With a mathematical model that simulates tubular and arteriolar dynamics, we tested whether an increase in the interaction between TGF and the myogenic mechanism can cause the transition from periodic to irregular dynamics. A one-dimensional bifurcation analysis, using the coefficient that couples TGF and the myogenic mechanism as a bifurcation parameter, shows some regions with chaotic dynamics. With two nephrons coupled electrotonically, the chaotic regions become larger. The results support the hypothesis that increased oscillator interactions contribute to the transition to irregular fluctuations, especially when neighboring nephrons are coupled, which is the case in vivo.

- renal autoregulation
- nonlinear dynamics
- tubuloglomerular feedback
- myogenic mechanism
- chaos

two mechanisms, tubuloglomerular feedback (TGF) and the myogenic mechanism, operate in each nephron of mammalian kidneys to regulate blood flow when arterial blood pressure changes (14–16, 22, 44). Each mechanism is potentially unstable and operates in a nonlinear regime: TGF oscillates at 0.02–0.04 Hz and the myogenic mechanism at 0.1–0.2 Hz. The TGF oscillation is the more pronounced. These mechanisms interact because they operate simultaneously on the contractile machinery of vascular smooth muscle cells. The interaction leads to modulation of the amplitude and frequency of the myogenic oscillation by TGF, a finding predicted by a mathematical model (26) and confirmed in experimental measurements (28).

The dynamics of the TGF-myogenic ensemble undergo a transition in rats with either a genetic or a renovascular form of hypertension to a state with characteristics of deterministic chaos (51). Lyapunov exponents are positive; the phase space attractor constructed from the data is low dimensional and has a noninteger correlation dimension (51); surrogate data analysis confirms the nonlinearity of the system (52).

The question at issue is the cause of the transition. Dilley and Arendshorst (8) and Leyssac and Holstein-Rathlou (23) reported that the feedback gain of TGF, a component of the coupling of TGF to the inherent dynamics of the afferent arteriole, is greater in hypertensive rats than in rats with normal blood pressure. In our model, an increase in TGF gain or in the parameter expressing the coupling, or both, will increase the strength of the interaction. In this study, we test the hypothesis that the strength of the TGF-myogenic interaction functions as a bifurcation parameter and, when increased, can cause a bifurcation to chaos. For the test, we use a mathematical model of a single nephron and its afferent arterioles (26), and of two such nephrons coupled by electrotonic vascular signal conduction (29). With an increase in the strength of the interaction, the single nephron model succeeds in predicting transitions to chaotic dynamics, but the regions with these dynamics in a one-dimensional parameter scan appear to be small and therefore inadequate to account for the experimental observations. The simulation of two nephrons coupled electrically, which is closer to the natural state than is the single nephron model, generates larger domains of chaos than the model of one nephron alone. The results obtained with the coupled nephron model therefore support the hypothesis that increased coupling between TGF and the myogenic mechanism can cause a change in the dynamics of nephron blood flow regulation similar to that found in rats with chronic high blood pressure.

## METHODS

### Model

The essential elements and mechanisms of the model are summarized in the causal loop diagram in Fig. 1. Signs at the arrowhead indicate the direction of change of a given variable in response to a change in the variable at the arrow's tail. The elements of the feedback loop that generate TGF and its oscillation comprise glomerular filtration (GFR), the axial mass transport of NaCl and water, and the epithelial transport of NaCl and water in various tubular segments. These elements are shown in Fig. 1, *right*. Epithelial transport of NaCl is dependent on axial flow rate in several of these segments, so that the concentration of NaCl in tubular fluid at the macula densa is strongly dependent on the axial flow rate of tubular fluid. The apical membranes of macula densa cells have transporters sensitive to ionic concentrations, chiefly Cl^{−}. The signal from the macula densa produced in response to variation of the NaCl concentration in tubular fluid propagates through a number of cells and finally affects the Ca^{2+} conductance of vascular smooth muscle cells. The change in Ca^{2+} conductance alters intracellular Ca^{2+} and modifies the length of the contractile mechanism, altering the arteriolar vascular resistance. The change in vascular resistance affects GFR, closing the loop and providing negative feedback.

The myogenic mechanism is represented by the loops in Fig. 1, *left*. The period of the oscillation is determined by instability in membrane potential caused by an interaction between Ca^{2+} and K^{+} currents. The magnitude of the oscillation in vascular diameter is determined by the level of intracellular Ca^{2+}. TGF affects the membrane Ca^{2+} conductance that is part of the myogenic mechanism in the model, allowing TGF to modulate both the amplitude and frequency of the myogenic mechanism. The magnitude of the TGF stimulus decreases with distance from the glomerulus, proceeding along the afferent arteriole toward the nearest branch point of the artery supplying blood to the nephron. A two-point approximation to this continuous decline gives the arteriolar model two distinct segments, each modeled separately and the one closer to the glomerulus receiving the greater TGF input. The arteriolar segments are coupled electrically. Details of the mathematical model of a single nephron and of a nephron pair were presented in previous work (26, 28). We used the same models in this study.

In the first of these papers, tubular pressure, flow, and NaCl concentration were modeled as functions of time and distance in a single nephron with a compliant epithelial wall that reabsorbs water and NaCl. The nephron consisted of a proximal tubule, a descending limb of Henle's loop, a thick ascending limb of Henle's loop, and an early distal tubule. Appropriate parameters were used for each tubular segment. The boundary conditions were GFR as the initial flow rate, a nonlinear expression for pressure at the end of the early distal tubule, reflecting the compliance of later segments of the tubule (24), and NaCl concentration in the glomerular ultrafiltrate. NaCl is the only tubular solute represented in the model, and its concentration is assumed to remain unchanged in the proximal tubule.

Each of two segments of the afferent arteriole was modeled with six ordinary differential equations whose dependent variables were membrane K^{+} conductance, intracellular Ca^{2+} concentration, membrane electrical potential difference, myosin light chain phosphorylation, length of the contractile mechanism, and length of two different sets of elastic elements, one in series with the contractile mechanism forming an ensemble, which is in parallel with the second elastic component. The length of the parallel elastic element was equal to the sum of the lengths of the series element and of the contractile mechanism; the circumference of the vascular lumen was calculated from this length. This cellular model of arteriolar dynamics is based on the work of Gonzalez-Fernandez and Ermentrout (9) which was designed to simulate vasomotion in cerebral arterioles.

The TGF signal was calculated from the tubular NaCl concentration at the end of the thick ascending limb with a logistic equation in a form used to fit experimental data; TGF output was used to modify membrane Ca^{2+} conductance to cause changes in smooth muscle contraction and to simulate the interaction between TGF and the myogenic mechanism. Both segments of the afferent arteriole generated spontaneous self-sustained oscillations because of the interaction between voltage-gated Ca^{2+} channels and Ca^{2+}- and voltage-dependent K^{+} channels. TGF stimulation of the arteriole produces maximal contraction at the point closest to the glomerulus; the vasoconstriction declines with distance from the glomerulus (29, 47); the fractional decrease data can be fitted with a single exponential. Two arteriolar segments were used as a two-point approximation to this decline; the TGF effect was greater in the segment close to the glomerulus and smaller in the farther segment. The two segments were coupled electrically with an ohmic conductance.

The interaction between TGF and the myogenic mechanism was expressed as: *g*_{Ca,j,c} = (1 + ζ_{j}·θ_{j})·*g*_{Ca}, where *g*_{Ca} is the native Ca^{2+} conductance, *g*_{Ca,j,c} is the Ca^{2+} conductance of the *j*th arteriolar segment when coupled to TGF, θ_{j} is TGF input to the *j*th arteriolar segment, and ζ is the parameter coupling the two mechanisms. We varied ζ to change the strength of the connection between TGF and the myogenic mechanism to test the hypothesis that ζ can act as a bifurcation parameter and explain the different dynamics observed in normal and hypertensive rats.

The two-nephron model consists of two versions of the single nephron model, each solved separately, but with the arteriolar segments farther from the glomerulus coupled electrically (29). These arteriolar segments were coupled to each other through an electrical node in the wall of the artery supplying the two nephrons with blood. In this study, cortical nephrons were simulated, and both nephrons of the pair were assigned the same length. Additional parameters needed for the two-nephron model were the conductance coupling the two nephrons through an electrical node, and a conductance of the node to ground. The model was otherwise used unmodified from the single nephron form.

#### Numerical methods.

The partial differential equations describing pressure, flow, and NaCl concentration in each renal tubule were solved using centered difference approximations and are second-order correct (40). The spatial step was 0.125 mm and the time step 1 ms. Reducing either the spatial step or the time step by a factor of 10 had no effect on the values of ζ required to produce bifurcations in the model results. The auxiliary equations, including the boundary conditions, were solved with the Newton Raphson method. All calculations were performed in double precision.

The solutions of the ordinary differential equations were obtained with Gear's variable step size method, using backward differentiation and numerically generated Jacobians. Because the solution of the partial differential equations required iterations at each time step, the time steps for the solution of the ordinary differential equations were synchronized with those used for the partial differential equations.

For each time step, the glomerular model was solved by using the afferent resistance and the tubular inlet pressure from the preceding time step as initial estimates. With the use of the calculated value for GFR for the inflow rate to the tubular system, the pressure, flow, and NaCl equations were solved iteratively at each time step of the whole system. Convergence was assumed when the fractional change of the Euclidian norm of the combined pressure and flow vectors was <10^{−3} and <10^{−5} for the NaCl concentration vector. The new calculated value for the NaCl concentration at the macula densa was applied to estimate the afferent arteriolar resistance. The procedure was repeated, using the new value for the afferent resistance and the tubular inlet pressure, until the system of equations converged. Convergence was assumed when the change in the Euclidean norm of the afferent resistance was <10^{−6} in successive iterations. Three or four iterations were typically needed to achieve convergence.

#### Analytic approach.

The model generates time series corresponding to solutions for tubular pressure, volume flow rate, and NaCl concentration at each spatial and time step and for each of the 12 ordinary differential equations used to simulate the afferent arteriole. We will present only nephron plasma flow rate, NaCl concentration in tubular fluid at the macula densa, and intracellular Ca^{2+} concentration in afferent arteriolar smooth muscle cells. These variables were chosen because, when combined, they form a phase space attractor with one variable, NaCl concentration at the macula densa, showing primarily TGF dynamics; a second, intracellular Ca^{2+} concentration in smooth muscle cells, representing dynamics mainly of the myogenic mechanism; and a third, nephron plasma flow rate, with the dynamics of both mechanisms.

The results of our simulations show limit cycle oscillations, quasiperiodicity, and chaos for different values of ζ, the parameter that couples TGF with the myogenic mechanism in the nephron model. In the phase space, each of these kinds of dynamics represents a specific type of attractor, i.e., invariant (steady) state for the system. A bifurcation diagram provides an overview of the different kinds of dynamics generated in the system. In such a diagram, the coordinates of the intersection between the trajectory of the system and a plane (a Poincaré section) in the space of state variables is mapped on the *y*-axis as a function of a parameter, in this case ζ.

The bifurcation diagrams are calculated by the following procedure. *1*) A Poincaré section is defined, which in our case is chosen to be *I*_{K} − 0.1 = 0, d*I*_{K}/d*t* > 0, where *I*_{K} is the K^{+} current in the vascular smooth muscle cell*. 2*) Starting with a value of ζ = 0.18, the state variables are initialized. *3*) The system is simulated for 3,000 s to allow the solution of the system to converge to a stationary state. *4*) During an extension of the simulation by 3,000 s, the intersections with the Poincaré section are calculated by using a linear interpolation between the states just before and after the intersection. *5*) The parameter ζ is increased by 1.8 × 10^{−4}, and the state variables are maintained at their current state (adiabatic initial conditions). *6*) The procedure is repeated from *step 3* until ζ reaches its final value of 0.4.

The model generated four main types of dynamics as ζ changed: a state in which only the myogenic oscillation is present in the blood flow time series and states with periodic, quasiperiodic, or chaotic dynamics. Chaotic systems have the property that they are sensitive to initial conditions of the state variables. This means that the evolution in time of two sets of initial conditions, differing initially by an arbitrarily small amount, on average diverge exponentially. For nonlinear periodic systems, in contrast, the orbits will approach and relax onto a single limit cycle. A system in a quasiperiodic state is nonperiodic but does not show sensitivity to initial conditions. In fact, the two orbits remain different by a value that at most grows proportionally with time, and not exponentially.

Lyapunov exponents are often used to classify dynamics. In the present case, however, the system is too complicated to calculate a complete spectrum of exponents by the first principle approach. Instead, we used the method that calculates the largest exponent from a time series (48). This method is based on the idea of attractor reconstruction by embedding the time series together with past versions of the same time series into a multidimensional space. Relative divergence of nearby points on the trajectory is averaged along the trajectory. The embedding dimension may be chosen arbitrarily large, since the estimated exponent will converge to a fixed value. However, the drawback is that the larger the embedding dimension the slower the convergence. For the present case, the time series is represented by nephron plasma flow, and embedding dimension was set to 12. The time series was 1,500 s long. A positive Lyapunov exponent is an indication of the sensitivity to initial conditions, which is a fundamental characteristic of chaotic dynamics, whereas periodic and quasiperiodic dynamics show a negative and zero-valued exponent, respectively.

## RESULTS

We wish to test whether systematic variation in ζ causes changes in the dynamics of the nephron. No changes were made in other parameters, or in the arterial pressure, which was 100 mmHg in all simulations reported here.

For different values of ζ, the model developed one of four different kinds of time-dependent behavior. At very low values, there was no TGF oscillation, and only the myogenic oscillation remained active. The three other patterns are shown in Fig. 2. At relatively low values of ζ, there was quasiperiodic behavior, which is characterized by the presence of two or more noncommensurate frequencies. Quasiperiodic motion has marginal stability, as indicated by the vanishing of the first two Lyapunov exponents. Figure 2, *top*, with ζ = 0.24, shows quasiperiodicity. The two oscillations synchronize at higher values of ζ, shown in Fig. 2, *middle*, for which ζ = 0.275. The larger oscillations are because of TGF. For each TGF cycle, there are five smaller peaks, each resulting from the self-sustained oscillation of the myogenic mechanism. The five peaks are not equally spaced, nor are they of equal magnitude in each TGF cycle. These variations among the myogenic peaks are the result of frequency and amplitude modulation (26, 28). The fluctuation pattern becomes irregular at higher values of ζ, shown in Fig. 2, *bottom*, which was generated using ζ = 0.35. The TGF fluctuations persist, but the magnitude varies from cycle to cycle and in no apparent order. There were three to five myogenic peaks for each TGF cycle, and the number of myogenic peaks per TGF cycle varied in no apparent order in successive cycles. The time series in Fig. 2 are taken from the output of one of the nephrons in the two-nephron model. The two nephrons in that model had identical properties and parameters, and the time series of the two were identical. The time series from the one-nephron model showed minimal differences from those shown in Fig. 2 and are therefore not shown here.

Each type of motion shown in Fig. 2 forms a characteristic attractor in phase space. A later analysis of the results makes use of a Poincaré section, which is the intersection of a phase space trajectory with a plane in phase space, and the phase space diagrams of the model results are therefore shown in Fig. 3. In each case, the phase space diagram is formed from the simulation results shown in the corresponding panels of Fig. 2, using the predicted time series of single nephron blood flow rate, NaCl concentration in tubular fluid at the macula densa, and intracellular Ca^{2+} in an arteriolar smooth muscle cell. Figure 3, *top*, shows the results with ζ = 0.24. The pattern is that of a two-dimensional torus, with the large motion resulting from TGF and the smaller motions attributable to the myogenic oscillation dispersed around the TGF trajectory. If simulated to infinity, the orbit will densely fill the torus, without repeating itself. The quasiperiodic behavior may appear periodic in Fig. 2, *top*, but the phase space attractor in Fig. 3 shows that the orbit runs on a torus. A closer look at the time series in Fig. 2, *top*, also confirms the aperiodic nature of the dynamics.

Figure 3, *middle*, shows the trajectory formed by a system with a stable oscillation. This tracing represents 600 s of simulated time, the same as in Fig. 3, *top* and *bottom*. The combined motion of the two oscillations retraces itself more than 14 times during this simulation without change from one repetition to the next, which is typical of phase-locked limit-cycle oscillations. Such oscillations are specific modes on the torus surface in which there is a rational ratio of the two periodicities. In our case, five myogenic oscillations are completed during each TGF cycle. The ability to synchronize different oscillatory modes is a characteristic feature of nonlinear systems. The myogenic trajectories are not equally spaced around the TGF loop because of frequency modulation by TGF, and they are not of equal amplitude because of amplitude modulation by TGF (26, 45).

Figure 3, *bottom*, shows the phase space diagram formed by the system at a higher value of ζ. The attractor no longer lies on the surface of a torus, and it does not retrace itself during the 600-s duration of the simulation. A chaotic attractor cannot run on a two-dimensional torus because two-dimensional dynamics are insufficient for chaos to occur and its attractor is therefore more complicated. The myogenic oscillations appear to differ in magnitude at different locations in phase space, suggesting that amplitude modulation persists in this chaos-like state. The time series in Fig. 2, *bottom*, also shows amplitude modulation of the myogenic oscillation by TGF.

To provide a quantitative description of the bifurcations shown by the model, we calculated one-dimensional bifurcation diagrams, which show the bifurcations through which chaos-like dynamics are developed as ζ is increased. The bifurcations taking place on the route from limit cycle to chaos-like behavior involve a torus bifurcation followed by saddle-node and homoclinic bifurcations. A torus birth bifurcation is a transition in which a periodic orbit loses its stability as two complex eigenvalues (called Floquet multipliers) cross out of the unit circle in the complex plane. In the present case, this happens when the newly born cycle (TGF in this case) modulates the existing cycle (from the myogenic mechanism) to form a torus in phase space. A saddle-node bifurcation is characterized by the collision of a stable and an unstable limit cycle during which both limit cycles cease to exist and the system seeks other attractors in the phase space. Bifurcations of this type are involved in the transitions between quasiperiodic and phase-locked periodic dynamics on the surface of the torus.

Homoclinic bifurcations, on the other hand, are involved in the destruction of the torus and the transition to chaos as the interaction between the two modes becomes too strong (1). A homoclinic bifurcation occurs, for instance, when a periodic orbit expands to touch the unstable manifold of an unstable equilibrium point. After the bifurcation, the periodic orbit no longer exists.

At low ζ values, only the myogenic mechanism oscillates. In the bifurcation diagram for the single nephron model (Fig. 4*A*), this mode is represented by a single curve. The TGF mechanism does not oscillate at this value of ζ because the signal from the macula densa is too weak. At ζ_{1} = 0.221, destabilization of the myogenic limit cycle takes place through a torus birth bifurcation. Oscillations are now present in both mechanisms.

In the region after ζ_{1}, there is an infinitely dense sequence of intervals, called windows or resonance tongues, with various mode lockings between the TGF and myogenic oscillations. Windows are terminated on both sides by saddle-node bifurcations and the transition to quasiperiodic dynamics. At ζ_{2} = 0.267, a wide window with 5:1 mode locking is born through a saddle-node bifurcation. The 5:1 mode locking is represented as five simultaneous curves, since the orbit makes five windings (myogenic oscillations) for each TGF cycle. This mode is present in the normotensive state. It is terminated at ζ_{3} = 0.309 through another saddle-node bifurcation, where the regular dynamics transform into chaos-like dynamics. For increasing ζ from this point, alternating regions show chaos-like and 5:1 and 4:1 dynamics.

The bifurcations in the two-nephron model (Fig. 4*B*) take place in much the same way as the single nephron model, but at different values of ζ. The torus bifurcations take place at ζ_{1} = 0.221, which as before leads to narrow windows with mode lockings and quasiperiodicity. The saddle-node bifurcation at ζ_{2} leading to 5:1 mode locking takes place at ζ = 0.249, a lower value than for the single nephron model, and is terminated at higher values, yielding a broader region with 5:1 dynamics. Most of the periodic windows that appear at larger ζ in the single nephron model are absent, and the chaos-like motion is present over the interval ζ = [0.31, 0.355]. At ζ_{4} = 0.355, a homoclinic bifurcation cause the emergence of 4:1 mode locking, which at higher ζ bifurcates to 5:1 mode locking through homoclinic bifurcations.

Figure 5 shows the Lyapunov exponent as a function of ζ. With the single nephron results (Fig. 5*A*), the exponent has negative values at the mode-locked regions and negative or almost zero exponents at ζ values where quasiperiodic dynamics occur. The narrow regions with positive exponents are chaotic. The values of the Lyapunov exponent with the two-nephron model (Fig. 5*B*) present a similar pattern but reflect the different values of ζ at which bifurcations occur.

## DISCUSSION

Arterial blood pressure, the major input to the systems regulating renal blood flow, fluctuates irregularly (13, 25). The pattern of fluctuations in the ultradian frequency band is 1/*f*, signifying that the logarithm of the spectral power varies inversely with the logarithm of the frequency; small fluctuations are more frequent and large ones less so. The fluctuations are large enough to alter distal tubule flow rate, and therefore to disturb the epithelial transport mechanisms responsible for regulation of the composition and volume of body fluids. Renal autoregulation serves as a filter to reduce the impact of pressure fluctuations at frequencies <100 mHz (21, 41). Because the blood pressure is not steady, analysis of the time-dependent behavior of autoregulation is essential; steady-state constructs cannot, in principle, produce an adequate understanding of how autoregulation protects flow-dependent transport processes so that they can respond appropriately to hormonal control important for regulation of the body fluids. Moreover, the kidney develops in an environment with time varying arterial pressure, hormones, and other factors that affect it. The organization of renal structure and function must therefore reflect a response to these ongoing challenges, responses that cannot be understood with knowledge only of steady-state function. Oscillations, modulation effects, nephron synchronization, and bifurcations to chaos, all of which have been observed experimentally, arise out of time-dependent processes.

The delayed negative feedback control associated with TGF is responsible for periodic oscillations of tubular flow rate, GFR, and of tubular fluid NaCl concentration at the macula densa in time series of tubular pressure from rats with normal blood pressure (16). The arterioles supplying blood to these nephrons undergo periodic vasomotion and cause another, higher-frequency oscillation in tubular pressure and nephron blood flow (50). The periodic oscillations become irregular in rats with chronic high blood pressure, and the time series have properties of chaos (15, 51).

In this study, we tested whether changes in the strength of coupling between the TGF and myogenic mechanisms change the dynamics of blood flow control. We used a spatially extended model of tubular pressure, tubular flow rate, and tubular NaCl concentration in a compliant tubule to address the question. Boundary conditions for the system equations were GFR as the initial flow rate, the tubular pressure at the furthest extent of the distal tubule we simulated, and the concentration of NaCl in the glomerular filtrate. Because NaCl is the only solute whose epithelial transport is simulated in the model, its concentration is assumed invariant through the proximal tubule. The differential equation describing the change in NaCl concentration is therefore solved only from the transition between these tubular segments to the end of the model's tubular system.

GFR is heavily dependent on hydrostatic pressure in glomerular capillaries, and this pressure is dependent on the vascular resistance of the afferent arteriole that supplies blood to the glomerulus. The vascular resistance is affected by several factors: the length dependence of the contractile mechanism, the operation of a self-sustained oscillation in arteriolar radius caused by an instability in membrane electrical potential difference, input from TGF, and interactions among these factors. These components of the system act to maintain GFR constant over the range of arterial blood pressures occurring during a day (13, 25), and they have nonlinear properties. Finally, both the simulated tubule segments and those beyond the simulation show a significant mechanical compliance in response to changes in the tubular pressure (24). Hydraulic resistance to fluid outflow from the tubule therefore declines as flow increases. These boundary conditions provide nonlinearities to the system that contribute to the bifurcations we have simulated.

The model succeeds in simulating a number of measures of renal blood flow performance (26). Time-averaged results of GFR, of filtration fraction, of tubular flow rate, of tubular hydrostatic pressure, and of NaCl concentration in tubular fluid at the macula densa correspond to those found in experimental measurements. The model also autoregulates GFR and renal plasma flow in response to changes in mean arterial pressure at flow rates characteristic of experimental responses (26). The dynamic predictions of the model include the presence of two oscillations at measured amplitudes and frequencies, sustained oscillations over a large range of arterial blood pressures and consistent with experimental studies of renal blood flow regulation, and modulation of the amplitude and frequency of the myogenic mechanism by TGF (26). Nonlinear interactions are found in measurements of whole kidney blood flow (6) and tubular pressure (28) and predicted in simulations with this model (26).

In this study, both the single nephron and two-nephron simulations undergo torus, saddle-node, and homoclinic bifurcations as the value of the coupling coefficient ζ is increased. This result is consistent with experimental findings (14, 15, 51). For the single nephron model, chaotic dynamics appear only for a small set of discontinuous values of the coupling parameter, whereas the two-nephron model behaves chaotically over a larger and more connected region. Nephrons are normally coupled to each other because of electrotonic signal propagation over the wall of the afferent arteriole (12, 29, 45), and the behavior of the two-nephron model is therefore more relevant to the kidney in vivo. The model also predicts that the bifurcation from quasiperiodicity to the 1:5 resonant oscillation in TGF occurs at lower values of ζ in the two-nephron model. Nephron coupling therefore increases the range of values over which both limit cycle oscillations and chaotic dynamics operate. Previously, we showed that nephron coupling improves the stability of the model with respect to the range of a parameter used to maintain the myogenic-to-TGF frequency ratio at 5:1 (29). Beach et al. (4) concluded that the TGF oscillation is more stable if nephrons interact than it is in single nephrons. The renal vascular structure is an example of a resource distribution network that can be a source of complex dynamics from both signaling and cardiovascular sources (27, 38, 39).

The model results suggest that an increased interaction between the TGF and myogenic oscillators can cause the bifurcation to the chaotic state. When the TGF oscillation is a limit cycle oscillator, there is synchronization between the two mechanisms at a 5:1 ratio, both in experimental results (45) and in the model. Increasing the value of ζ to produce chaos interferes with the synchronization between the mechanisms. We have shown previously that the interaction leads to synchronization and modulation effects (28). Clamping TGF input at a fixed value eliminated the time-varying modulation and allowed the myogenic mechanism to oscillate at a single frequency (26).

Using a simple nephron model lacking a representation of the myogenic mechanism (3), we have modeled transitions to chaos (18, 44). These earlier results showed period-doubling cascades as the bifurcation pathway from limit cycle oscillation to chaos. In this study, the pathway was through torus bifurcations. Period-doubling cascades have not been observed experimentally; tubules in the kidney are either oscillating or in a chaotic mode (15, 45, 49). The transition between the two states in vivo is therefore likely to be abrupt, which is characteristic of torus bifurcations. Abrupt transitions between dynamic states were found consistently in this study, as seen in Fig. 4. This distinction reinforces the conclusion that increased interaction between TGF and the myogenic mechanism can cause the bifurcation to chaos that occurs in animals with chronic hypertension.

TGF gain is increased in spontaneously hypertensive rats (8, 14). In our model, the gain is the product of the TGF output and the coupling parameter ζ. The coupling is assumed to act through the voltage-dependent Ca^{2+} conductance in the membrane of arteriolar smooth muscle cells. Strong experimental evidence suggests that the myogenic mechanism acts though membrane depolarization followed by Ca^{2+} influx through voltage-activated channels (7). Substantial evidence has accrued that TGF acts through activation of adenosine-sensitive A_{1} receptors (43, 46). The adenosine appears to be derived from interstitial dephosphorylation of ATP released from the macula densa cells (42). The exact localization of the adenosine A_{1} receptors mediating the TGF response is not known, but recent data suggest that adenosine receptors on vascular smooth muscle cells are critical for the TGF response (35). Activation of adenosine A_{1} receptors on the vascular smooth muscle cells of afferent arterioles activates phospholipase C, leading to increased inositol trisphosphate production with subsequent release of intracellular Ca^{2+} (10). The major effect of this initial and transient increase in intracellular Ca^{2+} seems to be an activation of a Ca^{2+}-dependent Cl^{−} channel that results in membrane depolarization followed by influx of extracellular Ca^{2+} through the voltage-activated Ca^{2+} channels, causing contraction of the smooth muscle cells (11). The central role of voltage-activated Ca^{2+} channels is further emphasized by the fact that the TGF response, the myogenic response, and autoregulation of renal blood flow all are abolished by blockers of voltage-activated Ca^{2+} channels (7, 30, 33, 34). Because entry of extracellular Ca^{2+} through voltage-activated Ca^{2+} channels seems to be a shared mechanism in both TGF and the myogenic mechanism, we chose to model the effect of TGF as a direct activation of voltage-activated Ca^{2+} channels. We used a single parameter to represent the entirety of the signaling process from the macula densa to the vascular smooth muscle cells. This is in effect equivalent to linearizing it, which appears to be justified because a more detailed description introduces a level of complexity inconsistent with all other aspects of the model, and also because no dynamic data are available to justify anything more.

The effects of TGF output and coupling parameter cannot be separated in our model because they multiply each other. One may therefore imagine that the predicted action of TGF represents signaling from the macula densa to receptors on smooth muscle cells, whereas the coupling coefficient expresses the sensitivity of intracellular coupling mechanisms to input to those receptors.

We have earlier used a form of the spatially extended model with a linear coupling of TGF to the diameter of the afferent arteriole rather than an explicit representation of the myogenic mechanism (17). This simplified model simulates the TGF oscillation but could not be made to bifurcate to chaos over a larger range of TGF gain values than has been observed experimentally, a finding consistent with the idea that interactions between the two oscillators are required in single nephron simulations.

We have also observed that nephrons communicate with each other using vascular signals initiated by TGF (29). The strength of this coupling is increased in hypertension (5, 47). We simulated a group of nephrons of different lengths, coupled by vascular signals (27). The nephron model (3) was much simpler than any we have used in this study and lacked a myogenic mechanism. Using parameter values that support a limit cycle oscillation with a single nephron version of this model, we found that many of the nephrons in this multinephron model operated in a chaos-like mode even at normal blood pressures. Increasing the coupling between nephrons increased the number of nephrons behaving chaotically. Layton et al. (20) and Pitman et al. (37), using their nephron model as the basis for a two-nephron configuration, found a range of values for a coupling term that generated chaos-like activity. Layton et al. (19) subsequently studied the effects of coupling two nephrons of different lengths on system dynamics, using their model of the thick ascending limb and TGF; nephron models were coupled linearly. They found that the two-nephron model exhibited limit cycle oscillations but could produce more complex dynamics under some combinations of nephron lengths. This result is quite similar to our finding with a multinephron model using a simpler nephron model than the one in this paper (27). Each of the simpler nephron models had a linear representation of vascular coupling between the macula densa and the afferent arteriole. The results of the current study show that the bifurcations of interest can also be caused in a single nephron model by an increase in value of the coefficient coupling TGF and the myogenic mechanism. All nephrons in this study were of identical lengths, a condition that, by itself, could not cause any of the bifurcations. The effects of the interplay between intranephron synchronization of TGF and the myogenic mechanisms, and internephron synchronization, remain to be studied.

Experimental results in two forms of experimental hypertension show irregular fluctuations in tubular and vascular dynamics (14, 15, 23, 51). The phase space attractor has been found to be low dimensional and to have a noninteger value, and the first Lyapunov exponent was positive (51). These are characteristics of datasets from chaotic systems. Surrogate data analysis confirmed the nonlinearity of the attractor (52). The ability to detect and measure chaos with these measures may be reduced when they are applied to experimental data. The principal reason for this uncertainty is that time series from real systems are affected by random processes, making it impossible to eliminate system and measurement noise. The animals used in those studies constitute complex systems from which it is technically impossible to eliminate all variation. Thus the experimental results can only be regarded as suggestive of chaos, and not a definitive demonstration.

An alternative to defining the dynamics generated by an experimental system is to simulate it, so that all sources of noise are eliminated and the results are strictly deterministic, the approach we used in this study. The models necessarily contain a number of approximations to the live kidney, however, and we know of no way to validate them completely. Moreover, the classic work on chaotic systems in the physics and mathematics literature generally uses third-order systems, whereas our models are much higher order, and the real kidney still higher. These considerations make it unlikely that a completely convincing case can be made that the observations of this study are due to chaos as rigorously defined.

The value of the approach we have followed, however, lies in the fact that chaotic systems have characteristic behaviors that differ from those of limit cycle oscillators, and in particular from those of linear systems, and these behaviors provide the basis for predictions about the performance of the real kidney. The task at hand then becomes one of predicting effects of nonlinear behavior in the kidney, and mapping these effects on to problems of physiological interest. The TGF and myogenic oscillations occur because those systems are nonlinear, at least second order, and contain dissipative terms, which are necessary conditions for a self-sustained oscillation to occur (2, 31, 32, 36). Whether systems with such properties oscillate is a matter of parameter values, and it would appear that the parameters in the kidney support these oscillatory dynamics. Because both TGF and the myogenic mechanism converge on a single contractile mechanism in arteriolar smooth muscle cells, an interaction is inevitable and synchronization is a likely outcome of that interaction.

Nephron synchronization is another effect of nonlinear interactions. This effect is well established in experimental measurements (12, 45, 49) and has been simulated with the same two-nephron model as used in this study (29) and with other simpler nephron models (4, 27, 37).

Nephron synchronization is impaired in hypertension (45). This effect may be a consequence of the bifurcation to the chaos-like state we have simulated in this study, or it may be a result of increased vascular signal coupling (5, 47), as we have also suggested (27), or both, and one change may be a response to the other.

### Perspectives and Significance

We simulated tubular pressure, flow, and NaCl concentration in spatially extended mathematical models of one and two nephrons. Each nephron model included an active myogenic mechanism and an interaction between TGF and the myogenic mechanism operating at the Ca^{2+} conductance of afferent arteriolar smooth muscle cells. Nephrons were coupled electrotonically in the two-nephron simulation. The tubular variables and nephron blood flow oscillate in normotensive rats but are irregular in hypertensive animals, in which TGF gain is increased. We varied the strength of coupling between TGF and the myogenic mechanism to test whether an increase in this parameter could cause the change in the dynamics, and whether the experimental results represented deterministic chaos. Solutions of both one- and two-nephron models developed the irregularity at higher coupling strengths, but the range of coupling strengths that could cause the irregularity was larger in the two-nephron model.

These results confirm that the observed increase in coupling between TGF and the myogenic mechanism can cause the irregularities seen in experimental data from hypertensive rats and that the change in dynamics represents a bifurcation to deterministic chaos. Nephrons are normally coupled electronically to their neighbor nephrons, so the two-nephron result is more likely to represent the system dynamics in vivo.

The oscillations in normotensive animals show extensive synchronization, but this phenomenon is decreased in hypertensive animals, reflecting the operation of a chaotic system. The presence of periodic oscillations is likely to lead to the formation of ensembles of synchronized nephrons, and the absence of these oscillations in hypertensive animals can be expected to reduce the size of such ensembles, and to reduce or eliminate the functional advantages such aggregation is likely to produce.

The strength of internephron coupling is also increased in hypertension, a change that models suggest can produce a bifurcation to chaos by itself. How the interaction between TGF and the myogenic mechanism affects the interaction between nephrons, and the effect of changes in both coupling strengths, remain to be studied.

## GRANTS

The work of D. J. Marsh was supported by National Institutes of Health Grant EB-003508 and by a grant from the Lündbeck Foundation of Copenhagen, Denmark. O. V. Sosnovtseva and N.-H. Holstein-Rathlou received support from Det Frie Forskingsråd. J. Laugesen, O. V. Sosnovtseva, E. Mosekilde, and N.-H. Holstein-Rathlou acknowledge support from the European Union via the European Network of Excellence Biosim, Contract No. LSHB-CT-2004-005137.

## DISCLOSURES

No conflicts of interest are declared by the authors.

- Copyright © 2010 the American Physiological Society