MODELING IN PHYSIOLOGY
Spectral properties of the tubuloglomerular feedback system

¹ Department of Mathematics, Duke University, Durham, North Carolina 27708-0320; ² Department of Mathematics, State University of New York, Buffalo 14214-3093; and ³ Department of Physiology and Biophysics, State University of New York, Stony Brook, New York 11794-8661

	ABSTRACT

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

A simple mathematical model was used to investigate the spectral properties of the tubuloglomerular feedback (TGF) system. A perturbation, consisting of small-amplitude broad-band forcing, was applied to simulated thick ascending limb (TAL) flow, and the resulting spectral response of the TGF pathway was assessed by computing a power spectrum from resulting TGF-regulated TAL flow. Power spectra were computed for both open- and closed-feedback-loop cases. Open-feedback-loop power spectra are consistent with a mathematical analysis that predicts a nodal pattern in TAL frequency response, with nodes corresponding to frequencies where oscillatory flow has a TAL transit time that equals the steady-state fluid transit time. Closed-feedback-loop spectra are dominated by the open-loop spectral response, provided that , the magnitude of feedback gain, is less than the critical value _c required for emergence of a sustained TGF-mediated oscillation. For  exceeding _c, closed-loop spectra have peaks corresponding to the fundamental frequency of the TGF-mediated oscillation and its harmonics. The harmonics, expressed in a nonsinusoidal waveform for tubular flow, are introduced by nonlinear elements of the TGF pathway, notably TAL transit time and the TGF response curve. The effect of transit time on the flow waveform leads to crests that are broader than troughs and to an asymmetry in the magnitudes of increasing and decreasing slopes. For feedback gain magnitude that is sufficiently large, the TGF response curve tends to give a square waveshape to the waveform. Published waveforms and power spectra of in vivo TGF oscillations have features consistent with the predictions of this analysis.

kidney; renal hemodynamics; nonlinear dynamics; mathematical model

	INTRODUCTION

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

THE POWER SPECTRUM, which quantifies the relative prevalence of frequency components in a signal, has emerged as an important tool in the analysis of experimental time series derived from renal hemodynamic variables. Spectral analysis helped establish that single-nephron oscillations in intratubular pressure, flow, and distal tubule chloride concentration, with frequency of 20-50 mHz, arise from an intrinsic instability in the tubuloglomerular feedback (TGF) loop (7); spectral analysis showed that significant spectral power was distributed over a range of frequencies in tubular flow that appears to exhibit deterministic chaos (22); and spectral analysis has been used to distinguish oscillations arising from TGF from those of intrinsic myogenic origin (3, 11, 23).

In this study we used a simple mathematical model to investigate the spectral properties of the TGF system. Previously, we have used the same model framework to elucidate the emergence of TGF-mediated oscillations (13, 17), to distinguish steady-state gain from instantaneous gain (14), and to characterize the nonlinear filter properties of the thick ascending limb (TAL) (15).

We first summarize the model equations and associated physiological parameters, and we review the TGF-mediated oscillation that may be exhibited by the model when , the magnitude of gain of the TGF pathway, exceeds a critical value _c. Then, using an open-feedback-loop configuration, we compute power spectra that characterize the low-pass filter of the model TAL and the transmission of the TGF signal from the macula densa (MD) to the afferent arteriole (AA). These results confirm the nodal structure identified in the companion study (15) and provide additional insight into the spectral properties of the TGF system components when operating in the absence of feedback.

The key new results of this study, however, are obtained from the closed-feedback-loop configuration. Power spectra, computed for selected values of feedback-loop gain magnitude , characterize the evolution of the TGF system as  increases through the regime that will not support sustained TGF-mediated oscillations,  < _c, and into the regime that does,  > _c. For  < _c, the power spectra are dominated by the open-loop spectral characteristics, with the fundamental resonant frequency of the TGF oscillation superimposed on the open-loop frequency response. For  > _c, spectral power is increased, and power spectra are composed of peaks corresponding to the fundamental resonant frequency of the TGF-mediated oscillation and its harmonics, which are explained via the principle of Fourier decomposition. The harmonics, which become more pronounced as gain magnitude increases, arise from distortions in TGF oscillations that are introduced by nonlinear properties of the TGF pathway, notably the integrative effect of TAL transit time and the constraints of the TGF response function. The distortions in the waveform of tubular flow arising from nonlinear transit time consist of crests that are broader than troughs and a slope asymmetry, in which the absolute magnitude of the ascending slope is less than the absolute magnitude of the descending slope. For TGF gain magnitude that is sufficiently large, the bounds of the TGF response function tend to give a square waveshape to the waveform. Finally, we observe that published waveforms and power spectra from in vivo measurements of TGF-mediated oscillations frequently exhibit features consistent with the predictions of this theoretical analysis.

Glossary Parameters

Co	Chloride concentration at TAL entrance (mM)
Cop	Steady-state chloride concentration at MD (mM)
k	Sensitivity of TGF response (mM1)
Km	Michaelis constant (mM)
L	Length of TAL (cm)
P	TAL chloride permeability (cm/s)
Qop	Steady-state SNGFR (nl/min)
Q	TGF-mediated range of SNGFR (nl/min)
r	Luminal radius of TAL (µm)
to	Steady-state TAL transit time (s)
Vmax	Maximum transport rate of chloride from TAL (nmol · cm2 · s1)
	Fraction of SNGFR reaching TAL
	Distributed delay interval at the JGA (s)
p	Pure delay interval at the JGA (s)

Independent variables

f	Frequency of flow oscillation (mHz)
t	Time (s)
x	Axial position along TAL (cm)

Specified functions

Ce(x)	Extratubular chloride concentration (mM)
(t)	Kernel function for distributed delay (dimensionless)

Dependent variables

C(x, t)	TAL chloride concentration (mM)
CMD(t)	Effective MD chloride concentration (mM)
F(CMD(t))	TAL fluid flow rate (nl/min)
Pf	Power spectral density
S(x)	Steady-state TAL chloride concentration (mM)
T(x, t)	Fluid transit time from TAL entrance (s)

	MATHEMATICAL MODEL

Model equations. The mathematical model for the TGF loop (14, 17) is given by the following system of coupled equations

(1)

(2)

(3)

Each equation is in nondimensional form, i.e., all variables and parameters have been normalized so that each is a dimensionless quantity (see APPENDIX A). Figure 1 provides a schematic representation of the model components. The space variable x is oriented so that it extends from the entrance of the TAL (x = 0), through the outer medulla, and into the cortex to the site of the MD (x = 1).

Equation 1 is a partial differential equation for the chloride ion concentration C in the intratubular fluid of the TAL of a short-looped nephron. At time t = 0, initial concentrations C(x, 0) (for x  [0, 1]) and C(1, t) (for t  (, 0)) must be specified. We impose the boundary condition C(0, t) = 1, which means that fluid entering the TAL has constant chloride concentration. The rate of change of that concentration at x  (0, 1] depends on processes represented by the three right-hand terms in Eq. 1. The first term is axial convective chloride transport at the intratubular flow speed F. The second is the transtubular efflux of chloride driven by active metabolic pumps situated in the tubular walls; that efflux is approximated by Michaelis-Menten kinetics, with maximum transport rate V_max and Michaelis constant K_m. The third term is transtubular chloride backleak, which depends on a specified fixed extratubular chloride concentration profile C_e(x) and on chloride permeability P.

Equation 2 describes fluid speed in the TAL as a function of the effective luminal chloride concentration C_MD at the MD (see below). This feedback relation is an empirical equation well established by steady-state experiments (20). The constant C_op is the steady-state chloride concentration obtained at the MD when F  1. The positive constants K₁ and K₂ describe, respectively, the range of the feedback response and its sensitivity to deviations from the steady state.

