## Abstract

Single-nephron proximal tubule pressure in spontaneously hypertensive rats (SHR) can exhibit highly irregular oscillations similar to deterministic chaos. We used a mathematical model of tubuloglomerular feedback (TGF) to investigate potential sources of the irregular oscillations and the corresponding complex power spectra in SHR. A bifurcation analysis of the TGF model equations, for nonzero thick ascending limb (TAL) NaCl permeability, was performed by finding roots of the characteristic equation, and numerical simulations of model solutions were conducted to assist in the interpretation of the analysis. These techniques revealed four parameter regions, consistent with TGF gain and delays in SHR, where multiple stable model solutions are possible: *1*) a region having one stable, time-independent steady-state solution; *2*) a region having one stable oscillatory solution only, of frequency *f*_{1}; *3*) a region having one stable oscillatory solution only, of frequency *f*_{2}, which is approximately equal to 2*f*_{1}; and *4*) a region having two possible stable oscillatory solutions, of frequencies *f*_{1} and *f*_{2}. In addition, we conducted simulations in which TAL volume was assumed to vary as a function of time and simulations in which two or three nephrons were assumed to have coupled TGF systems. Four potential sources of spectral complexity in SHR were identified: *1*) bifurcations that permit switching between different stable oscillatory modes, leading to multiple spectral peaks and their respective harmonic peaks; *2*) sustained lability in delay parameters, leading to broadening of peaks and of their harmonics; *3*) episodic, but abrupt, lability in delay parameters, leading to multiple peaks and their harmonics; and *4*) coupling of small numbers of nephrons, leading to multiple peaks and their harmonics. We conclude that the TGF system in SHR may exhibit multistability and that the complex power spectra of the irregular TGF fluctuations in this strain may be explained by switching between multiple dynamic modes, temporal variation in TGF parameters, and nephron coupling.

- renal hemodynamics
- hypertension
- mathematical model
- nonlinear dynamical system

the tubuloglomerular feedback (TGF) system is a key regulator of single-nephron glomerular filtration rate (SNGFR) and of water and electrolyte delivery to the distal nephron. In the 1980s, experiments in rats by Leyssac and colleagues (16, 31) demonstrated that nephron flow and related variables may exhibit regular oscillations with a period of ∼30 s. Figure 1, *A1* and *A2*, illustrates their findings. Experiments have shown that these regular oscillations are TGF mediated (14, 17, 31). TGF, a negative feedback loop with time delays, has a resonant frequency: mathematical models have indicated that if the feedback is sufficiently strong (i.e., the feedback has sufficiently large gain magnitude), then the stable behavior of the TGF system is a self-sustaining regular oscillation and not a time-independent steady state (SS^{1}) (18, 26, 29). Such a regular oscillation is called a limit-cycle oscillation (LCO). When nephron parameters change, and thereby change TGF loop delays and gain, a bifurcation may occur; i.e., there may be a parameter-dependent change from one type of stable solution (e.g., a SS) to another type of stable solution (e.g., a LCO).

Spontaneously hypertensive rats (SHR) exhibit irregular TGF-mediated oscillations that appear to have characteristics consistent with deterministic chaos (16, 52); a typical recording of these irregular oscillations is illustrated in Fig. 1*B1*. Despite extensive investigation, the origins and consequences of the irregular oscillations remain obscure. Mathematical models that include spatially distributed NaCl transport along the thick ascending limb (TAL) have not predicted the irregular oscillations (18, 26). In contrast, single- or multinephron models that use a less realistic TAL representation, but a more extensive vascular representation, have exhibited irregular oscillations and other complex behaviors such as antiphase coupling, period-doubling, and the emergence of chaotic oscillations by a period-doubling route. However, the chaotic oscillations, in general, are only seen in parameter regimes that lie beyond measured values of TGF gain (3, 20, 42, 44, 45).

We have recently developed a systematic method for investigating, in the context of our model of the TGF system, the impact of NaCl backleak from the interstitium into the TAL. In the absence of backleak, NaCl concentration in the model TAL has been shown to depend on TAL transit time only (28, 39), but the addition of backleak introduces a significant degree of spatial dependence (26, 27). Our investigation, described herein, revealed unexpectedly rich dynamics in a parameter regime that appears to correspond to that reported in SHR. This finding motivated us to reexamine the following question: to what extent can the irregular TGF oscillations in SHR be explained by the intrinsic dynamic characteristics of the TGF system in a single nephron or in a small number of coupled nephrons? Our strategy is to identify unusual behaviors in SHR and key features of the complex power spectra that are computed from the irregular oscillations. Then, we will use our mathematical model to identify factors that may contribute to that spectral complexity.

Many of the key features of TGF-mediated oscillations in SHR, which we seek to explain, are illustrated in Fig. 1*B2*. This figure, a power spectrum, shows the frequencies, and relative strengths, of the sinusoidal components that, when summed together, make up the long-time record of the waveshape shown in Fig. 1*B1*. In the spectrum, the numbers above the peaks show the frequencies of the sinusoidal components (in mHz). Careful analysis of this spectrum suggests that elements of order reside amidst the appearance of chaos. First, several strong peaks are found in the range (25–50 mHz) that is typical of the fundamental frequency of TGF-mediated oscillations in normotensive rats (14, 16, 17). Second, some of the peaks appear to be closely spaced doublets or triplets; candidates are seen at 13.5 and 16 mHz; 18.5 and 22 mHz; 27 and 30 mHz; 30 and 33 mHz; 30, 33, and 36.5 mHz; 40 and 44 mHz; 49, 51.5, and 54 mHz; and 72 and 76 mHz. Third, some peak frequencies may be harmonics of lower fundamental frequencies (i.e., higher frequencies that are integer multiples of a lower frequency). Examples of frequencies that appear to be in harmonic relationship are 13.5, 27, 40, 54, and 67 mHz; 16, 33, 49, and 65 mHz; 18.5, 36.5, and 54 mHz; 22, 44, and 65/67 mHz; 33 and 65/67 mHz; and 36.5 and 72 mHz. Another published SHR power spectrum (Fig. 1*E* in Ref. 52) shows two major peaks that appear to be in harmonic relationship: one apparent fundamental oscillation at ∼23 mHz (somewhat lower than typical in control rats) and a second oscillation at ∼47 mHz (somewhat higher than typical in control rats). Fourth, peaks in the SHR power spectrum are often broad compared with those in spectra from normotensive rats. Examples in Fig. 1*B2* include peaks in the following ranges: 11.5–17.5 mHz; 17–26 mHz; 28–38 mHz; 38–45 mHz; and 45–56 mHz. Finally, more spectral power is found at very low frequencies (less than ∼20 mHz) in SHR than in the Wistar-Kyoto (WKY) rat controls (compare Fig. 1, *A2* and *B2*; see also Fig. 1*E* in Ref. 52 and ⇓⇓⇓⇓⇓⇓Fig. 8 in Ref. 54). Fluctuations in this frequency range are likely to arise from sources extrinsic to the TGF system, such as shifts in the neurohumoral systems that regulate renal vascular and tubular function.

An unusual behavior we seek to explain is that the regular oscillations (i.e., LCO) found in controls have not been reported in SHR; rather, nephrons in SHR appear to exhibit either steady flows or highly irregular oscillations (16, 52). Thus the development of irregular oscillations in SHR does not seem to follow the period-doubling route to deterministic chaos (36). Finally, a significant characteristic of SHR is that they have higher TGF gain (11, 32, 46) and stronger internephron coupling than normotensive controls (7, 49). Indirect evidence suggests that SHR may also have elongated TGF system delays because of reduced proximal tubule flow relative to controls (32, 46). Holstein-Rathlou and colleagues (19) have hypothesized that differences in key physiological parameters, especially those that have been identified as bifurcation parameters (e.g., TGF gain or nephron-to-nephron coupling strength) may explain, at least in part, the irregular oscillations in SHR.

In this study, we examine four factors that may contribute to spectral complexity. The first is the intrinsic complex dynamics of our model TGF system in a parameter regime of high gain consistent with measurements in SHR (11, 32, 46). In this region, we will show that TGF exhibits multistability, in that perturbations can produce LCO at one of two frequencies, a frequency of *f*_{1} and a second frequency of *f*_{2} that is approximately equal to 2*f*_{1} (however, *f*_{2} is not a harmonic of *f*_{1}). Moreover, parameter changes can produce a bifurcation directly from a SS to the LCO of frequency *f*_{2}. The second factor is alterations in TGF system time delays, which are determined by nephron flow rates, nephron dimensions, and the time required for propagation of the TGF response through the juxtaglomerular apparatus (JGA). A previous analysis of our model predicted that TGF time delays are important determinants of TGF dynamics (26). In this study, we present a more extensive bifurcation analysis in which steady-state TAL flow rate (or, alternatively, TAL volume) is a bifurcation parameter. The results suggest that elongated time delays, and especially TAL transit time, can introduce substantial spectral peaks below the ∼30-mHz TGF fundamental seen in normotensive rats (see below, Fig. 5*B*).

