Multistability in tubuloglomerular feedback and spectral complexity in spontaneously hypertensive rats

【摘要】  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).

Fig. 1. A1 : experimental pressure record from proximal tubule of Wistar-Kyoto (WKY) rat showing a spontaneous limit-cycle oscillation (LCO). A2 : power spectrum for LCO in A1, in arbitrary units, showing relative strength of frequency components. B1 : corresponding pressure record from a spontaneously hypertensive rat (SHR) showing highly irregular oscillation. B2 : power spectrum for oscillation in B1. Labels of peaks in A2 and B2, added by us, give peak frequencies (in mHz). Adapted from Ref. 16 with permission from Blackwell Publishing.

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.

Fig. 2. Schematic representation of model configuration. Model represents 3 essential elements of the turbuloglomerular feedback (TGF) pathway: thick ascending limb (TAL; cylinder), delay at the macula densa (MD; box on the right ), and TGF response function (box on the left ). Symbols are identified in the GLOSSARY. Perturbations and coupling enter as modifications to TAL flow rate: when coupling is introduced, F i ( t ) is replaced by i ( t ); transient perturbations are administered by adding ( t ) to F i ( t ) or to i ( t ).

Fig. 3. Diagrams showing predictions of the characteristic equation, as a function of delay and gain magnitude, for the cases of zero ( A1 ) and nonzero ( B1 ) Cl - backleak into TAL; and bifurcation diagrams indicating observed behavior of model solutions for zero ( A2 ) and nonzero ( B2 ) backleak. SS, stable, time-independent, steady state; f 1, region having f 1 -LCO; f 2, region having f 2 -LCO; f 1 or f 2, region having either an f 1 -LCO or an f 2 -LCO; gray disk, likely parameter region of normotensive rats. Points in A2 and B3 labeled "A3" and "B3," respectively, refer to panels A3 and B3, which show corresponding LCO in model single-nephrom glomerular filtration rate (SNGFR) and TAL tubular fluid Cl - concentration at the MD. (Gray lines indicate SS solutions.) For the same values of gain and delay, the zero-backleak case has as its only solution an f 1 -LCO, whereas the nonzero-backleak case can also have an f 2 -LCO at nearly twice the frequency of the f 1 -LCO.

Fig. 4. A1 : diagram showing predictions of the characteristic equation for nonzero backleak in the (1/ )- plane. Inset : mutually perpendicular positive axes for the parameters, 1/, and; the origin is designated as "0." A2 : corresponding bifurcation diagram indicating observed behavior of stable solutions. B1 : SNGFR response to transient perturbation. B2 : power spectrum for a SS, both corresponding to point B in A2. C1 and C2 : f 1 -LCO and corresponding power spectrum for point C in A2, respectively. D1 and D2 : f 2 -LCO and corresponding power spectrum for point D in A2, respectively.

Fig. 5. Bifurcation diagrams in the -(1/ ) plane for = 5 ( A ) and = 7 ( B ). Arrows indicate paths of parameters that were varied to produce results shown in Figs. 6 and 7. Dot marks ( o, 1/ o ).

Fig. 6. Power spectra for SNGFR oscillations with variation in and for = 5, as indicated in Fig. 5 A. A1 and A2 : spectra with linear and logarithmic ordinate, respectively, for varying (solid curves) and (dashed curves) at constant speeds along paths indicated in Fig. 5 A. Gray curves, power spectra for stationary parameters corresponding to solid dot in Fig. 5 A. A3 and A4 : spectra with linear and logarithmic ordinate, respectively, for abruptly varying (solid curve) along path indicated by arrow in Fig. 5 A, between intervals of resting at the ends of the arrow. Gray curves, power spectra for stationary parameters corresponding to solid dot in Fig. 5 A.

Fig. 7. Power spectra for SNGFR oscillations with dynamic variation in and for = 7, as indicated in Fig. 5 B. B1 and B2 : spectra with linear and logarithmic ordinate, respectively, for varying (solid curves) and (dashed curves) at constant speeds along paths indicated in Fig. 5 B. Gray curves, power spectra for stationary parameters corresponding to solid dot in Fig. 5 A. B3 and B4 : spectra with linear and logarithmic ordinate, respectively, for abruptly varying (dashed curve) along path indicated by arrow in Fig. 5 B, between intervals of resting at the ends of the arrow. Gray curves, power spectra for stationary parameters corresponding to solid dot in Fig. 5 B.

Fig. 8. A1 : SNGFR oscillations in one of two coupled nephrons and corresponding power spectra with linear ( A2 ) and logarithmic ( A3 ) ordinates. B1 : SNGFR oscillations in one of three coupled nephrons and corresponding power spectra with linear ( B2 ) and logarithmic ( B3 ) ordinates.

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 2 i S '( 1 )