Equation 3 represents time delays in the feedback pathway between the luminal fluid chloride concentration at the MD, C(1, t), and an effective MD concentration C_MD(t), which is used to calculate the flow response that is modulated by smooth muscle of the AA. In quasi-steady state, Eq. 2 provides a good description of the TGF feedback response. However, dynamic experiments (1) show that a change in MD concentration does not significantly affect AA muscle tension until after a pure delay time _p, and then the effect is distributed in time so that a full response requires additional time, with, say, greatest weight in the time interval [t  _p  , t  _p], where  is a second delay parameter. To simulate the pure delay followed by a distributed delay, we introduced the convolution integral in Eq. 3 to describe the effective signal received by the AA at time t (17). We require that the kernel satisfy (t  s  /2) ds = 1, so that a constant value for C(1, t  _p) will be unchanged by the distributed delay of Eq. 3. For this study we assume that is given by

(4)

a function that makes use of one oscillation of a cosine curve, centered at the origin and scaled so that the function is nonnegative and continuously differentiable for all u. With this function, a step change in C results in a sigmoidal increase in C_MD over a nondimensional time interval of  (cf. Eq. 3).

A steady-state solution to Eqs. 1-4 may be obtained by setting F = 1 for 1 unit of normalized time (the transit time of the TAL at flow speed 1), starting at t = 0, to give the steady-state operating concentration C_op = C(1, 1) at the MD. If one specifies that C(1, t) = C_op for t  (, 1), then the input flow to the TAL, F, is fixed at 1 for all previous time. The feedback loop can then be closed at t = 1. If the system remains unperturbed, the system solution will not vary in time. We denote the resulting steady-state TAL concentration profile C(x, 1) by S(x).

Model parameters. A summary of parameters and variables, with their dimensional units as commonly reported, is given in the Glossary. The values of model parameters are given in Table 1; the criteria for their selection and supporting references were given in (13). The extratubular concentration is given in nondimensional form by C_e(x) = C_o(A₁ exp(A₃x) + A₂), where A₁ = (1  C_e(1)/C_o)/(1  e^A₃), A₂ = 1  A₁, and A₃ = 2, and where C_e(1) corresponds to a cortical interstitial concentration of 150 mM. Graphs of C_e and the steady-state luminal profile S were given in figure 1 of Ref. 13. The steady-state operating concentration C_op was calculated numerically using the TAL dimensions and transport parameters, with steady flow F = 1 in Eq. 1.

View this table: [in this window] [in a new window]   Table 1.   Parameters

Rather than a pure delay of 4 s in signal transmission at the juxtaglomerular apparatus (JGA), as we used in Ref. 13, a pure delay _p of 2 s was followed by a transition interval  of 3 s, providing a distributed delay in approximate agreement with experiments (1).

The steady-state TAL transit time t_o is a key parameter that plays a prominent role in this study. This transit time is the interval required for a water molecule to travel up the TAL, from the TAL entrance to the MD, at the steady-state flow rate, assuming plug flow; it is equal to the TAL volume divided by the steady-state TAL flow rate, i.e., t_o = r²L/(Q_op), and it corresponds to one unit of normalized time.

Bifurcation locus. For simplified versions of the model given by Eqs. 1-4, we have previously shown that, for some parameter ranges, the time-independent steady-state solution is unstable, and subsequent to a perturbation, the solution may take the form of stable, sustained oscillations (13, 14, 17). This parameter-dependent behavior probably arises from a Hopf bifurcation.

A bifurcation may occur when , the magnitude of the instantaneous gain of the feedback response, exceeds a critical value _c (13). The instantaneous gain, investigated in detail in Ref. 14, corresponds to the maximum reduction in single-nephron glomerular filtration rate (SNGFR) resulting from an instantaneous shift of the TAL flow column toward the MD, under the assumption that the SNGFR response is also instantaneous. The instantaneous gain is given by  = K₁K₂S'(1), where K₁K₂ is a measure of the strength of the feedback response, and S'(1) (a negative quantity) is the slope of the steady-state chloride concentration profile at the MD. (In a negative feedback loop, the feedback gain is negative by convention; thus we use the phrasing "gain magnitude" when referring to .)

The critical gain magnitude _c and the associated critical angular frequency _c can be determined from the model's characteristic equation, given in APPENDIX B. The parameters from Table 1 lead to critical gain _c  3.24 and critical frequency f_c = _c / (2)  45.9 mHz, in dimensional units. The critical frequency f_c predicts the frequency of the oscillations arising from the bifurcation.In this study, we assume that all parameters are fixed except for the sensitivity of the TGF feedback response, k; the gain magnitude  depends on the sensitivity through the equation K₂ = kC_o/2 (see APPENDIX A). Thus, by changing sensitivity k, we change . The baseline value of k used previously by us in Ref. 13 is 0.24 mM¹, which results in   3.14, a value below, but near, the critical value _c  3.24. We have previously hypothesized that the nearness of physiologically plausible values of  to _c may explain the tendency of some short-looped nephrons to exhibit sustained TGF-mediated oscillations while other nephrons do not (13).

View larger version (15K): [in this window] [in a new window]   Fig. 1.   Schematic representation of tubuloglomerular feedback (TGF) loop model. Thick ascending limb (TAL) is modeled by Eq. 1 as a rigid flow tube. Luminal chloride concentration C depends on the space coordinate x and the time t; Ce, extratubular chloride concentration. Actions of the glomerulus, proximal tubule, and descending limb, represented by the boxes, are modeled by phenomenological relations, given as Eqs. 2-4. C(1, t), luminal chloride concentration at the macula densa (MD); CMD, delayed signal; and F, TAL flow rate. In this study, the "Input" is broad-band forcing, which is added to the feedback response when the feedback loop is closed. "Output" is subjected to spectral analysis.

Power spectra. The fast Fourier transform (FFT) was used to obtain power spectra from numerical solutions of the model equations. A power spectrum of a signal, roughly speaking, is a graph that shows the relative strengths of the oscillatory components of differing frequency that combine to make up the signal. A power spectrum is usually represented as "power spectral density" (19); the abscissa in such a power spectrum is frequency, and the ordinate is the sum of the squares of the two Fourier coefficients corresponding to that frequency (see Eqs. C2 and C3 in APPENDIX C). We will denote the power spectral density corresponding to frequency f by P_f . Although a power spectrum is graphed as a continuous curve, in practice a computed spectrum is discrete, and there tends to be leakage from frequency "bins" corresponding to large Fourier coefficients into adjacent frequency bins (19). APPENDIX D summarizes the methods used to compute numerical solutions to the model equations and to compute power spectra based on those solutions. High resolution in space and time is required to compute solutions that are sufficiently accurate to provide acceptable spectral resolution (see DISCUSSION and APPENDIX D).

We will investigate the spectral properties of TGF system components by superimposing a broad-band forcing on a variable of the system; this variable is considered to be the input of the system. A broad-band forcing is a signal that includes many oscillatory components, of nearly uniform amplitude, over a range of frequencies. The effect of the forcing on an output of the system can be used to understand the action of the system components as a function of input signal frequency.

	RESULTS

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

Power spectra for open TGF loop. We have previously noted that both our model for the TAL and our model for the distributed delay contain a low-pass filter (13, 15-17), i.e., qualitatively speaking, low-frequency oscillations pass through these model components with little reduction in amplitude, but the amplitudes of high frequency oscillations are attenuated. To distinguish the contributions of the two filters to the TGF loop and to test the adequacy of our methods, we computed power spectra for each of the two filters separately, and then we computed the spectrum for the filters operating in series. In all these cases, the TGF loop was open, i.e., TGF-regulated flow did not enter the TAL.

The three power spectra arising from model components are shown in Fig. 2. In each case there was an input signal, containing broad-band forcing, and an output signal, as specified below; each power spectrum was computed from the output signal. The spectra in Fig. 2 (and also Fig. 3) were normalized by dividing by the power spectrum of the input signal; the power spectrum of that input signal, normalized by itself, is the constant value P_f = 1. 

View larger version (24K): [in this window] [in a new window]   Fig. 2.   Open-feedback-loop power spectra for distributed delay only, TAL only, and TAL with distributed delay, demonstrating the accuracy of the methods employed, through agreement with analytical predictions (see text). The low-pass filter of the distributed delay does not much affect the spectral properties of the model TAL for frequencies of 300 mHz and below. Nodes marked with asterisks correspond to values marked with asterisks in Table 2.