The third factor is vascular coupling between nephrons. Experimental studies have shown that nephrons sharing a common feed artery are functionally coupled, in that TGF responses in one nephron will influence the vascular resistance of closely coupled nephrons, leading to synchronous oscillations (7, 14, 22, 53). Finally, our simulations show that gradual shifts in system time delays within a bifurcation region lead to broadened spectral peaks, whereas abrupt changes can produce closely spaced doublets, multiple peaks, and frequency-shifted harmonic series.

Taken together, the results of our study provide a possible explanation for irregular TGF-mediated oscillations in SHR that does not require the postulation of deterministic chaos. Rather, the irregular oscillations in SHR and the associated complex power spectra may be consequences of intrinsic multistability of the TGF system in the parameter regime of the SHR, switching between oscillatory modes, strong nephron coupling between a few nephrons, and temporal variation in key parameters.

## MATHEMATICAL MODEL

### GLOSSARY

#### Independent Variables

*t*- Time (s)
*x*- Axial position along TAL (cm)

#### Specified Functions

- C
_{e}(*x*) - Extratubular [Cl
^{−}] (mM) - ψ
_{δ}(*t*) - Shape function for distributed delay (dimensionless)

#### Dependent Variables

- C
_{i}(*x*,*t*) - TAL [Cl
^{−}] (mM) - C(
*t*) - Effective MD [Cl
^{−}] (mM) - F
_{i}(*t*) - TAL fluid flow without coupling (nl/min)
_{i}(*t*)- TAL fluid flow with coupling (nl/min)
- S(
*x*) - Steady-state TAL [Cl
^{−}] (mM)

#### Parameters

*A*_{o}- Base-case TAL cross-sectional area, π
*r*^{2}(μm^{2}) - c
- Base-case TAL circumference, 2π
*r*(μm) - C
_{o} - [Cl
^{−}] at TAL entrance (mM) - C
_{op} - Steady-state [Cl
^{−}] at MD (mM) - F
_{op} - Base-case steady-state TAL flow rate (nl/min)
*k*_{i}- Sensitivity of TGF response (1/mM)
*K*_{M}- Michaelis constant (mM)
*L*- Length of TAL (cm)
*P*- TAL chloride permeability (cm/s)
- Q
_{op} - Steady-state SNGFR (nl/min)
- ΔQ
- TGF-mediated range of SNGFR (nl/min)
*r*- Luminal radius of TAL (μm)
*t*_{o}- Base-case steady-state TAL transit time (s)
*V*_{max}- Maximum transport rate of chloride (nmol·cm
^{−2}·s^{−1}) *K*_{1}- ΔQ/2Q
_{op}(dimensionless) *K*_{2}*k*_{i}C_{0}/2 (dimensionless)- α
- Fraction of SNGFR reaching TAL (dimensionless)
- γ
_{i} - Feedback gain magnitude, −
*K*_{1}*K*_{2i}*S*′(1) - γ
_{c} - Critical feedback gain magnitude (dimensionless)
- δ
- Distributed delay interval at the JGA (s)
- η
_{i} - TAL transit time scaling (dimensionless)
- τ
_{i} - τ
_{pi}+ δ/2, effective delay interval at the JGA (s) - τ
_{pi} - Discrete delay interval at the JGA (s)
- φ
_{ij} - Coefficient for coupling between nephrons
*i*and*j*

### Formulation

The mathematical model, which is similar to that previously described (26, 29), represents the TGF system in a short-looped rat nephron. The TAL is modeled as a rigid flow tube, and the actions of the afferent arteriole (AA), glomerulus, proximal tubule, and descending limb are represented by phenomenological relations. The only solute represented in the model TAL is the chloride ion (Cl^{−}), which is believed to be the principal signaling agent for TGF activation along the macula densa (MD) (41). The model has been expanded to represent coupling, through the preglomerular vasculature, with the TGF systems of other model nephrons (each labeled by an index *i*). Also, the model now includes a parameter η_{i}, for each model instance, that scales TAL cross-sectional area. Model variables and parameters, with their dimensional units as commonly reported, are given in the glossary. A schematic diagram for the *i*th model TGF system is given by Fig. 2. The spatial variable *x*, which designates position along the TAL, is oriented so that it extends from the TAL entrance (*x* = 0), through the outer medulla, and through the cortex to the site of the MD (*x = L*).

In the *i*th model the TGF system, TAL tubular fluid chloride concentration C_{i}(*x*,*t*), TAL tubular fluid flow , coupling with other nephrons, and a TGF signal delay at the JGA are represented by the following equations (1) (2) (3) (4) and (5) *Equation 1* expresses Cl^{−} conservation along the TAL; the first term on the right corresponds to axial intratubular solute advection; the two terms inside the large pair of parentheses correspond to active solute transport characterized by Michaelis-Menten-like kinetics and transepithelial diffusion characterized by backleak permeability *P*. At the entrance of the TAL (*x* = 0), we assume that the Cl^{−} concentration varies little; thus C_{i}(0,*t*) is set to a constant C_{o}.

*Equation 2* gives _{i}(*t*) as the sum of a nephron's own TGF-mediated TAL tubular flow response F_{i}(*t*) and a term from the other model TGF systems (as in Ref. 38); each summand in that term is a weighted difference between another system's TGF response F_{j}(*t*) and the time-independent steady-state flow αQ. The weight, or coupling coefficient, φ_{ij} is a nonnegative constant parameter that scales the influence of the *j*th nephron on the *i*th nephron.

*Equation 3*, the feedback response for an individual nephron, is an empirical relationship that has been well established by steady-state experiments (41). The constant C_{op} is the time-independent steady-state TAL tubular fluid Cl^{−} concentration at the MD when _{i}(*t*) assumes a steady-state value of αQ;C(*t*) represents the effective TAL tubular fluid Cl^{−} concentration after a discrete and distributed time delay described by *Eq. 4*. The final term in *Eq. 3*, ρ_{i}(*t*), represents transient perturbations that were introduced to test the stability of solutions.

Dynamic experimental results (5) indicate that a change in MD concentration does not significantly affect AA muscle tension until after a discrete (or “pure”) delay time τ and that the subsequent effect is temporally distributed so that a full response requires additional time, with the greatest weight in a time interval [*t*−τ−*δ*, *t*−τ], where δ is a second delay parameter. The convolution integral in *Eq. 4* describes the distributed delay signal received by the AA at time *t*; C_{MDi} is an “effective” chloride concentration at the MD that serves to transmit that delay to the TGF response represented in *Eq. 3*. *Equation 5* provides a weighting for a sigmoidal distribution of the distributed delay; a sigmoidal distribution is consistent with extant experimental data (Fig. 6 in Ref. 5).

For every choice of parameters within the physiological regime, the model equations have a time-independent steady-state solution; however, this solution may be stable (hence, a SS solution), or unstable, in which case an oscillatory solution is stable. The steady-state Cl^{−} profile along the TAL, denoted S(*x*), can be found by setting to the steady-state TGF operating point flow (F_{op}, which equals αQ) and solving the time-independent form of *Eq. 1* obtained when the term containing the time derivative is set to zero; the resulting equation is (6) where S(0) = C_{o}, the chloride concentration entering the TAL, is the initial value of S. Because *Eq. 6* does not depend on η_{i}, neither does S(*x*); therefore S(*x*) is the same for each model TGF system *i*.

The dimensionless parameter η_{i}, which scales TAL cross-sectional area, has an important alternative interpretation: one can consider η_{i} to be a simultaneous rescaling of the steady-state flow rate and of the transport parameters *V*_{max} and *P*. Thus one can rewrite *Eq. 1* as (7) which can be considered to have new values for flow rate, maximum transport rate, and backleak permeability: (*t*)/η_{i}, *V*_{max}/η_{i}, and *P*/η_{i}, respectively. Notwithstanding the reinterpretation of η_{i}, *Eqs. 7* and *1* have exactly the same solution.The equation for the steady-state solution corresponding to *Eq. 7* is (8) Because the η_{i}'s in *Eq. 8* cancel, *Eqs. 8* and *6* are exactly the same equation and have exactly the same solution S(*x*).

Thus one can interpret η_{i} in two ways. One can consider η_{i} to be a TAL volume scaling (relative to the base case of η_{i} = 1) that represents a uniform inflation (η_{i} > 1) or deflation (η_{i}< 1) of the TAL lumen all along the *i*th TAL, relative to a reference cross-sectional area *A*_{o}; if η_{i} = 1, the cross-sectional area is considered to be *A*_{o}. The intratubular perimeter length corresponding to the reference area *A*_{o} is denoted c; that perimeter length is considered to be fixed, so that transport capacity per unit TAL length remains constant when η_{i} is varied. Thus although the cross section may vary from simulation to simulation (or vary as a function of time), total TAL transport capacity is unchanged. This interpretation of η_{i} would be appropriate to describe the effects of a short-term change in intrarenal pressure that changes tubular volume without changing TAL transport capacity.