c Critical feedback gain magnitude (dimensionless)

Distributed delay interval at the JGA (s)

i TAL transit time scaling (dimensionless)

i p i + /2, effective delay interval at the JGA (s)

p i 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

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 MD i 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

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

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

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 ).

Table 1. Base-case parameter values

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

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 x 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 x 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 x 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 x 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 x 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.

Table 2. Parameters for coupled nephrons

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 8 B1-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.

Fig. 9. A1 : SNGFR oscillations in 1 of 3 coupled nephrons after abrupt 10% increase in 1/ i for each nephron, relative to the case of Fig. 8 B; corresponding power spectra with linear ( A2 ) and logarithmic ( A3 ) ordinates are based on the entire time series, both before and after increase in 1/ i. In A2 and A3, the peaks labeled with larger, Roman font correspond to newly introduced peaks, relative to Fig. 8, B2 and B3. B1 : SNGFR oscillations for the case of A1 but having, in addition, a sustained sinusoidal oscillation in 1/ i having a frequency of 10 mHz and an amplitude of 1/(10 i ), shown as corresponding power spectra with linear ( B2 ) and logarithmic ( B3 ) ordinates. In B2 and B3, the peaks labeled with larger, Roman font correspond to newly introduced peaks, relative to A2 and A3.

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 transmis

【参考文献】
  Arendshorst WJ and Beierwaltes WH. Renal and nephron hemodynamics in spontaneously hypertensive rats. Am J Physiol Renal Fluid Electrolyte Physiol 236: F246-F251, 1979.

Arendshorst WJ and Beierwaltes WH. Renal tubular reabsorption in spontaneously hypertensive rats. Am J Physiol Renal Fluid Electrolyte Physiol 237: F38-F47, 1979.

Barfred M, Mosekilde E, and Holstein-Rathlou N-H. Bifurcation analysis of nephron pressure and flow regulation. Chaos 6: 280-287, 1996.

Brännström K, Morsing P, and Arendshorst WJ. Exaggerated tubuloglomerular feedback activity in genetic hypertension is mediated by ANG II and AT 1 receptors. Am J Physiol Renal Fluid Electrolyte Physiol 270: F749-F755, 1996.

Casellas C and Moore LC. Autoregulation and tubuloglomerular feedback in juxtamedullary glomerular arterioles. Am J Physiol Renal Fluid Electrolyte Physiol 258: F660-F669, 1990.

Casellas D, Dupont M, Bouriquet N, Moore LC, Artuso A, and Mimran A. Anatomic pairing of afferent arterioles and renin cell distribution in rat kidneys. Am J Physiol Renal Fluid Electrolyte Physiol 267: F931-F936, 1994.

Chen YM, Yip KP, Marsh DJ, and Holstein-Rathlou N-H. Magnitude of TGF-initiated nephron-nephron interaction is increased in SHR. Am J Physiol Renal Fluid Electrolyte Physiol 269: F198-F204, 1995.

Chon KH, Raghavan R, Chen YM, Marsh DJ, and Yip KP. Interactions of TGF-dependent and myogenic oscillations in tubular pressure. Am J Physiol Renal Physiol 288: F298-F307, 2005.

Daniels FH and Arendshorst WJ. Tubuloglomerular feedback kinetics in spontaneously hypertensive and Wistar-Kyoto rats. Am J Physiol Renal Fluid Electrolyte Physiol 259: F529-F534, 1990.

Daniels FH, Arendshorst WJ, and Roberds RG. Tubuloglomerular feedback and autoregulation in spontaneously hypertensive rats. Am J Physiol Renal Fluid Electrolyte Physiol 258: F1479-F1489, 1990.

Dilley JR and Arendshorst WJ. Enhanced tubuloglomerular feedback activity in rats developing spontaneous hypertension. Am J Physiol Renal Fluid Electrolyte Physiol 247: F672-F679, 1984.

Dilley JR, Stier CT Jr, and Arendshorst WJ. Abnormalities in glomerular function in rats developing spontaneous hypertension. Am J Physiol Renal Fluid Electrolyte Physiol 246: F12-F20, 1984.

Ditlevsen S, Yip KP, and Holstein-Rathlou N-H. Parameter estimation in a stochastic model of the tubuloglomerular feedback mechanism in a rat nephron. Math Biosci 194: 49-69, 2005.

Holstein-Rathlou N-H. Synchronization of proximal intratubular pressure oscillations: evidence for interaction between nephrons. Pflügers Arch 408: 438-443, 1987.