For the separate TAL filter, Eq. 1 was solved with input given as steady-state flow F = Q_op plus small-amplitude broad-band forcing. The output was given by Eq. 2, but as F(C(1, t)), i.e., the flow was determined by the TGF response function, but without the delays that would be introduced by Eq. 3. The sensitivity k was adjusted to provide a gain magnitude of  = 1. In this case, we obtained the spectrum in Fig. 2 marked "TAL only"; the general trend of decreasing amplitude as a function of frequency indicates that the TAL operates as a low-pass filter. However, the spectrum exhibits local minima, which correspond to nodes, and local maxima, which correspond to antinodes. The nodal structure is explained in the companion study (15), which shows that the range of NaCl excursions at the MD depends on the fluid transit time through the TAL, with nodes corresponding to frequencies where the oscillatory flow has a transit time to the MD that equals the steady-state TAL fluid transit time. In Table 2 we list selected nodes predicted by analytical techniques in the companion study (15) and the corresponding nodes found in this study through numerical calculation of model solutions and subsequent spectral analysis. There is excellent agreement up to 500 mHz, but there is increasing divergence as frequency increases, arising from the approximate nature of the numerical calculations. Nonetheless, there is agreement with error less than 1.3% from 0 mHz through at least 1800 mHz.

View this table: [in this window] [in a new window]   Table 2.   Nodal frequencies of TAL (in mHz)

In APPENDIX D we explain that the power spectrum value for "TAL only" at f = 0, given by P₀  0.8190, indicates a steady-state gain of the TAL that is within 0.3% of a value computed previously by other means (14). This close agreement provides confirmation that the power spectrum has been correctly computed and scaled.

The curve marked "Distributed delay only" in Fig. 2 shows the spectral response of the distributed delay of Eq. 3; in this case, the input C(1, t  _p) was replaced by broad-band forcing with mean value zero, and C_MD was taken as the output. This spectrum also indicates a low-pass filter, and this spectrum also exhibits nodes, at frequencies of (2 + n)/3 per second, for n = 0, 1, 2, ... These nodal frequencies arise because the integral in Eq. 3 vanishes when the integrand has the form (t  s  /2) × sin(2(2 + n)s/3 + ), with given by Eq. 4, for any phase shift . A different kernel function would, of course, yield a different spectral structure, since it would be composed of different Fourier components, but physiologically reasonable choices of are likely to have little qualitative effect below 500 mHz (see APPENDIX C). In this spectrum for "Distributed delay only," the response at f = 0 is P₀ = 1, because a constant signal is transmitted undiminished.

The thick shaded curve marked "TAL with distributed delay" in Fig. 2 is the power spectrum of the TAL and distributed delay acting in series. The input was the same as for the case of "TAL only," and the output was obtained from Eq. 2, via the distributed delay of Eq. 3. Thus, this spectrum is the power spectrum for the open-feedback-loop configuration of the TGF model. We see from this spectrum that the only effect of the distributed delay below 500 mHz is to attenuate the spectral power of the TAL spectrum, especially above 300 mHz.

In the companion study (15) we found that the TAL low-pass filter, in the absence of chloride backleak, exhibited 1/f scaling for frequencies larger than about 64 mHz, i.e., the amplitude of chloride excursions at the TAL decreased inversely with frequency. When the spectra for "TAL only" and "Distributed delay only" are graphed on a log-log plot (not reproduced here), the resulting plots are linear for sufficiently large frequencies, which indicates that both spectra exhibit 1/f scaling. The TAL exhibits 1/f scaling, with chloride backleak present, in about the same range as in the absence of backleak. However, the distributed delay exhibits 1/f scaling only above 300 mHz; consequently, the spectrum produced by the two filters in series exhibits 1/f scaling only above 300 mHz. Thus the combined action of the two filters does not exhibit 1/f scaling in the range that includes much of the frequency domain of TGF and myogenic autoregulation. In particular, the scaling predicted by the model cannot be the source of 1/f scaling observed in experimental records of blood pressure in rat for frequencies ranging from 0.01 mHz up to 3 mHz (10).

Power spectra for closed TGF loop. Figure 3 gives power spectra for the closed-feedback-loop configuration illustrated in Fig. 1, corresponding to increasing values of instantaneous gain magnitude , which is a measure of feedback strength. For each spectrum, the input was small-amplitude broad-band forcing added to flow entering the TAL; the output, which was subjected to spectral analysis, was the TAL flow predicted by the feedback response (see Fig. 1). In the model, the spectral characteristics of TAL flow are the same as those of SNGFR, since by assumption TAL flow is a fixed fraction of SNGFR (see APPENDIX A and Ref. 13).

View larger version (33K): [in this window] [in a new window]   Fig. 3.   Power spectra for increasing magnitudes of instantaneous gain magnitude . Shaded curves, power spectral density (PSD) for the open-feedback-loop case; solid curves, PSD for closed-feedback-loop case. Dashed line, PSD of the small-amplitude input perturbation. Numbers beside the curves represent frequencies of local extrema in PSD for the closed-feedback-loop case (except in F, where extrema for both cases are shown). Open-feedback-loop curve in A is a portion of the shaded curve in Fig. 2. Note that for gain  less than the critical gain magnitude, c  3.24, power spectra are dominated by the spectral structure of TAL NaCl transport. For  > c, the spectra are dominated by the fundamental frequency of the stable TGF-mediated oscillation and its harmonics.

In Fig. 3, dashed lines are the (normalized) power spectrum of the input; solid curves are the closed feedback loop power spectra; the shaded curves, for comparison, are the spectra for open feedback loops (exhibited as "TAL with distributed delay" in Fig. 2). In Fig. 3, A-E, the labels along the curves are frequencies corresponding to nodes or antinodes of the closed-loop spectrum; in Fig. 3F, the extrema of both spectra are labeled.

Figure 3, A-D, shows the development of a resonant frequency, emerging from the open-loop spectrum, and increasing from ~34 to ~46 mHz as the gain magnitude  is increased from 1 to 3.24. A second resonant frequency emerges from the open-loop spectrum at ~89 mHz; this frequency can be identified by detailed analysis of the characteristic equation (Eq. B1 in APPENDIX B). The antinodes in Fig. 3 are slightly to the left of the frequency midpoints between nodes.

As noted in the section describing the MATHEMATICAL MODEL, a bifurcation may occur at the critical value of gain magnitude, _c  3.24. Consequently, for  exceeding _c, the small applied perturbation elicits sustained TGF oscillations at a frequency near the critical frequency f_c  45.9 mHz. This has dramatic effects on the power spectrum, as shown in Fig. 3, E and F. Two features are particularly noteworthy. First, the power at all frequencies is greatly increased, as indicated by the upward shift in the power spectral density curve. This additional power does not arise from the small-amplitude broad-band forcing; rather, it arises from sustained flow oscillations of large amplitude. Indeed, additional numerical studies employing a transient perturbation, but no broad-band forcing, produced spectra that are almost identical to those in Fig. 3, E and F. Additional studies also showed that the transition to spectra qualitatively like Fig. 3, E and F, occurred for < 3.26, confirming, in conjunction with the result of Fig. 3D, an abrupt transition to sustained oscillatory flow localized within ~0.6% of _c.

The second noteworthy feature is the emergence of a series of harmonics of the fundamental resonant frequency of the TGF system, ~46 mHz. As  increases, the higher frequency harmonics become stronger, as can be seen through comparison of E with F of Fig. 3. The harmonics in Fig. 3, E and F, arise from the action of the nonlinear elements in the TGF system in shaping the large-amplitude oscillations.

Waveshape distortion in open TGF loop. If all elements of the TGF system were linear, then oscillations in key variables would be pure (i.e., single-frequency) sine waves, and the high-frequency components that represent distortion from a pure sine wave (i.e., the harmonics) would not be present in Fig. 3, E and F. The nonlinear elements in our model include the filter and transport characteristics of the TAL (Eq. 1) and the TGF feedback relationship (Eq. 2). The pure and distributed delays in the feedback pathway (which enter through Eqs. 3 and 4) are linear elements; indeed, the effect of the distributed delay on a sinusoidal component is to attenuate its amplitude without changing its frequency.