Alternatively, one can consider η_{i} to be a scaling of the TAL flow rate in which the transport characteristics of the TAL have been also scaled in such a way that the same steady-state concentration profile S(*x*) is preserved. This situation might arise from a long-term adaptation of the TAL to a change in flow rate, as, for example, from a long-term change in either SNGFR or fractional absorption from the proximal tubule.

Both interpretations of η_{i} can be unified under the concept of the steady-state transit time interval for the TAL; this interval is the time required for a water molecule to travel up the TAL, from the TAL entrance to the MD, at a constant speed, assuming plug flow. For our steady-state base-case parameters, that transit time, denoted *t*_{o}, equals the TAL volume, *A*_{o}*L*, divided by the steady-state flow F_{op}. However, for a change in η_{i}, under either of the interpretations advanced above, the steady-state transit time becomes η_{i}*t*_{o} and thus η_{i} may be considered to be a dimensionless scaling of the base-case steady-state TAL transit time *t*_{o}.

### Parameters

Base-case values for model parameters are given in Table 1. Extratubular concentration is specified by C_{e}(*x*) = C_{o} [*A*_{1} exp(−*A*_{3} *x*) *+ A*_{2}], where *A*_{1} = (1 − C_{e}(*L*)*/*C_{o})/(1 − exp(−*A*_{3})), *A*_{2} = 1 − *A*_{1}, and *A*_{3} *= 2*, and where C_{e}(*L*) corresponds to a cortical interstitial concentration of 150 mM. Graphs for C_{e}(*x*) and the steady-state profile S(*x*) were given in Fig. 1 of Ref. 26. The steady-state concentration C_{op} corresponds to S(*L*).

In the base case, we consider a discrete delay τ of 2 s to be followed by a transition interval of 3 s, providing an overall delay τ_{i} = τ+ δ/2 of 3.5 s, in approximate agreement with experiments (5). The overall delay τ_{i} was varied by varying τ, whereas δ was considered to be a fixed parameter in all cases (see supplementary material, item 1; all supplementary material in this paper can be found at http://ajprenal.physiology.org/cgi/content/full/00048.2005.DC1).

### Analysis of Model Behavior

We used two methods to investigate model behavior: *1*) analysis of the model's characteristic equation and *2*) direct computation of numerical solutions of the model equations. The characteristic equation contains information that allows systematic predictions of generic model behavior as a function of model parameters. By “generic behavior,” we mean, e.g., whether the model's long-time behavior will approximate a time-independent, stable steady state (i.e., a SS) or a sustained oscillation. A numerical solution of the model equations allows one to determine model behavior for a specific set of parameters. These two methods are complementary, in that the first provides helpful, but not always certain, predictions for many parameter combinations, whereas the second is a more nearly certain predictor for behavior, but only for a particular set of parameters. Moreover, some features of the predictions by both methods should be entirely consistent, so each method provides a validation of the other. The methods that were used to compute numerical solutions from the model equations are summarized in the appendix; the appendix also contains a synopsis of the procedures that were used to confirm or identify bifurcation boundaries.

For this study, we used only the characteristic equation for a single, uncoupled nephron (φ_{i,j} = 0 for all *j* in *Eq. 2*). We have previously described how such a characteristic equation is obtained (26). Briefly, the nonlinear model equations (*Eqs. 1*–*5*) are combined into a single nonlinear equation and linearized to obtain a linear partial differential equation (PDE). An expression representing potential steady-state and oscillatory solutions is substituted into that linear PDE, and from the resulting equation, the conditions for the potential solutions that satisfy the linearized equation may be derived. Those conditions are embodied in the characteristic equation, which is given by (9) The index *i* has been dropped from the model parameters, because it is not relevant for the uncoupled case; also, all variables and parameters appearing in *Eq. 9* have been nondimensionalized as described in the appendix. In particular, τ_{p} and δ have been rendered dimensionless through division by the base-case steady-state TAL transit time *t*_{o}. However, η, already dimensionless, is unchanged. The additional parameter γ appearing in *Eq. 9* is the feedback-loop gain magnitude, a dimensionless quantity that is the product of two factors (26, 27): *1*) the dimensionless slope of the TGF response curve at the steady-state operating point (F_{op}, C_{op}), which is given in nondimensional form by *K*_{1}*K*_{2} (see appendix); and *2*) the magnitude of the nondimensional slope of the TAL tubular fluid Cl^{−} concentration profile alongside the MD [i.e., −S′(*L*), but expressed in nondimensional form]. The two factors that make up γ depend on the model parameters. In this study, *K*_{1} and −S′(*L*) are assumed to have fixed values; to obtain particular values of γ, *K*_{2} is varied.

A solution of the characteristic equation is a complex number λ that depends on the values assigned to the model parameters. Because λ is a complex number, it can be written in the form where Re λ is the real part of λ and Im λ is the imaginary part of λ. The real and imaginary parts of λ are indicators of the strength and angular frequency, respectively, of a solution to the linearized model equations. A solution to the characteristic equation is not unique; rather it is a number appearing in an infinite series of solutions λ_{1}, λ_{2}, λ_{3}, …. However, as shown in results, only λ_{1} and λ_{2} appear likely to have a tangible impact on TGF system behavior. Because a solution λ of the characteristic equation may reasonably be considered to be an eigenvalue of a particular differential operator (34), we will refer to such solutions as eigenvalues (see supplementary material, item 2).

When backleak permeability *P* is set to zero, the characteristic equation given in *Eq. 9* has a simpler form, given by (10) The two forms of the characteristic equation (*Eqs. 9* and *10*) were solved by the systematic application of a method previously described (appendix c in Ref. 27).

The characteristic equation makes predictions about the stable solutions of the linearized model equations, as a function of parameter values (26). If, for a particular set of parameters, Re λ_{n} < 0 for all indices *n*, then the only stable solution is an SS. If, on the other hand, Re λ_{n} > 0 for one, and only one index *n*, the only stable solution is a regular oscillation having frequency, for that index *n*, of ∼(Im λ_{n})/(2π*t*_{o}). For cases where Re λ_{n} > 0 for multiple values of the index *n*, there may be multiple stable oscillatory solutions. In such case, the nature of the solutions must be investigated directly by solving the model equations to determine potential solution behaviors, behaviors that may vary for differing initial conditions. The value of the gain γ that corresponds to the transition from a case where Re λ_{n} < 0, for all indices *n*, to a case where, for at least one index *n*, Re λ_{n} > 0, is a called a critical gain, denoted γ_{c} (thus for that *n*, Re λ_{n} = 0 at the point of transition). For the base-case parameters, γ_{c} = 3.24 where Re λ_{1} = 0. (For additional information regarding the applicability and limitations of the characteristic equation, see supplementary material, item 3).

## RESULTS

### Backleak Permeability May Be an Important Bifurcation Parameter

As we now demonstrate, for the case of nonzero TAL backleak permeability, the qualitative behavior of our TGF system model is very different from that previously inferred by us from the case of no backleak permeability (26).

The diagrams in Fig. 3, *A1* and *B1*, for a parameter region in the the τ-γ plane, are based entirely on our analysis of the characteristic equations *Eqs. 10* and *9*, for no backleak (*P* = 0) and for nonzero backleak (*P* = 1.5 × 10^{−5} cm/s), respectively. In both cases, steady-state TAL transit time was fixed at its base-case value, thus η = 1. Throughout results, τ (or τ_{i}), which characterizes the TGF signal delay at the MD, designates the sum of the nondimensional values of τ_{p} (or τ) and δ (which correspond to the durations of pure delay τ_{p} and a distributed delay δ, as reported in Ref. 5). The gain magnitude γ (or γ_{i}) characterizes the strength of the TGF response (see mathematical model), as in experimental studies (15, 47). Solid curves in these diagrams mark loci where the real part of an eigenvalue λ_{n} is zero and the sign of the real part changes from negative to positive, as indicated by the inequalities in the figures. Dashed curves mark loci where the real part of an eigenvalue overtakes the real part of another eigenvalue. Thus, for example, for values of τ and γ that fall above the curves marked “Re λ_{2} = Re λ_{1},” Re λ_{2} exceeds Re λ_{1}.

Although, as will be apparent below, the diagrams in Fig. 3, *A1* and *B1*, provide a good guide to the generic solutions of our model equations, they are nonetheless based on the linearized forms of those equations, and detailed model behavior must be investigated and confirmed by means of simulations in which numerical solutions to the unaltered model equations are computed. Figure 3, *A2* and *B2*, which portray the same region of the τ-γ plane as do Fig. 3, *A1* and *B1*, indicate the behaviors of model solutions based on numerical solutions, which were computed as described in the appendix. Figure *3*, *A2* and *B2*, which are bifurcation diagrams, show the true model behaviors, which are substantially consistent with the indications of Fig. 3, *A1* and *B1*.

Figure 3*A2*, for the no-backleak case, corresponds to Fig. 3*A1*. For this case, only two generic stable behaviors are observed: if the values of τ and λ fall below the curve labeled “Re λ_{1} = 0,” then, for any initial conditions, or for any transient perturbation of a steady-state solution, the model solutions always more and more closely approximate the SS solution. However, for values of τ and γ above the curve Re λ_{1} = 0, the only stable solution that was observed was a regular oscillation that converged to a fixed, periodic trajectory; the period and amplitude of that oscillation vary as a function of τ and γ. This oscillation, which is an LCO, has a frequency *f*_{1} that approximates the frequencies [viz., (Im λ_{1})/(2π*t*_{o})] that can be found by solving the characteristic equation along the locus where Re λ_{1} = 0. Because the LCO can be characterized by the frequency *f*_{1}, we call it an “*f*_{1}-LCO”; this LCO may be considered to be the eigenmode associated with the eigenvalue λ_{1}. The *f*_{1}-LCO solutions for SNGFR and for MD Cl^{−} concentration, at the point labeled “*A3*” in Fig. 3*A2* (τ = 0.135, γ = 8.00), are shown in Fig. 3*A3*. A large number of numerical simulations, using sundry perturbations (which included the transient addition of oscillations of various frequencies to TAL flow), revealed no other stable solution for the region above the loci of Re λ_{1} = 0.

However, when TAL backleak permeability *P* was increased from 0 to *P* = 1.5 × 10^{−5} cm/s, the diagram obtained by solving the characteristic equation (*Eq. 9*), and the observed model behavior, both changed dramatically, as shown in Fig. 3, *B1* and *B2*, which are analogous to Fig. 3, *A1* and *A2*. In Fig. 3*B1* the curves for the conditions Re λ_{2} = 0 and Re λ_{3} = 0 cross the curve for Re λ_{1} = 0. Our numerical simulations revealed that such crossings can have important consequences for the nature of the stable solutions that can occur when the real part of at least one eigenvalue is positive. As in Fig. 3*A1*, a bifurcation locus Re λ_{1} = 0 marks the transition from an SS to an *f*_{1}-LCO. However, this bifurcation locus is now restricted to the portion of the diagram that lies to the right of the crossing of Re λ_{1} = 0 and Re λ_{2} = 0. Moreover, numerical solutions revealed a curve, labeled “Empirical Boundary,” that lies above the curve for Re λ_{2} = Re λ_{1} and that marks the boundary between the *f*_{1}-LCO and a different type of dynamic behavior: in the parameter region between the Empirical Boundary and the curve for Re λ_{2} = 0, there is a potential stable LCO solution that corresponds to λ_{2}; we call that solution an “*f*_{2}-LCO.” The numerical solutions showed that an *f*_{2}-LCO has a frequency approximately twice that of an *f*_{1}-LCO. In the *top left* where only Re λ_{2} > 0, which is labeled by “*f*_{2},” the only stable solution is an *f*_{2}-LCO. Moreover, above the curves for the Empirical Boundary and for Re λ_{1} = 0, two stable LCO, of differing frequencies, can be elicited, depending on initial conditions. One of these solutions is an *f*_{1}-LCO, and the other is an *f*_{2}-LCO; we call this behavior an “*f*_{1}, *f*_{2}-LCO.” Thus this parameter region, where the real parts of both λ_{1} and λ_{2} are both positive, and where Re λ_{2} exceeds Re λ_{1}, exhibits bistability with respect to the *f*_{1} and *f*_{2} solutions. We did not observe any *f*_{3}-LCO. The region in which Re λ_{3} > Re λ_{2} is completely contained within the region where Re λ_{2} > 0, and thus the relationship of the real parts of these eigenvalues is analogous to that for the real parts of λ_{2} and λ_{1} in the case of *P* = 0 shown in Fig. 3*A1*. The corresponding solution behaviors for these analogous cases are consistent in that *f*_{2}-LCO were not observed in the parameter region of Fig. 3*A2*, and *f*_{3}-LCO were not observed in Fig. 3*B2*.

The gray disk, which has its center at values of τ ≈ 0.2228 and γ ≈ 3.24 (and that thus lies on the curve for Re λ_{1} = 0), represents an approximate region for typical parameters in normotensive rats, which, based on our previous studies, appears to be near the bifurcation boundary from the SS to *f*_{1}-LCO (26, 30).

LCO in SNGFR and tubular fluid Cl^{−} concentration at the MD that were computed for points *A3* and *B3* are shown in Fig. 3, *A2* and *B2*, respectively. The delay and gain for both points were set to (τ, γ) = (0.135, 8.00). In the zero backleak case, SNGFR and Cl^{−} concentration, shown in Fig. 3*A3*, both oscillate with an *f*_{1}-LCO frequency of ∼50.1 mHz. However, for the nonzero backleak case, the analogous point *B3* lies in the region for *f*_{1}, *f*_{2}-LCO; thus two different stable solutions are possible. In Fig. 3*B3*, we exhibit the *f*_{2}-LCO, which oscillates with a frequency of ∼101 mHz.

Three considerations may clarify these findings. First, in the regions where LCO are the stable solutions, the frequencies of those solutions vary as a function of τ, γ, and η. Moreover, within a region where a particular type of LCO is observed, its frequency is much more sensitive to variations in time delays, and thus to changes in τ and 1/η, than to changes in the gain γ. Thus, for example, an *f*_{1}-LCO has a frequency along a vertical line (i.e., a fixed value of τ) in Fig. 3, *A2* or *B2*, that closely approximates the frequency predicted by the characteristic equation at the point where that line crosses the bifurcation curve for Re λ_{1} = 0. An analogous result was found for *f*_{2}-LCO.

Second, although an *f*_{2}-LCO has a frequency about twice that of an *f*_{1}-LCO (provided that the values of τ and η in the *f*_{2} case do not differ or differ little from the corresponding values in the *f*_{1} case), these frequencies are not in precise harmonic relationship. Thus even though the frequencies *f*_{1} and *f*_{2} correspond to the first two eigenmodes for resonate standing waves within the TAL, they are differentially influenced by their respective bifurcations and by nonlinearities inherent in the full model equations.

Finally, for ease of comparison among our TGF investigations, we have used a uniform set of base-case model parameters (26–30, 35). These parameters, which were based on experimental evidence (26), yield an *f*_{1}-LCO at ∼45 mHz (29), which is in the upper frequency range of experimentally measured LCO oscillations in rats (19). However, our model framework is generic, in the sense that the results depend on dimensionless lumped parameters (e.g., τ, γ, and η), and differing combinations of those parameters predict oscillations of varying frequencies. Indeed, as shown in Fig. 5 (gray dashed lines) and exemplified in numerical results of Figs. 6 and 7, *f*_{1}-LCO frequencies can change substantially with changes in TAL transit time, which is characterized by η, and these frequencies encompass the physiological range.

### TAL Transit Time May Be an Important Bifurcation Parameter

We varied TAL transit time by varying η, which may be considered to scale either the TAL cross-sectional area (and hence TAL volume) or the base-case steady-state TAL flow F_{op} (see mathematical model). To investigate the impact of η in the context of the previously identified bifurcation parameters τ and γ, we extended the planar parameter region in Fig. 3 to three dimensions, with axes oriented as shown in the *inset* in Fig. 4*A1*. Figure 4, *A1* and *A2*, shows diagrams corresponding to a portion of the (1/η)-γ parameter plane for fixed feedback delay τ = 0.2228 (dimensional value, 3.5 s), and for *P* = 1.5 × 10^{−5} cm/s. These diagrams are analogous to, and share a qualitative similarity with, the τ-γ diagrams given by Fig. 3, *B1* and *B2*; the interpretation of Fig. 4, *A1* and *A2*, is therefore similar to the interpretation of Fig. 3, *B1* and *B2*.

Below the curves in Fig. 4, *A1* and *A2*, the only stable solution is a time-independent steady state (labeled “Re λ_{n}< 0” and “SS” in Fig. 4, *A1* and *A2*, respectively). In Fig. 4*A2*, the parameter regions labeled “*f*_{1}”, “*f*_{2},” and “*f*_{1} or *f*_{2}” have the same model solution properties as described for the analogously labeled regions in Fig. 3*B2*. In addition, however, Fig. 4*A2* contains a distinct region, labeled by “*f*_{3},” where Re λ_{3} is sufficiently large that an oscillation having a frequency of ∼3*f*_{1} is stable (note that in Fig. 4*A1*, the loci for Re λ_{3} and Re λ_{2} cross, so there is a region in which only *f*_{3}-LCO are stable, analogous to the region in Fig. 3*B2* where only *f*_{2}-LCO are stable).

These results indicate that η is also a bifurcation parameter, and that 1/η has effects similar to those of τ (compare Figs. 3*B2* and 4*A2*). Both τ and η can be considered delay parameters, because τ characterizes the delay in TGF signal transmission at the JGA whereas η is proportional to TAL transit time. However, we have previously shown, in the simpler model context of zero TAL backleak, that the bifurcation behavior depends on the ratio of the delay at the JGA to the TAL transit time (26); in the present context, that ratio corresponds to the ratio τ/η. This explains the strong similarities between the effects of changes in τ and of changes in 1/η.

### Multiple Stable Solutions and Multistability May Introduce Spectral Complexity

Figure 4, *B1* and *B2*, shows simulated oscillations of SNGFR and the associated power spectrum corresponding to *point B* in Fig. 4*A2*, where (1/η, γ) = (2/3, 9/2). (The power spectra reported in this study were produced by the same method as described in Ref. 29; however, no broadband perturbations were used in the present study.) Because *point B* lies below the curves, where the real parts of all eigenvalues are negative, the stable solution is a SS: a solution that is transiently perturbed will always return more and more closely to that steady state, provided that no further perturbations are applied. Such a return is an evolving process of successive approximation, like that exemplified in Fig. 4*B1*. The model's steady-state SNGFR is 30 nl/min, and the power spectrum for that steady state shows that all frequency components are zero (all power spectra presented herein are based on the deviation from the steady-state SNGFR). *Point C*, where (1/η, γ) = (2/3, 11/2), corresponds to a single stable solution, the *f*_{1}-LCO. Simulated SNGFR oscillations and the associated power spectrum for *point C* are shown in Fig. 4, *C1* and *C2*. The oscillation frequency is 35 mHz, and because the oscillation is nonsinusoidal it has harmonics (29), the first of which is at 70 mHz. *Point D*, at (1/η, γ) = (2/3, 13/2), lies in the bistable region of *f*_{1}, *f*_{2}-LCO. Simulated SNGFR oscillations (for the *f*_{2}-LCO) and the associated power spectrum are shown in Fig. 4, *D1* and *D2*. The frequency of this oscillation is 66 mHz, which is nearly 2 × 35 mHz. It is noteworthy that three very different stable model solutions exist for values of gain that span about 2/9 of the reported experimental range of 1.5 to 9.9 (15). Indeed, examination of Fig. 4*A2* (near *points B*, *C*, and *D*) shows that the three differing solution types (viz., SS, *f*_{1}-LCO, and *f*_{2}-LCO) can be obtained for an even smaller range of γ than was used in the example presented here. Hence, in this parameter region, even small stochastic variations in model parameters can result in switching between oscillatory modes.

It is noteworthy that *point D* in Fig. 4*A2* is in a bistable “*f*_{1} or *f*_{2}” region having two stable solutions, viz., an *f*_{1}-LCO and an *f*_{2}-LCO (an analogous region was shown in Fig. 3*B2*). The *f*_{1}-LCO at *point D* is very similar to the *f*_{1}-LCO found at *point C* (see Fig. 4*C1*); the *f*_{2}-LCO at *point D* was shown in Fig. 4*D1*. The two stable solutions at *point D* are an instance of multistability. Whether the *f*_{1}-LCO or the *f*_{2}-LCO emerges as the stable mode depends on initial conditions (see text describing Fig. 7 below); moreover, perturbations can cause switching from one mode to another. Because perturbations and parameter variations can occur within short time intervals both multistability and multiple stable solutions occurring within a limited parameter regime may result in multiple strong peaks in power spectra.

### Parameter Nonstationarity Can Introduce Spectral Complexity

Figure 5, *A* and *B*, shows τ-(1/η) bifurcation diagrams for γ = 5 and γ = 7, respectively; in both cases, *P* = 1.5 × 10^{−5} cm/s. These diagrams correspond to a portion of the τ-(1/η) plane, which is perpendicular to the γ-axis (see inset in Fig. 4*A1*).

In Fig. 5*A*, the region above the curve labeled by “Re λ_{1} = 0” corresponds only to stable solutions exhibiting an *f*_{1}-LCO; the region below the curve corresponds only to the time-independent steady-state solution, SS. The solid dot, at (τ_{o}, 1/η_{o}) = (0.2228, 0.7282), lies just within the region having an *f*_{1}-LCO solution. The three gray dashed curves mark the τ-(1/η) loci for which the frequencies of *f*_{1}-LCO are predicted by the characteristic equation to be ∼20, ∼35, and ∼50 mHz.

For fixed gain (γ = 5), time-dependent variations were introduced in the values of 1/η and τ. In the “η case,” τ was fixed at τ_{o}, but η was varied linearly, over an interval of 4,096 s (≈68 min), at a uniform rate of change, from 80 to 120% of η_{o}. Subsequently, in the “τ case,” η was fixed at η_{o}, but τ was varied linearly, over the same time interval, at a uniform rate of change, from 120 to 80% of τ. In Fig. 6, *A1* and *A2*, for linear and logarithmic ordinate scales, respectively, power spectra for these simulations (solid curves for η and dashed curves for τ) and for the stationary parameter case (gray curves) of (τ, η) = (τ_{o}, η_{o}), are shown. In the simulations with variation in η or τ, peaks in the power spectra are broadened; i.e., SNGFR oscillations contain a wider range of frequencies, relative to the stationary case. Peaks in the power spectra of the η case are slightly to the right of stationary case, because the higher frequencies correspond to larger 1/η's, which in turn correspond to parameter values farther away from the SS region; thus larger 1/η corresponds to oscillations of larger amplitude. In contrast, the peaks in the power spectra of the τ case are slightly to the left of the base case. In this case, the lower frequencies correspond to larger τ's and thus to oscillations of larger amplitude. (Note the good agreement between the frequencies in the power spectra, which were computed from the full model, and the frequencies predicted by the characteristic equation, as shown by the gray dashed curves in Fig. 5*A*.)

In another simulation, also for γ = 5, τ was fixed, but η was set to be 80% of η_{o} for approximately half of the time interval; then, in two nondimensional time units (2*t*_{o} = ∼31.42 s), η was increased linearly to 120% of η_{o}; and for the remainder of the time interval, η = 120% of η_{o}. The analogous simulation was also conducted for τ. Power spectra obtained for the η case and for the stationary case (τ, η) = (τ_{o}, η_{o}) are shown in Fig. 6, *A3* and *A4*, in linear and logarithmic ordinate scales, respectively. Two large peaks, rather than one broadened peak, were obtained for the η case; these peaks, at 43 and 37 mHz, correspond to η = 80% of 0 and 120% of η_{o}, respectively. The 43-mHz peak is larger because oscillations have larger amplitude for smaller η's. In contrast to the η arrow, which lies approximately perpendicular to the “∼35 mHz” curve, the τ arrow lies approximately parallel to the ∼35 mHz curve. Thus the two peaks obtained by varying τ (not shown) lie much closer to each other, compared with the η case.

The experiments described above were repeated for γ = 7; the bifurcation diagram for this case is shown in Fig. 5*B*, and corresponding results are in Fig. 7. Because the base parameters (dot, Fig. 5*B*) lie within the region of *f*_{1}-LCO or *f*_{2}-LCO, one of these two stable solutions had to be chosen; we chose the *f*_{2}-LCO of ∼70 mHz (gray curves, Fig. 7). However, because the *f*_{1}, *f*_{2}-LCO region is bistable, the mode corresponding to an *f*_{1}-LCO was also excited, as can be seen in the logarithmic power spectrum (gray curve, Fig. 7*B2*). As in the previous experiments, η and τ were varied linearly, at a constant rate of change, corresponding to paths along the arrows indicated in Fig. 5*B*, over a interval of ∼68 min. The tails of both arrows lie within the *f*_{1} region, whereas the heads lie within the bistable region. As shown in Fig. 7, *B1* and *B2*, *f*_{1}-LCO were obtained for both cases. As in the γ = 5 case, the peaks (∼37 mHz) were broadened, because the oscillations involve a range of frequencies, with the η peak wider than the τ peak. It is noteworthy that if η and τ were varied along different arrows such that the parameters began in the *f*_{2} region and moved into the bistable region, or if an *f*_{1}-LCO in the bistable region was perturbed using a sinusoidal wave having a frequency approximating the frequency of the *f*_{2}-LCO, the resulting solutions were *f*_{2}-LCO (results not shown).

In the next set of simulations, τ was set to be 120% of τ_{o} for approximately half of the time interval; then, in two nondimensional time units, τ was decreased at a constant rate to 80% of τ_{o}; and for the remainder of the interval, τ = 80% of τ_{o}. The power spectrum associated with the resulting oscillations has three prominent peaks (Fig. 7, *B3* and *B4*). The 35-mHz peak corresponds to the *f*_{1}-LCO obtained for τ = 120% of τ_{o}, which lies within the *f*_{1} region; the 70-mHz peak is its first harmonic. The 67-mHz peak, which is also excited, corresponds to the *f*_{2}-LCO obtained for τ = 80% of τ_{o}, which lies within the bistable region. A goal of this numerical experiment was to exhibit distinct peaks at approximately the fundamental frequency and the first harmonic, similar to the relationships seen in some power spectra from physiological experiments (see discussion). Thus it is desirable in the simulation for the fundamental frequency to remain relatively constant as the parameter is varied. Because, as we have seen, the same variations in τ (as a percentage of the base case) introduce smaller changes in the fundamental frequency than do variations in η, results are shown in Fig. 7, *B3* and *B4*, for the τ case rather than for the η case.

### Nephron-to-Nephron Coupling May Introduce Spectral Complexity

We investigated the effects of coupling a model nephron having parameters in the bistable (“*f*_{1} or *f*_{2}”) region with nephrons in the *f*_{1} region. The bifurcation conditions for a system of coupled nephrons differ somewhat from those for a single nephron (38). However, for simplicity, the bifurcation conditions for a single nephron were used in this study as a guide for predicting coupled behavior. First, two nephrons, designated *1* and *2*, were coupled. The parameters for the nephrons are given at the top of Table 2; the coupling parameter φ_{ij} was set to 0.2 for both nephrons. Figure 8, *A1*–*A3*, show oscillations in SNGFR for *nephron 1* and corresponding power spectra, in linear and logarithmic scales. In this nephron, the frequencies of ∼27 and ∼60 mHz were most excited. These frequencies correspond to the fundamental frequency of *nephron 2*, which is in the *f*_{1} region, and the *f*_{2}-LCO of *nephron 1*, which is in the bistable region. The fundamental frequency of *nephron 1* (33 mHz) and the first harmonic of *nephron 2* (54 mHz) were also excited.

In another simulation, the TGF systems of three model nephrons, designated *1*, *2*, and *3*, were coupled. The parameters for the nephrons are given at the bottom of Table 2; the coupling parameter φ_{ij} was set to 0.2 in all couplings among the nephrons. With these parameters, *nephron 1* is in the bistable region, whereas *nephrons 2* and *3* are in the *f*_{1}-LCO region. Figures 8B1-8B3 show oscillations in SNGFR for *nephron 1* and corresponding power spectra. In this case, the most excited frequencies are 27 mHz (*f*_{1}-LCO of *nephron 3*), 42 mHz (*f*_{1}-LCO of nephron 2), and 62 mHz (*f*_{2}-LCO of *nephron 1*). The *f*_{1} frequency of nephron 1 (34 mHz) and the first harmonics of *nephrons 2* and *3* (84 and 55 mHz) were also excited. The *f*_{1}-LCO's in *nephrons 2* and *3* do not synchronize because their delays are sufficiently different (23, 38). These results predict that irregular oscillations in SNGFR and complexities in corresponding power spectra can be obtained via internephron coupling.

### Combined Effects of Bifurcations, Coupling, and Parameter Lability

The results summarized in Fig. 8 show that a large measure of irregularity can be introduced into our model's TGF-mediated oscillation by internephron coupling among nephrons when both *f*_{1}-LCO and *f*_{2}-LCO are present. However, the degree of spectral complexity in Fig. 8, *B1* and *B2*, is much less than that shown in the experimental spectrum reproduced in Fig. 1*B2*. We now show that when a modest degree of temporal variation is introduced in η, model power spectra with features similar to Fig. 1*B2* can be produced.

We began with the model configuration for three coupled nephrons that was investigated above. Temporal variation was introduced in two stages. First, similar to the cases shown in Fig. 6, *A3* and *A4*, we introduced an abrupt change in 1/η_{i} (*i* = 1, 2, 3); this change shifted the LCO frequencies of each of the three coupled nephrons. Second, we introduced a sustained 10-mHz sinusoidal oscillation in 1/η_{i}, a perturbation that is consistent with the substantial spectral power that has been reported below 20 mHz in SHR (15, 52, 54).

The abrupt change in 1/η_{i} was performed as previously described for Fig. 6, except that 1/η_{i} was increased by 10% in each of the three coupled nephrons. Such a shift might arise from a sustained systemic decrease in blood pressure that results in decreased volume within the renal capsule. The results are shown in Fig. 9, *A1*–*A3*. Peaks previously identified in Fig. 8, *B2* and *B3*, are in lower-case italics; new peaks are in larger Roman font. Note that the number of significant peaks is essentially doubled, because the flow in each nephron first operated at the LCO frequency shown in Fig. 8*B2* and then at a shifted frequency. Thus *nephron 3* produced a spectral peaks of its *f*_{1}-LCO at both 27 and 29 mHz, and *nephron 2* produced analogous peaks at 42 and 45 mHz. *Nephron 1* produced spectral peaks of its *f*_{2}-LCO at 62 and 67 mHz. The *f*_{1} frequency for *nephron 1* was sympathetically excited at both 34 and 37 mHz.

The sustained sinusoidal oscillation of 10 mHz was then introduced in addition to the abrupt change; results are shown in Fig. 9, *B1*–*B3*. Peaks previously identified in Fig. 8, *B2* and *B3*, and Fig. 9, *A2* or *A3*, are in lower-case italics; new peaks are in larger Roman font. When the sinusoidal oscillation was introduced, the peaks were broadened (especially noticeable at the bases of the peaks), and their amplitude was significantly reduced (note differing vertical scales in Fig. 9, *A2* and *B2*). However, the major peaks identified in Fig. 9*A2* were all present in Fig. 9*B2*, although some were slightly shifted (e.g., from 34 to 35 mHz). Moreover, the addition of the sinusoidal oscillation not only introduced, as expected, a peak at 10 mHz, but other new peaks were also introduced. These new peaks may arise from simple mechanisms. For example, the peak at 59 mHz appears to be the first harmonic of the *f*_{1}-LCO of *nephron 3*, which has a frequency, to the nearest half-mHz, of 29.5 mHz. The peak at 32 mHz may be the result of the composite subharmonics of both *f*_{2}-LCO frequencies, 62 and 67 mHz. Oscillations at 65 mHz may be a harmonic of the 32-mHz oscillation. Alternatively, these peaks may arise through other mechanisms that have been recently identified in our models of coupled nephrons, e.g., spectral splitting (38) or heterodyning (23).

These final examples demonstrate that substantial spectral complexity can be introduced by the combined action of bifurcations, lability in key model parameters, and internephron coupling. The complexity includes characteristics identified in the experimental spectrum of Fig. 1*B2*, including strong peaks in the range 25–50 mHz, doublet and triplet peaks, and harmonics. The degree of complexity exemplified in Fig. 9*B2* suggests that the introduction of additional parameter lability could produce spectra as complex as that shown in Fig. 1*B2*.

## DISCUSSION

The aim of this investigation is to attempt to account for the complex spectral signature of TGF-mediated LCO in SHR by means of a minimal model of the TGF system in a single nephron or in a small group of nephrons. We previously developed the mathematical model and have used it to explain the emergence of oscillations in nephron flow in terms of feedback loop gain and delays (26, 27, 37); to characterize the filter properties of the TAL and the spectral properties of the feedback loop (28, 29); and to investigate the potential physiological significance of TGF oscillations in normal rats (30, 35).

In our earlier studies, we investigated TGF system stability by analysis of a characteristic equation derived under the assumption of zero backleak of NaCl in the model TAL (26, 38), or by analysis of a more general characteristic equation, derived for nonzero backleak, which was solved for isolated parameter values (25, 28–30). Analysis of the characteristic equation allowed us to predict the stable dynamic state of the TGF system as a function of its fundamental physical and biological parameters. Our TGF model and its characteristic equation have been a useful guide to TGF system behavior in the parameter regime typical of mature, normotensive rats, which do not exhibit irregular oscillations. For the present investigation, we developed a method to obtain general numerical solutions of the characteristic equation for the case of nonzero backleak; this facilitated a comprehensive analysis that included the effects of both TAL NaCl backleak and TAL transit time on potential dynamic behaviors of the TGF model system.^{2}

### Model Predictions

The most important prediction of the present study is that the TGF system may exhibit multiple stable solutions (SS, *f*_{1}-LCO, and *f*_{2}-LCO) in a parameter regime near that of normotensive rats (see Fig. 3*B2*). Moreover, a modest increase in TGF gain can place the TGF operating point within the region where multistabilty is predicted. In this region, small fluctuations in gain, TAL flow, or system time delays can produce switching between differing stable oscillatory modes. This surprisingly rich dynamic model behavior arises from the inclusion of NaCl backleak in the TAL. Backleak introduces a degree of spatial dependence into our TGF model because of the influence on net transport of the increase in interstitial NaCl concentration between the cortex and the medulla.

By means of numerical simulations suggested by analysis of our model's characteristic equation, we have shown that the representation of three phenomena in our model TGF system is sufficient to produce a degree of TGF system spectral complexity similar to that measured in SHR. In each case, the phenomena are reasonable and/or consistent with experimental observations. The phenomena are the following: *1*) SHR exhibit higher TGF gain and may have longer TGF transit times or shorter TGF response times than normotensive controls (4, 9, 15, 32, 46); *2*) there is a modest degree of vascular coupling among a small number of neighboring nonidentical nephrons (7, 14, 22, 49); and *3*) some system parameters exhibit time variability.

Our analysis predicts that elevated TGF gain in SHR places the TGF operating point in a bistable region where both *f*_{1}-LCO and *f*_{2}-LCO are possible (see Fig. 3, *B1*–*B3*). Furthermore, elongation of the TAL transit time in SHR, as represented in Fig. 4 for 1/η < 1, shifts the operating point to the left into a parameter regime where the regions supporting SS, *f*_{1}-LCO, and both *f*_{1}-LCO and *f*_{2}-LCO are so closely spaced that switching between modes is likely given even modest time variation in bifurcation parameters such as TAL flow or TGF gain. Alternatively, a reduction in feedback time delay across the MD would reduce the bifurcation parameter τ and would thus also shift the TGF operating point to an analogous region, as illustrated in Fig. 3*B2*. More generally, decreases in the quotient τ/η, by means of a simultaneous decrease in τ and increase in η, will tend to shift the operating point toward the bistable region.

These bifurcation phenomena could account for three common features of TGF dynamics in SHR. The first is the presence of multiple fundamental frequencies and harmonic series in SHR spectra (see analysis of Fig. 1*B2* in our introduction; also see Fig. 1, *E* and *F*, in Ref. 52). Our results show that this could be explained by abrupt shifts in parameter values, especially TAL flow rate, that lead to switching between *f*_{1}-LCO and *f*_{2}-LCO, or by internephron coupling. The second is the observation that SHR do not exhibit *f*_{1}-LCO; only steady flows or complex oscillations are seen (16, 52). Our model predicts that many nephrons in SHR may normally operate in the parameter region where *f*_{2}-LCO can emerge directly from steady flows. Finally, our results suggest that increased TAL transit time in SHR could account for the somewhat lowered TGF fundamental frequency in SHR, relative to normotensive rats (16, 52).

### Applicable Findings From the Experimental Record

An important aspect of our hypothesis is that the TGF system in SHR rats frequently operates in a parameter regime where *f*_{1}-LCO and *f*_{2}-LCO are predicted. This requires, compared with control rats, higher TGF gain and differences in TGF system time delays in SHR. Although there is convincing evidence for higher TGF gain in SHR (11, 32, 46), there is less evidence concerning TGF time delays.

Daniels and Arendshorst (9) analyzed the dynamic TGF response to step perturbations in loop of Henle flow. They found significantly shorter (∼33%) TGF response time delays in SHR compared with WKY control rats. Furthermore, although not explicitly analyzed, the initial slope of the TGF response in SHR appears to be markedly steeper than in WKY rats (see Fig. 5 in Ref. 9). However, the measured reduction in TGF response time delay includes *1*) propagation of the flow perturbation through the loop of Henle, *2*) the time delay associated with TGF signal transmission across the MD, and *3*) the AA response time. If the response time delay in SHR were attributable mostly to delays of *types 2* and *3*, then our parameter τ would be much reduced, by perhaps as much as ∼50% or more (a decrease of 50% in τ would diminish its normalized base-case value, 0.22, to 0.11), whereas η, which is proportional to transit time, would be decreased by a smaller percentage of its base-case value. If SHR gain doubled from a plausible (normotensive) value of ∼3.5 (30) to ∼7, the resulting SHR values of τ, η, and γ would be sufficient, in many cases, to move parameters from our estimated normotensive range [around (τ, 1/η, γ) = (0.22, 1.0, 3.5)] into the bistable region or into the *f*_{2}-LCO region, despite increased 1/η; see Fig. *5B*.

An elongation of TAL transit time, which would decrease 1/η, could arise from either increased TAL volume or decreased TAL flow. We could not find direct experimental data regarding TAL volume in SHR. However, kidney weights in SHR, when normalized for body weight, have been found to be similar to (48), or larger than (40), kidney weights in WKY control rats, which suggests similar or larger TAL volumes between strains, if TAL volume scales with kidney weight. The available data concerning nephron flows in SHR are mixed. In young SHR, SNGFR and distal flows were found to be markedly reduced relative to WKY rats (11, 12), whereas the same laboratory found no differences between adult SHR and WKY rats (1, 2). In contrast, Takabatake (46) compared TGF responses in SHR and WKY rats and found uniformly lower mean early proximal flow rates in SHR than in WKY rats. (Because of the experimental design where six groups were compared, only the flows at loop perfusion rates of 20 and 40 nl/min were significant in that study.) Furthermore, indirect evidence from a number of TGF studies suggests lower TAL flows in SHR relative to normotensive rats. In those studies where TGF gain in SHR was found to be increased relative to controls (4, 11, 15, 46), the turning point of the TGF curve (which corresponds to the late proximal tubule perfusion rate at which TGF sensitivity is highest) is significantly lower (∼25–30%) in SHR than in normotensive rats. Separate studies by the same investigators reported either TGF-dependent irregular fluctuations in tubular pressure (16) or efficient TGF-dependent autoregulation in SHR (10, 12), behaviors that require high TGF gain. This suggests that TGF has been reset to a lower range of late proximal flows such that the TGF operating point is near the turning point of the TGF curve.

Published data also contain evidence that supports the hypothesis that the TGF system in SHR rats may exhibit multiple oscillatory modes or mode switching. Sosnovtseva et al. (43) used an elegant double-wavelet analysis to investigate frequency and amplitude modulation of the myogenic autoregulatory mechanism by the TGF system. That study contains results from records of tubular pressure oscillations in SHR that are consistent with a rapid, stable (∼80 s) shift from an *f*_{1}-LCO to an *f*_{2}-LCO after a brief excursion to what may be an excitation of an *f*_{3} oscillatory mode (see Fig. 4*D*, filled circles, in Ref. 43). Further evidence comes from inspection of temporal records of tubular pressure oscillations in SHR rats. Yip et al. (52) presented an SHR spectrum (Fig. 1*F* in their study) that has a number of well-defined peaks, including a major peak at 23 mHz and several major peaks near its first harmonic of 46 mHz, viz., peaks at 39, 43, and 47 mHz. The peaks near 46 mHz are clearly too large (in a spectrum having a linear-scale ordinate) to be harmonics of the 23-mHz oscillation (see, e.g., in their Fig. 1*D*, the harmonic of the 34-mHz LCO, which is hardly discernible at 67 mHz). Close examination of the corresponding time record (in their Fig. 1*B*) reveals the episodic appearance of well-defined LCO within the frequency band predicted for *f*_{2}-LCO. The same behavior can be seen in the data from a Goldblatt hypertensive rat (Fig. 1*C* in Ref. 52) and in the SHR time series reported by Holstein-Rathlou and Leyssac that is included in the present study as Fig. 1*B2*.

Based on the foregoing considerations, we believe that available experimental evidence lends substantial support to our hypothesis that the parameters of SHR frequently reside within a bistable (*f*_{1} or *f*_{2}) region, although we acknowledge that the evidence is less than conclusive. Additional experimental data are needed, especially concerning possible alterations in TGF system time delays in SHR relative to other strains. Also needed are careful analyses of experimental data, by means of modern signal processing methods, to identify TGF mode switching in SHR.

### Complexity in SHR

Our model's prediction of multiple TGF oscillatory modes in SHR provides another source of complexity in TGF dynamics. Our results suggest that coupling between nephrons exhibiting different TGF oscillatory modes, as exhibited in the model simulations in Fig. 8, *A1*–*B3*, can result in two (or more) dominant peaks in the power spectra. Such differences in oscillatory mode can arise in SHR from small differences in nephron time delays, which can arise from small differences in nephron volumes and flows. Furthermore, our simulations indicate that coupling between nephrons with differing oscillatory frequencies may result in mixed-mode oscillations, in which both the *f*_{1}-LCO and *f*_{2}-LCO are significantly excited (Fig. 8, *A1*–*A3*). Because the *f*_{2}-LCO is not a simple integer multiple of the fundamental frequency, the waveforms are very irregular and are similar in many respects to in vivo recordings. Additional complexity is introduced when three nephrons with somewhat differing parameters are coupled (Fig. 8, *B1*–*B3*). Hence, we hypothesize that the irregular oscillations observed in hypertensive rats (16, 52) could arise from internephron coupling between two or more nephrons, with at least one of them having parameters in an *f*_{2}-LCO regime.

The representation of nephron coupling in our model is based on physiological evidence which indicates that individual nephrons are not functionally autonomous but can be influenced by neighboring nephrons. It has been shown in rats, rabbits, mice, and humans that the majority of nephrons are organized in pairs or triplets, and that their AAs arise from the same cortical radial artery (6). Such coupled nephrons can exhibit synchronous oscillations (14, 21), and internephron coupling appears to be stronger in SHR (7, 49). A recent modeling study predicted that coupling can significantly reduce the gain threshold for LCO emergence (38).

Our findings that internephron coupling can introduce substantial spectral complexity builds on and extends previous simulation studies by Holstein-Rathlou and colleagues (3, 20, 42, 44, 45). They used a series of multinephron models in simulations designed to investigate the possible presence of deterministic chaos in this system as the phenomenon responsible for the emergence of irregular, chaos-like oscillations in SHR. However, their predictions differ in two important ways from in vivo observations in SHR. First, the emergence of irregular oscillations in their simulations follows the period-doubling route to deterministic chaos. Such systems exhibit LCO, and they undergo a series of period-doubling bifurcations as gain increases, but these behaviors do not appear to be typical of SHR (16). Second, the irregular oscillations only emerge in parameter regimes beyond those measured in SHR. The differences between the findings of Holstein-Rathlou and colleagues (3, 20, 42, 44, 45) and our simulation results may be explained by differences in the TAL models used.In their multinephron models, they employed a representation of the TAL that lacked spatial dependence. Our results suggest that this simplification eliminates much of the rich dynamics exhibited in our simulations when using parameter values closer to measured values.

Nevertheless, it is clear that our TGF models of a single nephron or small numbers of coupled nephrons, when employed with constant parameters, can provide no explanation for those distinct frequencies in the experimental power spectra that cannot be attributed to *f*_{1}-LCO or *f*_{2}-LCO, their harmonics, or their respective doublets. However, our simulations show that the introduction of time-varying parameters can introduce much additional complexity in the TGF power spectrum, as illustrated in Figs. 5, 6, 7, and 9. The most dramatic effect was seen with variation in η, which represents changes in TAL flow or volume. A clear example of such a change in tubular flow in SHR can be seen in the elegant study by Walstead and Yip (50), who measured fluctuating proximal tubular flow in SHR with a nonobstructive video technique. In their data records (see, e.g., Fig. 1 in Ref. 50), abrupt shifts in proximal flow are present, as are multiple strong peaks in the associated power spectrum.

Finally, our results complement those of a recent modeling study by Ditlevsen et al. (13), which predicted that rapid and high-amplitude stochastic fluctuations in TGF sensitivity can produce irregular TGF oscillations. Although we did not explicitly include stochastic perturbations in the present study, random fluctuations in flows and parameter values are undoubtedly present in vivo. In our model of TGF in SHR, even small fluctuations in gain and time delays would be sufficient to result in switching of oscillatory modes and highly irregular oscillations. The difference in the behaviors of these two TGF models is likely related to the empirical third-order linear system employed by Ditlevsen et al. to represent the dynamic contribution of the TAL. That third-order system, which has been used in a series of TGF studies (3, 20, 42, 44, 45), may not display the complex dynamics exhibited by our simple, but much more realistic, spatially distributed model of the TAL.

In traditional power spectra that are computed from long-time data records, the spectral features that are produced by temporal variation in TGF system parameters will accumulate by a process of superposition and produce a spectrum that is, in many respects, similar to a time-exposure photograph. Hence, slow variation in some parameters will tend to broaden spectral peaks, whereas an abrupt shift in a parameter will result in the superposition of peaks from the preshift spectrum with peaks that represent the state of the system after the shift. Examples of both slow parameter variation and abrupt parameter shifts appear to be present in SHR spectra, which exhibit both broadened peaks and closely spaced peaks of similar amplitude in spectra having a linear-scale ordinate. However, the low-frequency peaks in the SHR spectra (below ∼15 or ∼20 mHz) are difficult to reconcile with any reasonable shift of system parameters. These slow oscillations in SHR most likely arise from extrarenal factors (i.e., renal sympathetic activity) that influence the TGF system. In this regard, several recent studies (8, 54) have demonstrated that renal hemodynamics in normotensive rats are sufficiently stationary to permit analysis by time-invariant spectral methods. In contrast, SHR exhibit oscillations that are sufficiently nonstationary to require the application of time-varying spectral methods, several of which have been developed in recent years (33, 43, 51). These powerful methods show the evolution of the frequency spectrum over time, thereby making it possible to separate the intrinsic dynamics of the TGF system from the dynamics of other systems that influence TGF.

Taken together, our simulations provide the basis for an alternative explanation of the irregular fluctuations in the TGF system in SHR. Rather than emerging from a nonlinear dynamical system capable of exhibiting deterministic chaos (52), the complex fluctuations may instead arise mostly from an intrinsic multistability of the TGF system. This multistability, in conjunction with significant coupling of a small number of neighboring nephrons and temporal variations in key parameters, can result in the complex power spectra similar to those found in SHR. Understanding hemodynamic disregulation in SHR will require experimental verification of the key predictions of our TGF model, reanalysis of data records with time-varying spectral methods, and elucidation of the identity and dynamics of the processes that impinge upon the TGF system.

## APPENDIX

### Normalization of Model Equations

For generality in the derivation of the characteristic equations (*Eqs. 9* and *10*), and for simplicity in our numerical simulations, we used forms of *Eqs. 1*–*5* that had been normalized to render them nondimensional. The normalization was accomplished by dividing the dimensional variables and parameters by dimensional reference values. Thus normalized variables and parameters, indicated by a tildesymbol, *x̃ = x/L*, *t̃* = *t*/*t*_{o}, s̃ = s*/t*_{o}, C̃_{i}(*x̃,t̃*) = C_{i}(*x,t*)*/*C_{o}, C̃(*t̃*) = C(*t*)*/*C_{o}, _{i}(*t̃*) = (*t*)*/*F_{op}, F̃_{i}(*t̃*) = F_{i}(*t*)*/*F_{op}, C̃_{e}(*x̃*) = C_{e}(*x*)*/*C_{o}, *Ã*_{o} = *A*_{o}*/A*_{o}, c̃= c*/*c(where c= 2π*r*), *Ṽ*_{max} = *Ṽ*_{max}/(C_{o} F_{op}*/*(c_{Ao}*L*)), *K̃*_{M} *= K*_{M}*/*C_{o}, *P̃ = P/*(F_{op}/(c*L*)), *K*_{1} =ΔQ/(2Q), *K*_{2i}*= k*_{i}C_{o}/2, δ̃ = δ*/t*_{o}, τ̃_{i} *= τ*_{i}/*t*_{o}, τ̃_{pi} = τ_{pi}/*t*_{o}, ρ̃_{i}(*t̃*) = ρ_{i}(*t*)/F_{op}, and ψ̃_{δ̃}(s̃) = ψ_{δ}(s)/(1/*t*_{o}). When these normalized variables are substituted into the model equations, and the tilde symbols are dropped, the normalized forms of *Eqs. 1*–*3* are given by (A1) (A2) (A3) The forms of *Eqs. 4* and *5* are unchanged by normalization.

### Numerical Methods and Procedures

The solutions to the model equations (*Eqs. 1*–*5*, although normalized as described above) were closely approximated by numerical solutions. The numerical solutions were computed, as previously described (24, 29), by solving finite-difference equations that correspond to the model equations: a spatially second-order ENO method was used in conjunction with Heun's method to obtain a method that was second order in both space and time. The spatial discretization used 640 subintervals, which yielded a space step of Δ*x* = (0.5 cm)/640; the time step was Δ*t*= 1/320 s.

Numerical solutions of the model equations were used to confirm the loci of transitions from SS to LCO that are predicted by the characteristic equation. The model equations were first solved for stable solutions at two points on differing sides of a predicted bifurcation locus; if this test indicated a bifurcation, successive tests were conducted on each side of the locus to localize the bifurcation site. This procedure was repeated at several such points along each predicted locus curve. A bifurcation locus curve that was not precisely predicted by the characteristic equation (e.g., the “Empirical Boundary” on Fig. 3*B2*) was first located approximately by random sampling within the region of the suspected bifurcation. Then the curve was localized, to numerical accuracy, by testing densely spaced points along lines nearly perpendicular to the bifurcation curve. By means of these procedures, the diagrams given by Figs. 3, A2 and *B2*, 4*A2*, and 5, *A* and *B*, were constructed.

## GRANTS

This research was principally supported by National Institute of Diabetes and Digestive and Kidney Diseases Grant DK-42091 (to H. E. Layton). Additional support was provided by National Science Foundation Grant DMS-0340654 (to A. T. Layton).

## Acknowledgments

We are grateful to Michael C. Mackey and to Michael C. Reed for helpful suggestions.

Portions of this work were presented at Experimental Biology 2000 (*FASEB J* 14: A135, 2000; abstract no. 119.9) and ASN Renal Week 2002 (*J Am Soc Neph* 13: 333A, 2002; abstract no. SA-P0343).

## Footnotes

↵1 The acronym “SS” will be taken to mean a steady state that is both time independent and stable. By “time independent,” we mean that no model variables vary as a function of time; by “stable,” we mean that, after any transient perturbation of the steady state, the model solutions, as a function of time, will more and more closely approximate the SS model solution.

↵2 As early as 1991, we had observed that a second oscillatory mode could be transiently excited (for ∼300 s) in our model (see Fig. 13 in Ref. 26); however, in the absence of a general analysis of the characteristic equation, we did not know how to interpret that finding. We now recognize that we had excited the

*f*_{2}mode near the empirical boundary but within the region where the*f*_{1}-LCO is the only stable mode (see Fig. 3*B2*in the present study). Thus this prior result demonstrated that an unstable mode can be excited for a significant interval before returning to a stable mode.The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 the American Physiological Society