Holstein-Rathlou N-H. A closed-loop analysis of the tubuloglomerular feedback mechanism. Am J Physiol Renal Fluid Electrolyte Physiol 261: F880-F889, 1991.

Holstein-Rathlou N-H and Leyssac PP. TGF-mediated oscillations in proximal intratubular pressure: differences between spontaneously hypertensive rats and Wistar-Kyoto rats. Acta Physiol Scand 126: 333-339, 1986.

Holstein-Rathlou N-H and Marsh DJ. Oscillations of tubular pressure, flow, and distal chloride concentration in rats. Am J Physiol Renal Fluid Electrolyte Physiol 256: F1007-F1014, 1989.

Holstein-Rathlou N-H and Marsh DJ. A dynamic model of the tubuloglomerular feedback mechanism. Am J Physiol Renal Fluid Electrolyte Physiol 258: F1448-F1459, 1990.

Holstein-Rathlou N-H and Marsh DJ. Renal blood flow regulation and arterial pressure fluctuations: a case study in nonlinear dynamics. Physiol Rev 74: 637-681, 1994.

Jensen KS, Holstein-Rathlou N-H, Leyssac PP, Mosekilde E, and Rasmussen DR. Chaos in a system of interacting nephrons. In: Life Sciences: Chaos in Biological Systems, edited by Degn H, Holden AV, and Olsen LF. New York, Plenum, 1987.

Just A, Wittman U, Ehmke H, and Kirchheim HR. Autoregulation of renal blood flow in the conscious dog and the contribution of the tubuloglomerular feedback. J Physiol 506: 275-290, 1998.

Källskog Ö and Marsh DJ. TGF-initiated vascular interactions between adjacent nephrons in the rat kidney. Am J Physiol Renal Fluid Electrolyte Physiol 259: F60-F64, 1990.

Kesseler KJ. Analysis of Feedback-Mediated Oscillations in Two Coupled Nephrons (PhD thesis). Durham, NC: Duke University, 2004.

Layton HE and Pitman EB. A dynamic numerical method for models of renal tubules. Bull Math Bio 56: 547-565, 1994.

Layton HE, Pitman EB, and Knepper MA. A dynamic numerical method for models of the urine concentrating mechanism. SIAM J Appl Math 55: 1390-1418, 1995.

Layton HE, Pitman EB, and Moore LC. Bifurcation analysis of TGF-mediated oscillations in SNGFR. Am J Physiol Renal Fluid Electrolyte Physiol 261: F904-F919, 1991.

Layton HE, Pitman EB, and Moore LC. Instantaneous and steady-state gains in the tubuloglomerular feedback system. Am J Physiol Renal Fluid Electrolyte Physiol 268: F163-F174, 1995.

Layton HE, Pitman EB, and Moore LC. Nonlinear filter properties of the thick ascending limb. Am J Physiol Renal Physiol 273: F625-F634, 1997.

Layton HE, Pitman EB, and Moore LC. Spectral properties of the tubuloglomerular feedback system. Am J Physiol Renal Physiol 273: F635-F649, 1997.

Layton HE, Pitman EB, and Moore LC. Limit-cycle oscillations and tubuloglomerular feedback regulation of distal sodium delivery. Am J Physiol Renal Physiol 278: F287-F301, 2000.

Leyssac PP and Baumbach L. An oscillating intratubular pressure response to alterations in Henle loop flow in the rat kidney. Acta Physiol Scand 117: 415-419, 1983.

Leyssac PP and Holstein-Rathlou N-H. Tubulo-glomerular feedback response: enhancement in adult spontaneously hypertensive rats and effects of anaesthetics. Pflügers Arch 413: 267-272, 1989.

Marsh DJ, Sosnovtseva OV, Pavlov AN, Yip KP, and Holstein-Rathlou N-H. Frequency encoding in renal blood flow regulation. Am J Physiol Regul Integr Comp Physiol 288: R1160-R1167, 2005.

Oldson DR. Flow Perturbations in a Mathematical Model of the Tubuloglomerular Feedback System. (PhD thesis). Durham, NC: Duke University, 2003.

Oldson DR, Layton HE, and Moore LC. Effect of sustained flow perturbations on stability and compensation of tubuloglomerular feedback. Am J Physiol Renal Physiol 285: F972-F989, 2003.

Peitgen HO, Jürgens H, and Saupe D. Chaos and Fractals: New Frontiers of Science. New York: Springer-Verlag, 1992.

Pitman EB, Layton HE, and Moore LC. Dynamic flow in the nephron: filtered delay in the TGF pathway. Contemp Math 114: 317-336, 1993.