Examples of distortion by nonlinear elements of the TGF pathway are illustrated in Figs. 4 and 5. Columns A and B in Fig. 4 show open-feedback-loop responses to specified sinusoidal input flows. Columns C and D in Fig. 4 exhibit waveforms arising in the closed-feedback-loop case for  = 5 and  = 10, which both exceed _c. [Experimental studies indicate that steady-state, in vivo gain magnitude, which slightly underestimates instantaneous gain magnitude (14), ranges from 1.5 to 9.9 (6).]

View larger version (29K): [in this window] [in a new window]   Fig. 4.   Shaping of waveforms by components of the TGF pathway. In each column, QIN is an input single-nephron glomerular filtration rate (SNGFR), T is the TAL fluid transit time to the MD (in units of steady-state transit time to), C is the luminal chloride concentration at the MD, and Q is the output SNGFR (for the closed-loop cases, columns C and D, QIN = Q). Vertical dashed lines in rows 1-3 coincide with local extrema of transit time; horizontal shaded bars in row 1, which span the trough and the crest of the SNGFR waveform, correspond to maximum and minimum transit times, respectively. Wide-line shaded waveforms in rows 4-6 are sine waves, for comparison. See text for detailed explanation.

In each column of Fig. 4, row 1 shows the input SNGFR (Q_IN) as a function of elapsed time (the time-scale is at the base of Fig. 4). In the model, Q_IN is related to input TAL flow by a constant factor, F_IN = Q_IN. Row 2 of Fig. 4 gives TAL fluid transit time T from the TAL entrance to the MD, which is computed from Eq. 2 in the companion study (15) and which is expressed in units of the steady-state transit time t_o  15.7 s. Transit time is an important quantity since theoretical considerations indicate that chloride concentration at the MD depends largely on TAL transit time (15). Row 3 of Fig. 4 gives luminal chloride concentration at the MD. The vertical dashed lines in rows 1-3 of Fig. 4 coincide with a local maximum and a local minimum of the transit time; the shaded bars in row 1 of Fig. 4 indicate the corresponding transit-time intervals.

The thin solid curves in rows 4-6 of Fig. 4 give the TGF responses, expressed as SNGFR, arising from Eq. 2 (Q = F/), for the indicated magnitudes of gain  (C4, C6, D4, and D5 of Fig. 4 are intentionally left blank). The wide shaded curves in Fig. 4, rows 4-6, are sine functions, adjusted to match the frequency of the solid curves and to have a similar amplitude. They are provided so that distortion relative to sine curves will be apparent.

The frequency of the sinusoidal waveforms in A1 and B1 in Fig. 4 was specified at 45.9 mHz, the value of the critical frequency f_c, to permit comparison with the closed loop results of columns C and D. Column A of Fig. 4 has SNGFR excursions of amplitude Q_IN = 6 nl/min, so that SNGFR has lower and upper bounds of 27 and 33 nl/min, respectively. For an input oscillation of this small amplitude, the response of the system, at every stage, and for values of  up to 10, appears to be essentially linear, i.e., the resulting waveforms in rows 2-6 of column A of Fig. 4 all appear to be nearly sinusoidal. Nonetheless, there is a small degree of distortion, which is apparent in Fig. 4, A4-A6, and in the corresponding power spectra (see, e.g., Fig. 5A).

In Fig. 4A1 we observe that the maximum transit-time interval, indicated by the top left shaded bar, spans a trough in Q_IN, and the minimum transit-time interval, indicated by the other shaded bar, spans a crest. These transit-time intervals correspond to the maximum and minimum transit times marked by the dashed lines on Fig. 4A2. In Fig. 4A3, maximum and minimum transit times correspond to minimum and maximum MD chloride concentrations, respectively, consistent with the analysis of the companion study (15). For the special choice of the critical frequency f_c, the response waveform Q in Fig. 4A4 is in phase with the input Q_IN.

In column B of Fig. 4, Q_IN = 9, and oscillations in Q_IN are therefore bounded by minimum and maximum flow rates of 21 and 39 nl/min, respectively. As a consequence of these larger amplitude oscillations, two nonlinear features emerge. First, transit time T increases slowly relative to its rate of decrease, an asymmetry arising from the larger time interval corresponding to the maximum T, relative to the minimum T. This asymmetry, which may be seen by comparing the shaded bars in Fig. 4A1 to those in Fig. 4B1, and which is apparent in Fig. 4B2, is reflected in asymmetry in the concentration record in Fig. 4B3, which leads to the output waveform for Q in Fig. 4B4. This waveform has two particular features, arising from TAL transit, that distinguish it from the shaded sine wave: the wide crest of the wave, relative to the trough, and a rise to the crest that is slower than the fall from the crest. We call this slope asymmetry "slow up/fast down."

A second nonlinear element arising through the TGF response given by Eq. 2 is apparent in Fig. 4, B5 and B6. Because the TGF response has maximum amplitude of Q = 18 nl, the response is bounded below and above by 21 and 39 nl/min, respectively (see APPENDIX A); consequently, a "railing" effect occurs for values of MD concentrations that lead to extreme values of effective concentration C_MD through Eq. 3. In Fig. 4B5 we observe railing at the lower bound, corresponding to large MD concentration. In Fig. 4B6, for  = 10, we observe railing at both extremes, which tends to produce a square waveform.

Waveshape distortion in closed TGF loop. When the feedback loop is closed, as in columns C and D of Fig. 4, the nonlinear behavior of the system is compounded by the nonlinear feedback, leading to a broader crest and generally squarer waveshape in the waveform for TAL flow. (With the closure of the loop, the waveforms in Fig. 4, C1 and C5, coincide, as do the waveforms in D1 and D6.) However, with the increasing nonlinearity, the action of the TAL low-pass filter becomes apparent in the waveforms for transit time in Fig. 4, C2 and D2: because the filter integrates flow, the large slopes in Q_IN are reduced.

In Fig. 4C5, with the closure of the loop, the slow up/fast down characteristic is enhanced, relative to Fig. 4B5, although the distributed delay of Eq. 3 tends to reduce the degree of the effect that could be expected, given the pronounced fast up/slow down waveform in Fig. 4C3. The flattening in the crest of the waveform in Fig. 4C5 arises from the broad trough in Fig. 4C3; and small negative slope within the crest of Fig. 4C5 arises from the small positive slope in the trough of Fig. 4C3. Also, by comparison with the peaked maxima in Fig. 4C3, one sees that that railing has clipped the lower range of SNGFR flow in Fig. 4C5. Thus the case illustrated in column C of Fig. 4 represents the mixed effects of the nonlinear TAL and TGF responses.

In column D of Fig. 4, where the gain magnitude corresponds to the extreme physiological range, the clipping effect of railing dominates the transformation of MD concentration, in Fig. 4D3, to flow, in Fig. 4D6. Transit time in Fig. 4D2 has the largest maximum and the smallest minimum of the cases examined, leading to the most pronounced extrema in MD concentration, in Fig. 4D3. Compared to the other waveforms for flow in Fig. 4, the square waveform in Fig. 4D6 has the most pronounced deviation from the reference sine wave.

Effect of nonlinear elements on power spectrum. Figure 5 provides power spectra corresponding to some of the waveforms of Fig. 4. The thin solid curves in Fig. 5, A and B, are the spectra for the sinusoidal, and thus single-frequency, flow inputs Q_IN illustrated in Fig. 4, A1 and B1. The thick shaded curves in Fig. 5, A and B, are the power spectra of the TGF response, corresponding to the solid-line curves in Fig. 4, A4 and B4. Although the response illustrated in Fig. 4A4 appears to be substantially linear, the power spectrum for the response shows that the elements of the TGF pathway have introduced substantial spectral structure: the peaks in the gray curve in Fig. 5A correspond to the fundamental frequency of the input, plus a series of harmonics. When the amplitude of the input is increased, as in Fig. 4B1, the clearly nonlinear response observed in Fig. 4B4 corresponds to the more pronounced excitation of harmonics shown in Fig. 5B.

View larger version (46K): [in this window] [in a new window]   Fig. 5.   Power spectra of selected waveforms from Fig. 4. A-D correspond to columns A-D in Fig. 4. A and B: open feedback loop. C and D: closed feedback loop. A: even in the case of small-amplitude sinusoidal input QIN, significant harmonics may be generated in resulting SNGFR Q by the nonlinearities of the TGF pathway. B: harmonics become more pronounced for larger-amplitude input. C: power spectra of luminal chloride concentration C at the MD and SNGFR Q for a gain magnitude near the middle of the experimentally measured range. Harmonics apparent in Q are reduced in concentration C with each passage through the model TAL and then rebuilt by the TGF response. D: power spectra for a gain magnitude at the high end of the physiological range. Harmonics are more pronounced in this case, but still the model TAL tends to reduce their power. Labels above peaks are the fundamental frequency and its harmonics.

Figure 5, C and D, give power spectra of the MD chloride concentration C (thin solid curves), corresponding to Fig. 4, C3 and D3, and power spectra of SNGFR Q (thick shaded curves), corresponding to Fig. 4, C5 and D6. These spectra show the substantial increase in the power of the harmonics as the distortion from the sine waveform becomes more pronounced with increasing gain magnitude. Also, these spectra for the MD concentration show the action of the TAL low-pass filter in reducing the strength of the harmonics that are present in TAL flow. These harmonics are then reconstructed in the flow by the TGF response function, largely through the effect of railing.

The pattern of harmonic frequencies in Fig. 5 can be understood in terms of the Fourier components of a periodic wave. A superposition of sine curves, with varying frequencies and amplitudes, is required to represent a nonsinusoidal periodic oscillation, and if that oscillation has frequency f, then the sine curves must have frequencies of nf, n = 1, 2, 3, ..., since the oscillation is periodic. As increases, the waveforms of the TGF pathway become more distorted, with the curve segments connecting extrema becoming steeper; consequently, the high-frequency Fourier components make larger contributions to the waveform representation.

	DISCUSSION

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

We have used a mathematical model to investigate the spectral properties of the TGF pathway. For an open-feedback-loop configuration, the results of this study are consistent with the nodal TAL response pattern predicted in the companion study (15). For the closed-feedback-loop configuration, this study predicts that the spectral properties of TGF-regulated flow depend largely on whether the gain magnitude exceeds the critical gain required for the emergence of sustained TGF-mediated oscillations. The nodal structure of power spectra for subcritical gains is a consequence of the filter properties of the TAL; for sustained oscillations, power spectra exhibit a harmonic structure that arises from nonlinear properties of the model TGF pathway.

Although the results presented here are based on numerical calculations that employ a fixed parameter set (with the exception of TGF sensitivity, which is used to vary feedback gain), the study's qualitative conclusions are independent of the particular parameter choices in Table 1; indeed, the results and conclusions depend only on the structural characteristics of the model, through its dependence on steady-state transit time t_o, and on the combinations of parameters that generate the critical gain magnitude and critical frequency through the characteristic equation (Eq. B1). Thus the nodal patterns and harmonic frequency structure observed in Fig. 3 should arise for any choice of parameters in the physiological range, through a rescaling of the frequency axis.

Effects of idealizations in model formulation. Factors not included in the model formulation may affect the in vivo spectral characteristics of the TGF system. Several of these factors were considered in the companion study (15), including the elastic compliance of the tubular walls, axial diffusion of NaCl within luminal flow, spatial inhomogeneities in TAL luminal diameter and transport capacity, and the effect of time-varying luminal concentrations on transepithelial transport rate. Other factors, particular to the applicability of this study, include oscillations introduced by respiration, the dynamic properties of absorption by the proximal tubule and descending limb, and spectral characteristics contributed by spontaneous vasomotion of the renal vasculature. However, because rat respiration has a frequency of ~1 Hz, this factor is not likely to significantly affect tubular flow spectral characteristics at frequencies below 500 mHz (see APPENDIX C). Little research has been conducted on the dynamic properties of glomerulotubular balance (GTB), but existing experimental studies suggest that GTB is robust for flow variations within the physiological range (21). The effect of spontaneous vasomotion awaits further investigation.

A final factor that may impact the applicability of the results is the phenomenological characterization of the delay in feedback response given by Eqs. 3 and 4. However, theoretical considerations, developed in APPENDIX C, indicate that the spectral structure introduced by the delay in feedback, in the range of applicability of the model, will be insensitive to the precise mathematical characterization of the delay, provided that the characterization has certain general features. The insensitivity is a consequence of the short time scale of the delay, relative to the steady-state TAL transit time.

Numerical methodology. In the course of this investigation, we found that great care must be exercised to obtain power spectra that faithfully represent the nonlinear features predicted by the mathematical model. Numerical solutions to model equations must be computed with sufficient accuracy to preserve the structure inherent in the model equations, and good frequency resolution must be attained in the power spectra that are computed from the numerical solutions. The computation of power spectra, based on given data, has been heavily studied (19); the methodology used in this study is summarized in APPENDIX D.

We considered the computation of accurate numerical solutions for dynamic TAL flow in Ref. 18, where we reviewed research which shows that numerical methods may produce approximate solutions to model equations that exhibit artifactual diffusion (which redistributes spectral power) and/or artifactual dispersion of Fourier components (which displaces propagation speed of spectral components as a function of frequency). In this study we used a low-diffusion, low-dispersion method, combined with high resolution in space and time, to faithfully represent the high-frequency Fourier components in the concentration profiles of the TAL and thus preserve spectral structure (see APPENDIX D). The results for test cases were confirmed by comparison with the analytical results in the companion study (15), which predict a regular nodal pattern (cf. figure 1 in Ref. 15 and Fig. 2 in this study).

Three previous model studies that used similar formulations for the TAL (8, 9, 11) appear to have not detected the regular pattern of harmonics predicted by this study, and the waveforms exhibited in Ref. 8 do not exhibit marked nonsinusoidal features. The discrepancy between these studies and our results may be due, at least in part, to a highly dispersive numerical method that was used in the previous studies, coupled with low spatial and temporal resolution.

Comparison with published experimental data: waveforms. This model study predicts specific patterns of waveform distortion in tubular flow, which are associated with specific spectral characteristics. Are the salient characteristics of these waveforms and associated power spectra observable in vivo? The answer to this question speaks directly to the adequacy of our model to represent essential features of the TGF system. Moreover, if waveform distortion is observable in vivo, the results of this model study have important implications for the interpretation of power spectra derived from experimental records.

To determine whether the features of waveform distortion shown in Fig. 4 were present in vivo, we examined published tracings of TGF oscillations. The first feature we sought to identify is a broadening of the crest of the flow waveform, relative to the trough, which retains a more pointed shape, as seen in Fig. 4C5. This feature is most obvious at gain magnitudes of  ~ 4-5, before the waveform is large enough to be significantly constrained by the limits of the TGF response function. The second feature was a difference in the magnitude of ascending and descending slopes of the flow waveform, which was also seen in Fig. 4C5. These features, which arise from the inverse relationship between fluid speed and TAL transit time, are particularly evident in the MD chloride concentration waveforms in Fig. 4 of this study and in figure 3 of the companion study (15), before the distributed delay of AA response has acted to smooth the curves. Note that when examining concentrations, the trough, rather than the crest, is broadened, and the ascending slope magnitude is larger than descending magnitude. (In figure 3 of the companion study (15), the flatness of the trough of the curve corresponding to f = 0.5/t_o may be accentuated by solute backleak, which impairs the capability of the TAL to reduce chloride concentration at low flow rates.)

Eight experimental time records for single-nephron pressure, flow, or MD concentration were examined; the results are summarized in Table 3. Every record showed differences, in the majority of the displayed periods, in ascending and descending slope magnitudes, with the fall in the recorded variable being discernibly more rapid than the rise (except for the reverse case in distal chloride concentration). Four of the eight records exhibited broadening of the crests of the waveform, whereas the excursions through the minima were sharp. Because the broadening of the crests of the oscillatory flow waveform may be accentuated by increasing TGF gain (see Fig. 4), its appearance in some records but not others may be attributable to internephron heterogeneity.

View this table: [in this window] [in a new window]   Table 3.   Waveshape summary