Pitman EB, Zaritski RM, Kesseler KJ, Moore LC, and Layton HE. Feedback-mediated dynamics in two coupled nephrons. Bull Math Biol 66: 1463-1492, 2004.

Pitman EB, Zaritski RM, Moore LC, and Layton HE. A reduced model for nephron flow dyanmics mediated by tubuloglomerular feedback. In: Membrane Transport and Renal Physiology, The IMA Volumes in Mathematics and Its Applications, edited by Layton HE and Weinstein AM. New York: Springer, 2002, vol. 129, p. 345-364.

Saleh FH and Jurjus AR. A comparative study of morphological changes in spontaneously hypertensive rats and normotensive Wistar Kyoto rats treated with an angiotensin-converting enzyme inhibitor or a calcium-channel blocker. J Pathol 193: 415-420, 2001. <a href="/cgi/external_ref?access_num=10.1002/1096-9896(2000)9999:9999

Schnermann J and Briggs JP. Function of the juxtaglomerular apparatus: control of glomerular hemodynamics and renin secretion. In: The Kidney: Physiology and Pathophysiology (3rd ed.), edited by Seldin DW and Giebisch G. Philadelphia, PA: Lippincott Williams & Wilkins, 2000, p. 945-980.

Sosnovtseva OV, Pavlov AN, Mosekilde E, and Holstein-Rathlou N-H. Bimodal oscillations in nephron autoregulation. Phys Rev E Stat Nonlin Soft Matter Phys 66: 061909-1-061909-7, 2002.

Sosnovtseva OV, Pavlov AN, Mosekilde E, Holstein-Rathlou N-H, and Marsh DJ. Double-wavelet approach to study frequency and amplitude modulation in renal autoregulation. Phys Rev E Stat Nonlin Soft Matter Phys 70: 031915-1-031915-8, 2004.

Sosnovtseva OV, Postnov DE, Mosekilde E, and Holstein-Rathlou N-H. Synchronization of tubular pressure oscillations in interacting nephrons. Chaos, Solitons & Fractals 15: 343-369, 2003.

Sosnovtseva OV, Postnov DE, Nekrasov AM, Mosekilde E, and Holstein-Rathlou N-H. Phase multistability of self-modulated oscillations. Phys Rev E Stat Nonlin Soft Matter Phys 66: 036224-1-036224-9, 2002.

Takabatake T, Ushiogi Y, Ohta K, and Hattori N. Attenuation of enhanced tubuloglomerular feedback activity in SHR by renal denervation. Am J Physiol Renal Fluid Electrolyte Physiol 258: F980-F985, 1990.

Thomson SC and Blantz RC. Homeostatic efficiency of tubuloglomerular feedback in hydropenia, euvolemia, and acute volume expansion. Am J Physiol Renal Fluid Electrolyte Physiol 264: F930-F936, 1993.

Uyehara CFT and Gellai M. Impairment of renal function preceded establishment of hypertension in spontaneously hypertensive rats. Am J Physiol Regul Integr Comp Physiol 265: R943-R950, 1993.

Wagner AJ, Holstein-Rathou NH, and Marsh DJ. Internephron coupling by conducted vasomotor responses in normotensive and spontaneously hypertensive rats. Am J Physiol Renal Physiol 272: F372-F379, 1997.

Walstead C and Yip K-P. Acute arterial hypertension inhibits proximal tubular reabsorption in normotensive rat but not in SHR. Am J Physiol Regul Integr Comp Physiol 286: R726-R733, 2004.

Wang H, Siu K, Ju K, Moore LC, and Chon KH. Identification of transient renal autoregulatory mechanisms using time-frequency spectral techniques. IEEE Trans Biomed Eng 52: 1033-1039, 2005.

Yip KP, Holstein-Rathlou N-H, and Marsh DJ. Chaos in blood flow control in genetic and renovascular hypertensive rats. Am J Physiol Renal Fluid Electrolyte Physiol 261: F400-F408, 1991.

Yip KP, Holstein-Rathlou N-H, and Marsh DJ. Dynamics of TGF-initiated nephron-nephron interactions in normotensive rats and SHR. Am J Physiol Renal Fluid Electrolyte Physiol 262: F980-F988, 1992.

Zou R, Cupples WA, Yip KP, Holstein-Rathou NH, and Chon KH. Time-varying properties of renal autorgulatory mechanims. IEEE Trans Biomed Eng 49: 1112-1120, 2002.

作者单位：1 Department of Mathematics, Duke University, Durham, North Carolina; and 2 Department of Physiology and Biophysics, State University of New York, Stony Brook, New York