A typical measured TGF oscillatory waveform of late proximal tubule pressure, from a study by Holstein-Rathlou et al. (11) (their figure 1, bottom right), is reproduced in Fig. 6B; an earlier study has shown that the waveforms of flow and pressure have similar shapes (7). In Fig. 6A we exhibit model SNGFR, for  = 4; the waveform has been scaled temporally to have the same period as the experimentally measured waveform, 28.4 s, which corresponds to a frequency of 35.2 mHz. In Fig. 6A, the dashed curve, appearing in the fourth through the sixth oscillation, is a sine wave, provided for comparison with the nonsinusoidal model waveform. In Fig. 6, the vertical dashed lines, appearing in the seventh through the tenth oscillations in both A and B, are drawn to coincide with the points on the model waveform with local maximum slope magnitude. Comparison of the waveforms in Fig. 6, A and B, indicates a close correspondence between the wave shape predictions of the model and the in vivo measurement: both waveforms exhibit broad crests and a smaller ascending slope magnitude relative to descending magnitude.

View larger version (17K): [in this window] [in a new window]   Fig. 6.   Comparison of model waveform of SNGFR (A) with experimental waveform of proximal tubule pressure (B). Model waveform, computed for gain magnitude  = 4, has been scaled temporally to have the same period as the experimental waveform. Dashed curve in A, given for the interval extending from 100 to 180 s, is a sine wave; vertical dashed lines mark times that correspond to maximum slope magnitudes in the model waveform. B: typical TGF-mediated oscillations in proximal tubule pressure. Both waveforms exhibit wide crests, relative to troughs, and small ascending slope magnitudes, relative to descending slope magnitudes. [B is reprinted from Holstein-Rathlou et al. (11).]

From this examination of experimental results, we conclude that in vivo TGF oscillations exhibit features consistent with the predictions of our model. However, the mechanisms that operate in the model may not be the sole causes of waveform distortion in vivo. For example, the slope asymmetry could also result from differences in the time constants of the turn-on and turn-off transitions of the TGF mechanism acting across the JGA, a feature not represented in the model; indeed, experiments suggest such a hysteresis in TGF responses to manipulations of Henle's loop flow between zero and high flows (figure 13 of Ref. 20). However, these results are not strictly comparable to our simulation studies, since the TGF on-transition will be principally determined by the washout of the TAL, while the off-transient will be dominated by the kinetics of TAL NaCl absorption, as the epithelium dilutes the stationary fluid column within the TAL. To a lesser extent, all commonly used experimental techniques will include, and may be influenced by, the dynamics of TAL absorption. Hence, there is no definitive evidence, at present, regarding asymmetry in the transmission of the TGF signal across the JGA; but if such asymmetry exists, it would reinforce the asymmetry arising from the transit-time dependence of TAL NaCl absorption.

Another alternative cause for waveform distortion is a displacement of the TGF response operating point to a site near the top of the response curve, which could result in the broadening of the waveform crest, but not the trough, by imposing a constraint on the flow increase allowed by the feedback response (in Eq. 2 we assume that the operating point is at the midpoint of the response curve). However, experimental evidence indicates that the TGF operating point is usually near the center of the TGF curve in normal, extracellular volume-replete rats (20), consistent with the model assumptions (13).

Comparison with published experimental results: power spectra. Regardless of the mechanism of waveform distortion, the TGF waveform in vivo has marked similarities to those shown in Figs. 4 and 6A, which suggests that the harmonics of the TGF fundamental should be present in power spectra computed from experimental data. In particular, power spectra from in vivo data should have some of the features of the spectra shown in Fig. 3, E and F, especially for low frequencies (for instance, <200 mHz), where confounding factors may be reduced.

The prediction that power spectra of TGF oscillations will be characterized by maxima at the fundamental frequency and its harmonics is supported by spectra in the experimental literature. Studies by Holstein-Rathlou and Marsh (7) and Yip et al. (22) appear to show a resonant TGF oscillation in rat proximal tubular pressure and its first harmonic, similar in this respect to our Fig. 3E. In figure 1B in Ref. 7, an oscillation in pressure of 35-44 mHz appears to have a harmonic in the range of 74-88 mHz. [The higher frequency components at 133-163 mHz in figure 1B of Ref. 7 may arise from intrinsic vascular oscillations of the AA (1, 23); also note that the plot in figure 1B of Ref. 7 uses a linear, rather than a logarithmic, ordinate, which may account for the absence of discernible harmonic structure in the power spectrum of distal chloride concentration given in the same figure.] In figure 1D in Ref. 22, an oscillation of 33-35 mHz appears to have a harmonic in the range of 66-70 mHz. In figure 1E of Ref. 22, a resonant frequency of ~23 mHz appears to have a harmonic at ~47 mHz.

Another study by Yip et al. (23), however, contains power spectra of rat efferent arteriolar flow that appear to exhibit a nodal pattern. For example, figure 5B (dashed curve) of Ref. 23 exhibits local minima at ~50, 97, 153, and 204 mHz and local maxima at 24, 66, 123, and 172 mHz, patterns that are very similar to the nodal structure observed in our Fig. 3, A-D. However, it appears that these experimental flow measurements already exhibit a resonant TGF oscillation of modest power (see figure 5A of Ref. 23), in which case comparison should be made with our Fig. 3, E-F, where there is better qualitative agreement in shape but somewhat less agreement in nodal pattern. A plausible explanation for the apparent discrepancy is that the nephron examined may have a subcritical gain magnitude, and its oscillations may be driven through coupling with a neighboring nephron that is spontaneously oscillating (5). Indeed, when we perturbed our model configuration with  = 3, using a large-amplitude signal with a sweeping frequency, from 0 to 1 mHz, we obtained a power spectrum that was in good qualitative agreement with the nodal structure in the closed-loop spectrum of our Fig. 3C and similar to the dashed-line spectrum reported by Yip et al. (23) in their figure 5B.

Finally, studies of whole kidney renal blood flow and arterial pressure by Cupples et al. (3) yielded complex power spectra that exhibit elements suggestive of both nodal patterns and TGF harmonics (see figure 1 in Ref. 3).

Based on these comparisons with experimental data, we conclude that sustained TGF-mediated oscillations in vivo can be sufficiently nonsinusoidal to exhibit substantial power in harmonics that lie above the fundamental frequency. The same is true for oscillations in nephrons with subcritical gain magnitude, since the TGF system may express the nodal structure of the TAL filter in response to perturbations. The harmonics and nodal structure may be confounding factors in studies of in vivo power spectra, where detailed analysis has been used to identify and quantify other oscillatory elements, e.g., the intrinsic myogenic oscillation and interactions between TGF and the myogenic oscillation (2, 23). High numerical resolution in computer simulations and careful consideration of experimental design and data analysis are needed in studies of TGF dynamics, because nonlinear characteristics of the system can be easily lost or obscured.

	APPENDIX A

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

Normalization of Equations

(A1)

(A2)

	APPENDIX B

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

The Characteristic Equation

(B1)

(B2)

	APPENDIX C

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

Spectral Properties of the Distributed Delay

The basis of these general considerations is that the time scale of the distributed MD delay (2-4 s; Ref. 1) is much shorter than that of the transit time of the TAL (15-20 s) or of the sustained oscillations that may be mediated by the TGF pathway (with period 20-33 s; Ref. 7). These time scales correspond to characteristic frequencies of 250-500 mHz for the MD delay, 50-60 mHz for the TAL transit time, and 30-50 mHz for the sustained oscillations. Thus spectral structure below 500 mHz (and particularly below 250 mHz) will be dominated by contributions from TAL transit time or TGF-mediated oscillations. It follows as a corollary that tubular compliance, estimated in Ref. 7 to have a characteristic time of 1 s (corresponding to 1000 mHz), will not significantly affect spectral structure below 500 mHz.

The kernel function represents the time course of the signal that modulates AA diameter. Experiments show that a step increase in MD chloride concentration produces a sigmoidal decrease in AA diameter like that shown in figure 6 of Ref. 1. To represent this sigmoidal transition, we make the following mathematical assumptions about the kernel function: the function is continuous on the real line; the function is nonnegative in an interval W = [/2, /2] of duration  and equal to zero outside the interval; and in W, the function increases monotonically from a value of zero to a maximum amplitude, near the center of the interval W, and then decreases to a value of zero (thus is symmetric, or nearly so, about the center of W). The form of the kernel specified in Eq. 4 is consistent with these assumptions.

A general kernel function, meeting these assumptions for on W, can be expressed as a Fourier series (4),

(C1)

where cos(2nu/) and sin(2nu/) are the basis functions, and where a_n and b_n are the Fourier coefficients, which are given by

(C2)

(C3)

Outside W,  0. Note that for n  1, the sine and cosine functions in Eq. C1 have periods /n, which evenly divide the interval for which may be nonzero.

Since we require that a constant input concentration C in Eq. 3 pass through the distributed delay without a change in value (i.e., C_MD must equal C if C is unchanging in time), the kernel function must have weight one on the interval W (i.e., (u)du = 1), which implies that a_o = 2/ for all choices of the kernel function (the definite integrals of all other terms of Eq. C1 vanish). A hypothetical choice of using only this constant term (i.e.,  = 1/, in W) is a worst case, because Fourier components of the input function of the form sin(2nu/ + ), for arbitrary phase shift  and for n = 1, 2, 3, ..., will be orthogonal to , and will therefore produce pronounced minima in the power spectrum at the frequencies n/. If  = 3 s, as in this study, these minima will fall at frequencies (in mHz) of 333, 666, 1000, 1333, ...

The choice of a constant function produces a linear transition in response to a step change in input C, but the choice of given by Eq. 4 produces a sigmoidal transition. The corresponding power spectrum will have minima with the same spacing, but the minima will start at 666 mHz, since only input components of the form sin(2nu/ + ) for n = 2, 3, 4, ..., will be orthogonal to . These minima occur in the curve labeled "Distributed delay only" in Fig. 2.

For the general kernel, taken to possess the assumed properties of the physiological kernel function, we expect that the cosine terms, which are even functions, will dominate the sine terms, which are odd functions, and that lower frequency terms will dominate high frequency terms. Thus, the largest Fourier coefficients will be a_o and a₁, as in Eq. 4, and consequently the most significant orthogonal cancellations will occur for input frequencies of n/, n = 2, 3, ..., and orthogonal cancellations arising from the input frequency 1/ will be small. If the distributed delay is spread on an interval of duration 4 s or less, then 2/ will be no smaller than 500 mHz, and the less significant frequency 1/ will be no smaller than 250 mHz. It follows that the conclusions of this study are unlikely to be affected by the choice of .

Now consider a more general formulation of the MD delay in which the formal distinction between the pure and distributed delays is removed. Let Eq. 3 be replaced by C_MD(t) = (t  s  /2) C(1, s) ds, with the interval W of duration  lengthened to include the duration of the pure delay, _p, and with (u) modified to be nearly zero for some portion of the left side of the interval W, for instance, [/2, /2 + _p]. Again, the kernel function can be represented as a Fourier series, with the same generic form as given in Eq. C1, but the basis functions of the Fourier series will not have periods that are multiples of the interval during which is significantly different from zero. Nonetheless, significant cancellations will only occur for the input functions that have periods that evenly divide the interval for which is significantly different from zero, and the consequences for power spectra will not differ from those already noted. Thus, the conclusions of this study are unlikely to be affected by the formal separation of the MD delay into pure and distributed components.

	APPENDIX D

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

Numerical Methods

Figures 2 and 3. The power spectra in Figs. 2 and 3 arise from an imposed perturbation of the form I(t) = (t)/I_o, where (t) =  cos (2 f_nt + (1)ⁿ/2), with f_n = t_onM/(N × 1 s), and where I_o = N × max{(t)|t  [0, p]}, with period p specified below. For Fig. 2, M = 2 and N = 2048; in this case, I(t) has dimensional period p = t_o/f₁ = 1024 s (~17 min) and a maximum amplitude of ~0.973 × 10⁴. For Fig. 3, M = 1 and N = 1024; in this case, I(t) has the same period, but a maximum amplitude of ~1.94 × 10³. The function I(t) is equivalent to a sum of sine functions with alternating phase shifts of 0 and , a construction which allows the perturbation to evolve gradually from I(0) = 0, thus avoiding the excitation of high-frequency modes that may lead to aliasing. Division of  by I_o yields a maximum amplitude of ~10³; this scaling ensures that we observe the linear response properties of the system, at least for gain magnitude below the critical bifurcation value _c.

In Fig. 2, for the case designated "Distributed delay only," C(1, s  _p) in the integrand of Eq. 3 was replaced by I(t), and C_MD(t) was considered to be the output O(t). For the case "TAL only," F was taken to be 1 + I(t), and O(t) was taken to be the dimensionless SNGFR, computed without any delay, i.e., O(t) = 1 + K₁ tanh[K₂(C_op  C(1, t))]. For "TAL with distributed delay," F was taken to be 1 + I(t), and O(t) = F(C_MD(t)), computed from Eqs. 1-4.

Equation 1 was solved using a second-order essentially nonoscillatory (ENO) scheme, coupled with Heun's method for the time advance. This algorithm yields solutions that exhibit second-order convergence in both space and time (12). The integral of Eq. 3 was evaluated by the trapezoidal rule. The numerical time and space steps were x = 1/640 and t = 1/320 = 3.125 × 10³ s, where t, here and below, is given in dimensional units. This high degree of numerical grid refinement is required for sufficiently accurate resolution of oscillations up to 1 Hz, and it provides good qualitative results up to 2 Hz (see Fig. 2 and Table 2). As shown in the companion study (15), flow oscillations produce standing waves in luminal chloride concentration along the TAL; consequently, each frequency component of the broad-band forcing having frequency greater than 1/t_o will produce one or more nodes, relative to the steady-state concentration profile, at sites along the length of TAL. To obtain valid information about the spectral properties of the model, the standing-wave components must be resolved by the numerical methods. For oscillations of 2 Hz, the nodes will be separated (see Ref. 15) by dimensional length (L/t_o)/(2 Hz)  0.0159 cm, or nondimensional length 0.0318, and the associated wavelength will be twice this length. Thus the 640 subintervals used for the TAL resolved the wavelength of this highest frequency component on a numerical grid of ~10 points.

Sampling of model output O(t) for Fig. 2 and for Fig. 3, A-D, began after one period of I(t). In Fig. 3, E -F, the perturbation I(t) allowed sustained oscillatory solutions to develop; in these cases, sampling began after two conditions were met: 1) the oscillations in F reached maximum amplitude, and 2) an integral number of periods of I(t) had elapsed. The output O(t) for all cases was sampled at 5 Hz (i.e., every 64t = 0.2 s) for 5 × 2048 points, corresponding to a real time interval of 2048 s, exactly twice the period of the perturbation I(t). The mean of the output O(t) was computed and subtracted from O(t) to prepare the data for spectral analysis.

The spectra displayed in Figs. 2-3 are estimates of power spectral density, called periodograms, which are computed from the discrete Fourier transform of the demeaned output O(t). In our implementation, the periodograms were computed via a FFT and a supplementary algorithm from Ref. 19, both adapted to double precision arithmetic. The supplementary algorithm minimizes spectral variance per data point by averaging periodograms obtained from overlapping data sets. We used four overlapping sets of 4096 points, and we chose the Welch window for the FFT. The periodograms obtained from O(t) were normalized through division by the periodogram of I.

The domain of the periodogram values corresponds to the frequencies f_n = 2f_N n/4096, where n = 0, 1, ..., 2048, and where the Nyquist frequency f_N is given by (2 × 64t)¹. Thus, domain values are spaced at intervals of about 1.221 mHz, from 0 to 2.5 Hz.

For n > 0, the ordinate values of the periodogram, P_n, approximate the sum of the squares of the Fourier coefficients (corresponding to the f_n) of the response to the input signal I(t). However, P₀ corresponds to the square of the constant term of the Fourier series, and should therefore equal the absolute value of the open-feedback-loop steady-state gain G_ss. For the parameters used in the back-leak case, we showed in Ref. 14 that |G_ss|/  0.9069. For the case "TAL only" and in each case shown in Fig. 3, we find that /  0.9050, which is excellent agreement, given that different numerical methods were used in the two studies. [Although a frequency component for zero frequency is not introduced by the perturbation I(t), a value of P₀ emerges in the spectra as a consequence of spectral leakage in the limit as frequency tends to zero.]

Figure 4. The steady-state TAL profile S(x) corresponding to flow F = 1 was computed from Eq. 1 via the ENO scheme. Normalized versions of the oscillations in Fig. 4, A1 and B1, were then introduced through F. The waveforms in rows 2-6 were recorded after at least one period of the flow oscillation to ensure that the initial profile S(x) had been expelled. The transit-time integral was evaluated by the trapezoidal rule. To elicit the oscillations in Fig. 4, columns C and D, TAL flow was initially perturbed by a square pulse (10% of steady-state flow) lasting for one transit time interval t_o. Waveforms were recorded after oscillations reached full amplitude.

Figure 5. The signals in Fig. 5 were processed as described for Figs. 2-3, including normalization by the spectrum of the broad-band perturbation used in Figs. 2-3, to allow magnitude comparisons among the figures. The spectra for MD concentration C (1, t), which corresponds to the signal that is magnified by TGF, were multiplied by 10² to permit comparison with the spectra for SNGFR.

Figure 6. The waveform in Fig. 6A, for  = 4, was computed like the waveforms in columns C and D of Fig. 4. The waveform in Fig. 6B was digitally scanned at 300 dots per inch from a reprint of Ref. 11 and stored as a PostScript file. To clearly exhibit the shape of the waveform, the image was stretched horizontally by a factor of about 2.5 (by scaling within the PostScript file), and the axes and legends were redrawn.

	ACKNOWLEDGEMENTS

We thank Chris Clausen for helpful discussions on the methodology of spectral analysis. We thank John M. Davies for assistance in preparation of Figs. 1-6.

	FOOTNOTES

This work was supported in part by National Institute of Diabetes and Digestive and Kidney Diseases Grant DK-42091.

A brief report on this work appeared in Proc. Int. Congr. Industr. Appl. Math. 3rd (Zeitschrift für Angewandte Mathematik und Mechanik, 76: 33-35, 1996).

Address for reprint requests: H. E. Layton, Department of Mathematics, Duke University, Box 90320, Durham, NC 27708-0320 (E-mail: layton{at}math.duke.edu)

Received 1 April 1996; accepted in final form 19 June 1997.

	REFERENCES

Top Abstract Introduction Results Discussion Appendix A Appendix B Appendix C Appendix D References

1.   Casellas, C., and L. C. Moore. Autoregulation and tubuloglomerular feedback in juxtamedullary glomerular arterioles. Am. J. Physiol. 258 (Renal Fluid Electrolyte Physiol. 27): F660-F669, 1990[Abstract/Free Full Text].

2.   Chon, K. H., Y.-M. Chen, V. Z. Marmarelis, D. J. Marsh, and N.-H. Holstein-Rathlou. Detection of interactions between myogenic and TGF mechanisms using nonlinear analysis. Am. J. Physiol. 267 (Renal Fluid Electrolyte Physiol. 36): F160-F173, 1994[Abstract/Free Full Text].

3.   Cupples, W. A., P. Novak, V. Novak, and F. C. Salevsky. Spontaneous blood pressure fluctuations and renal blood flow dynamics. Am. J. Physiol. 270 (Renal Fluid Electrolyte Physiol. 39): F82-F89, 1996[Abstract/Free Full Text].

4.   Farlow, S. J. Partial Differential Equations for Scientists and Engineers. Mineola, NY: Dover, 1993, p. 81-83.

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

6.   Holstein-Rathlou, N.-H. A closed-loop analysis of the tubuloglomerular feedback system. Am. J. Physiol. 261 (Renal Fluid Electrolyte Physiol. 30): F880-F889, 1991[Abstract/Free Full Text].

7.   Holstein-Rathlou, N.-H., and D. J. Marsh. Oscillations of tubular pressure, flow, and distal chloride concentration in rats. Am. J. Physiol. 256 (Renal Fluid Electrolyte Physiol. 25): F1007-F1014, 1989[Abstract/Free Full Text].

8.   Holstein-Rathlou, N.-H., and D. J. Marsh. A dynamic model of the tubuloglomerular feedback mechanism. Am. J. Physiol. 258 (Renal Fluid Electrolyte Physiol. 27): F1448-F1459, 1990[Abstract/Free Full Text].

9.   Holstein-Rathlou, N.-H., and D. J. Marsh. A dynamic model of renal blood flow autoregulation. Bull. Math. Biol. 56: 411-429, 1994[Medline].

10.   Holstein-Rathlou, N.-H., and D. J. Marsh. Renal blood flow regulation and arterial pressure fluctuations: a case study in nonlinear dynamics. Physiol. Rev. 74: 637-681, 1994[Abstract/Free Full Text].

11.   Holstein-Rathlou, N.-H., A. J. Wagner, and D. J. Marsh. Tubuloglomerular feedback dynamics and renal blood flow autoregulation in rats. Am. J. Physiol. 260 (Renal Fluid Electrolyte Physiol. 29): F53-F68, 1991[Abstract/Free Full Text].

12.   Layton, H. E., and E. B. Pitman. A dynamic numerical method for models of renal tubules. Bull. Math. Biol. 56: 547-565, 1994[Medline].

13.   Layton, H. E., E. B. Pitman, and L. C. Moore. Bifurcation analysis of TGF-mediated oscillations in SNGFR. Am. J. Physiol. 261 (Renal Fluid Electrolyte Physiol. 30): F904-F919, 1991[Abstract/Free Full Text].

14.   Layton, H. E., E. B. Pitman, and L. C. Moore. Instantaneous and steady-state gains in the tubuloglomerular feedback system. Am. J. Physiol. 268 (Renal Fluid Electrolyte Physiol. 37): F163-F174, 1995[Abstract/Free Full Text].

15.   Layton, H. E., E. B. Pitman, and L. C. Moore. Nonlinear filter properties of the thick ascending limb. Am. J. Physiol. 273 (Renal Physiol. 42): F625-F634, 1997[Abstract/Free Full Text].

16.   Pitman, E. B., and H. E. Layton. Tubuloglomerular feedback in a dynamic nephron. Commun. Pure Appl. Math. 42: 759-787, 1989.

17.   Pitman, E. B., H. E. Layton, and L. C. Moore. Dynamic flow in the nephron: filtered delay in the TGF pathway. Contemp. Math. 141: 317-336, 1993.

18.   Pitman, E. B., H. E. Layton, and L. C. Moore. Numerical simulation of propagating concentration profiles in renal tubules. Bull. Math. Biol. 56: 567-586, 1994[Medline].

19.   Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in Fortran: The Art of Scientific Computing (2nd ed.). New York: Cambridge Univ. Press, 1992, p. 542-551.

20.   Schnermann, J., and J. Briggs. Function of the juxtaglomerular apparatus: Control of glomerular hemodynamics and renin secretion. In: The Kidney: Physiology and Pathophysiology (2nd ed.), edited by D. W. Seldin, and G. Giebisch. New York: Raven, 1992, p. 1249-1289.

21.   Wilcox, C. S., C. Baylis, and C. S. Wingo. Glomerular-tubular balance and proximal regulation. In: The Kidney: Physiology and Pathophysiology (2nd ed.), edited by D. W. Seldin, and G. Giebisch. New York: Raven, 1992, p. 1807-1841.

22.   Yip, K.-P., N.-H. Holstein-Rathlou, and D. J. Marsh. Chaos in blood flow control in genetic and renovascular hypertensive rats. Am. J. Physiol. 261 (Renal Fluid Electrolyte Physiol. 30): F400-F408, 1991[Abstract/Free Full Text].

23.   Yip, K.-P., N.-H. Holstein-Rathlou, and D. J. Marsh. Mechanisms of temporal variation in single-nephron blood flow in rats. Am. J. Physiol. 264 (Renal Fluid Electrolyte Physiol. 33): F427-F434, 1993[Abstract/Free Full Text].