HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 288/5/R1143    most recent 00539.2004v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (14) Citing Articles via Google Scholar Google Scholar Articles by Marsh, D. J. Articles by Holstein-Rathlou, N.-H. Search for Related Content PubMed PubMed Citation Articles by Marsh, D. J. Articles by Holstein-Rathlou, N.-H.

RENAL HEMODYNAMICS AND CARDIORENAL INTEGRATION

Nonlinear interactions in renal blood flow regulation

¹Department of Molecular Pharmacology, Physiology, and Biotechnology, Brown University, Providence, Rhode Island; ²Department of Physics, Danish Technical University, Lyngby, Denmark; ³Department of Biomedical Engineering, State University of New York at Stony Brook; ⁴Department of Medical Physiology, Panum Institut, University of Copenhagen, Copenhagen N, Denmark

Submitted 9 August 2004 ; accepted in final form 21 January 2005

	ABSTRACT

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
We have developed a model of tubuloglomerular feedback (TGF) and the myogenic mechanism in afferent arterioles to understand how the two mechanisms are coupled. This paper presents the model. The tubular model predicts pressure, flow, and NaCl concentration as functions of time and tubular length in a compliant tubule that reabsorbs NaCl and water; boundary conditions are glomerular filtration rate (GFR), a nonlinear outflow resistance, and initial NaCl concentration. The glomerular model calculates GFR from a change in protein concentration using estimates of capillary hydrostatic pressure, tubular hydrostatic pressure, and plasma flow rate. The arteriolar model predicts fraction of open K channels, intracellular Ca concentration (Ca_i), potential difference, rate of actin–myosin cross bridge formation, force of contraction, and length of elastic elements, and was solved for two arteriolar segments, identical except for the strength of TGF input, with a third, fixed resistance segment representing prearteriolar vessels. The two arteriolar segments are electrically coupled. The arteriolar, glomerular, and tubular models are linked; TGF modulates arteriolar circumference, which determines vascular resistance and glomerular capillary pressure. The model couples TGF input to voltage-gated Ca channels. It predicts autoregulation of GFR and renal blood flow, matches experimental measures of tubular pressure and macula densa NaCl concentration, and predicts TGF-induced oscillations and a faster smaller vasomotor oscillation. There are nonlinear interactions between TGF and the myogenic mechanism, which include the modulation of the frequency and amplitude of the myogenic oscillation by TGF. The prediction of modulation is confirmed in a companion study (28).

tubuloglomerular feedback; myogenic mechanism; oscillations; vasomotion; computer simulation

We have previously used numerical simulations of mathematical models to understand the interaction between TGF and the myogenic mechanism (15, 17). The model of TGF is based on a reasonably realistic representation of tubular pressure and flow, tubular reabsorption and compliance, and of the glomerular filtration process. This model predicts oscillations in tubular flow and pressure, in renal blood flow, and in the concentration of Cl^– in tubular fluid in contact with the macula densa. The predictions are in good agreement with experimental measures of these variables.

Our previous simulations of the myogenic mechanism had been based on a model developed by Ursino and Fabbri (36). Except for the contractile process, their model makes use of empirical equations and is not based on the kind of membrane processes and intracellular signaling that are likely to be responsible for the arteriole's dynamics. Thus the only interactions between the two mechanisms that could be simulated were those directly dependent on vascular pressure in the arterioles. We have tested the combination of the Ursino model of the blood vessel and our model of TGF, but it fails to predict phase coupling or other forms of interactions. One may reasonably conclude that additional factors responsible for the interaction were not expressed in the earlier combined model.

We have therefore turned to another model of arteriolar vasomotion, derived by Gonzalez-Fernandez and Ermentrout (12) to describe the dynamics of cerebral arterioles. With suitable adjustments to the model to account for differences in the dynamical behavior of cerebral and renal afferent arterioles, we have substituted it for the Ursino model to explore specifics of the interaction.

The Gonzalez-Fernandez and Ermentrout model generates periodic autonomous vasomotion by simulating the interaction among Ca and voltage-sensitive K channels, voltage-gated Ca channels, and the membrane potential. The oscillation in Ca that it produces is used to vary the phosphorylation of myosin light chain, which in turn controls the activation of the contractile mechanism. The contractile mechanism is in series with elastic elements, and the ensemble is in parallel with additional elastic elements. The Gonzalez-Fernandez and Ermentrout model structure permits coupling to be introduced at a number of sites, but we have limited our initial work to exploring a coupling of TGF to voltage-gated Ca channels. We present the model in this paper, compare it to experimental results, and test its sensitivity to variation of parameters. In a companion paper, we show that the TGF oscillation modulates both the amplitude and the frequency of myogenic vasomotion, and we use the model to show why coupling of TGF to voltage-gated Ca channels can produce this result (28).

The presence of nonlinear interactions is an important concept in the analysis of the model results. Nonlinear interactions may be defined as those that arise either from the coupling of two nonlinear systems, where the coupling is either linear or nonlinear, or from the nonlinear coupling of two linear systems. Apart from harmonic oscillators, which are linear and unstable, the occurrence of limit-cycle oscillations in a system of equations requires that the system of equations be second-order or higher, and contain first-order dissipative terms that account for energy losses. The equations that describe both the TGF and the myogenic components of our model contain many nonlinearities. The hypothesized coupling in the model is also nonlinear. The reason for this interest is that nonlinear interactions can give rise to a number of system properties, including chaos, synchronization, and frequency modulation (32), which may be physiologically important, and which do not occur in linear systems. One manifestation of nonlinear interactions may be frequency and phase adjustments, which introduce additional peaks into the power spectrum. The model results have this characteristic. We used an analytical method called the bispectrum to test whether these peaks are the result of nonlinear interactions. The method is designed specifically for this purpose (31).

	MODELS

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
The tubule model consists of three partial differential equations with nonlinear boundary conditions. The dependent variables are pressure and flow along the tubule from the glomerulus to the macula densa, and tubular fluid NaCl concentration in the loop of Henle and early distal tubule; the independent variables are time and tubular length. Auxiliary conditions account for tubular reabsorption of solutes and water and for tubular compliance. The boundary condition at the outflow of the ascending limb calculates the tubular pressure as a function of tubular flow rate and is based on the observation that the tubular outflow resistance decreases with increasing flow rate, presumably because the increased tubular hydrostatic pressure acts to dilate the tubule (27). The boundary condition at the beginning of the tubule is glomerular filtration rate, whose value is derived from the solution of the standard model of Deen et al. (10). The glomerular capillary pressure and flow are determined by the combined action of the myogenic mechanism and TGF on afferent arterioles, which provide flow resistance to the arterial pressure. The boundary condition at the beginning of the loop of Henle sets the concentration of tubular fluid equal to the interstitial concentration.

The TGF model has been presented in earlier publications and is used here without further modification (15, 17). The equations and parameter values are presented in Appendix A, as is the GFR model of Deen et al. (10).

The Gonzalez-Fernandez and Ermentrout model (12) is based on the hypothesis that Ca and voltage-sensitive K channels and voltage-gated Ca channels interact to cause periodic action potentials. The oscillation in membrane Ca current causes an oscillation in Ca_i, which causes periodic phosphorylation of myosin light chain, and, in turn, contraction of smooth muscle. Six ordinary differential equations account for the fraction of open K channels, Ca_i, the membrane potential, the phosphorylation of myosin light chain, the length-dependent action of elastic elements, and the contraction of smooth muscle. The model calculates segmental vascular resistance as a single lumped variable varying with time. Arterial pressure is the input to the vascular resistance calculation. The local intravascular pressure also acts on the individual arteriolar segments to adjust the length of the contractile mechanism and the associated elastic elements. The resulting system is nonlinear and stiff (see Appendix A).

We made the following modifications to Gonzalez-Fernandez and Ermentrout's original formulation:

1) Three preglomerular vascular resistance segments were used: the first provides the fixed resistance of prearteriolar vasculature; the second is a site for the action of the myogenic mechanism, with some input from TGF; and the third segment receives the major input from TGF and also has a fully functional myogenic mechanism. Identical parameter values were used in the second and third segments, except that the TGF coupling term was one-third as large in the second segment as in the third, to express the length-dependent arteriolar response to the TGF signal (4, 6, 30).

2) The myogenic oscillation frequency was set to 0.12 Hz, corresponding to the experimentally observed value in renal afferent arterioles (43). Gonzalez-Fernandez and Ermentrout simulated the 1-Hz oscillation in cerebral arterioles. The adjustment was made by reducing the rate coefficient of K channel opening.

3) The reduced frequency of vasomotion required to simulate the observed frequency of vasomotion caused the arteriole's radius to increase and the vascular resistance to drop. To restore renal blood flow to experimentally measured values, we adjusted parameters in the equation that weight the contributions of elastic and smooth muscle action.

4) The reduced frequency of vasomotion also allowed Ca_i to drop further between action potentials than was the case with the parameter set used by Gonzalez-Fernandez and Ermentrout. Appropriate adjustments were made and are described below.

5) Experimental results have established the occurrence of retrograde propagation of vasoconstriction in the afferent arteriole (1, 6, 13, 14, 20, 34, 39, 42). The most probable mechanism is conduction of an electrical signal between the cells of the vascular wall (3, 13, 40). In preliminary studies with the model we observed that if the second and third segments act independently of one another and because the TGF effect is greater in the third than in the second segment, the two segments oscillate at different frequencies and become desynchronized. A small electrical coupling between the two permits them to become synchronized. This formulation is in accordance with experimental studies showing longitudinal electrical coupling between cells of the vascular wall (3, 13).

The fixed vascular resistance of the first segment, the variable resistances of the second and third segment, and fixed resistance of the efferent arteriole were used to calculate glomerular capillary pressure and plasma flow rate, which are needed for the solution of the glomerular model. The form of the equations derived by Gonzalez-Fernandez and Ermentrout was used with no other modifications.

	NUMERICAL METHODS

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
The partial differential equations describing tubular pressure, flow, and NaCl concentration were approximated as a set of finite difference equations using a centered difference scheme that is second-order correct (15, 17, 34, 38, 45). The spatial step was 0.025 cm, and the time step was 2.5 x 10^–4 s. The auxiliary equations, including the boundary conditions, were solved using the Newton-Raphson method. All computations were performed in double precision.

When solving partial differential equations by numerical methods, the numerical procedure may introduce artificial smoothing of the solution. The degree of smoothing depends among other things on the ratio of the spatial step size to the time step (x/t). In general, a larger value of x/t will be associated with a greater degree of smoothing. To test for the presence of artificial smoothing, we conducted five simulations in which the spatial step size was successively reduced by a factor of 2. The last of these simulations therefore had a step size of 7.8 x 10^–4 cm. All simulations, including the one using the standard 2.5 x 10^–2 cm spatial step, had oscillations generated by both TGF and the myogenic mechanism and that occurred at the same times in the simulations. There were minor variations in the amplitude, none systematic. The power spectra of the simulations each had 12 peaks at precisely the same frequencies as the others. The power at each peak frequency retained the same rank order in all simulations. We conclude that with the present choices of spatial step size and time step, the degree of smoothing due to the numerical procedure is negligible.

The solution of the arteriolar model of Gonzalez-Fernandez and Ermentrout for the second and third arteriolar segments was obtained with Gear's method (11). This method makes use of a variable step size combined with an explicit method such as the Adams-Moulton predictor-corrector and is well suited to stiff problems, such as this one. For the solution of the arteriolar equations, we used an initial value of 10^–6 s for the time step and solved the equations forward to the end of the current time step being used for the solution of the partial differential equations, as described in the previous paragraph. Because the solution of the partial differential equations required iterations at each time step, the initial conditions for the arteriolar equations were stored and were reused at the start of each iteration.

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

Parameters. Tables in the Appendix A contain the parameter values used in the simulations. The parameter values for the tubule and GFR models are in Tables 3–5, and are taken from (15) and (17), which contain the original literature citations and were not modified. The parameter values for the arteriolar model are contained in Tables 6–10 and are from (12), which contains the original literature citations. Changes in the parameter values are discussed in the text.

The implementation of the Gonzalez-Fernandez and Ermentrout model provides for the coupling of the myogenic and TGF mechanisms. To reduce the frequency of the myogenic oscillation from 1 Hz to 0.12 Hz, the parameter _n, the rate coefficient of K channel openings, was reduced from 2.7 to 0.3 s. This change reduced the frequency of autonomous vasomotion, which increased the time averaged radius and increased GFR and renal plasma flow. Changes made to compensate for the decreased frequency were to parameters affecting Ca_i: g_Ca, the membrane Ca conductance; k_Ca, the rate coefficient that governs Ca removal from the cell; and k, the rate coefficient that affects the rate of actin–myosin cross-bridge formation. We also adjusted the values of weighting parameters for the elastic and contractile components, w_e and w_m, used in equations A33 and A34, respectively, setting targets for the time averaged GFR and filtration fraction of 35 nl/min and 0.25, respectively, at a mean arterial pressure of 100 mmHg.

Analytical methods. Each solution of the model simulated 600 s, which required 240 s of machine time on a workstation with a 2.8-GHz processor and 1 Gb of memory. The runs were sampled at 4 Hz, yielding time series with 2,400 values. The last 2,048 of these values were used to calculate the mean plasma flow. Power spectra were computed using standard fast Fourier transform methods.

An earlier study used Volterra-Wiener kernels to show nonlinear interactions between TGF and the myogenic mechanism in whole kidney blood flow measurements made under conditions of quasi-white noise forcing of arterial pressure (7). The ability to produce such interactions is therefore an important test for the validity of the model.

In the traditional power spectrum, nonlinear interactions between oscillators with dissimilar frequencies can produce peaks that represent the sum, the difference, or both, of the frequencies of the primary mechanisms. Other peaks may represent primary activation modes, or harmonics. There are no means to differentiate among the several possible causes of such peaks (31). The presence of nonlinear interactions was therefore assessed by calculating the bispectrum, using a method described by Nikias and Petropulu (31). This method is well suited for the detection of nonlinear interactions, and is described in Appendix B. This appendix also explains why the method is useful for detecting nonlinear interactions.

	RESULTS

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
Figure 1 shows the predicted membrane K and Ca currents and the membrane potential in each of the two arteriolar segments. Currents directed out of the cell are positive; the negative Ca current spikes therefore represent an influx of Ca into the cell. Both the K and Ca currents oscillate regularly in each segment, establishing the basis for the fast myogenic oscillation. The ionic currents are slightly out of phase with each other in each segment and not of the same magnitude. The difference between them causes net changes in intracellular charge and therefore changes in membrane potential.

View larger version (42K): [in this window] [in a new window]   Fig. 1. Membrane K current, Ca current, and potential difference (PD) in the more proximal arteriolar segment (A) and the more distal segment (B).

Figure 2 compares the myogenic oscillations in membrane potential with the calculated variation in NaCl concentration in tubular fluid at the macula densa. There is periodic variation in NaCl concentration every 42.7 s, which is within the range of reported values(16, 17, 24) but longer than the value of 28.5 s that is generated by the TGF-tubule model with no active afferent arteriole (17). The periods of low NaCl concentration in Fig. 1 are correlated with periods of prolonged intervals between myogenic oscillations and reduced ionic currents; higher NaCl concentrations have the opposite effect. These correlations are consistent with the idea that the variation in membrane current flows are caused by the interaction between TGF and voltage-gated Ca channels.

View larger version (42K): [in this window] [in a new window]   Fig. 2. Membrane action potentials and the NaCl concentration in tubular fluid at the macula densa. The upper curve is the NaCl concentration, the lower one the PD.

View larger version (18K): [in this window] [in a new window]   Fig. 3. Model results: single nephron glomerular filtration rate (GFR) as a function of time. Arterial blood pressure: 100 mmHg.

View larger version (17K): [in this window] [in a new window]   Fig. 5. Model results: proximal tubule hydrostatic pressure as a function of time. Arterial blood pressure: 100 mmHg.

The top panel of Fig. 6 shows the power spectrum calculated from the results of Fig. 4. Dominant spectral peaks are at f₁ = 0.023 Hz, f₂ = 0.070 Hz; f₃ = 0.094 Hz, f₄ = 0.117 Hz, f₅ = 0.141 Hz, and f₆ = 0.168 Hz. The largest peak, f₁, is at 0.023 Hz, and can be identified with the operation of TGF. As shown in the accompanying paper (28), the TGF oscillation generated by the model is stationary and would therefore be expected to produce a single major peak in the power spectrum, possibly with some harmonics (28).

View larger version (15K): [in this window] [in a new window]   Fig. 6. Power spectra. A: calculated from data of Fig. 4 and shown in decibels. B: calculated from data of Fig. 4 and shown in cm6 x 10–12 s–3. C: calculated from a simulation with TGF fixed to eliminate the slow oscillation. Data are presented as in middle panel.

View larger version (20K): [in this window] [in a new window]   Fig. 4. Model results: single nephron plasma flow rate as a function of time. Arterial blood pressure: 100 mmHg.

The standard power spectrum is incapable, in principle, of distinguishing among spontaneously excited normal modes, coupled modes, or those due to harmonics. The bispectrum detects second-order nonlinear interactions between two harmonically related components, known as quadratic phase coupling (QPC), and these interactions result in the power at their sum and/or difference frequencies (e.g., f₃ = f₁ + f₂). Typical examples of these interactions are amplitude and frequency modulation. QPC can occur only among harmonically related components (three frequencies are harmonically related when one of them is the sum and/or difference of the other two). Thus, if the bispectrum plot shows a near-zero value at (f₁, f₂), then this suggests that any frequency modes present in the power spectrum at f₁, f₂, and f₁ + f₂, most likely represent spontaneously excited independent modes rather than quadratically coupled modes.

To gain greater insight into the processes responsible for generating the power spectrum of Fig. 6, we calculated the bispectrum of the plasma flow results of Fig. 4. The output is three-dimensional, with frequency on two coordinates, and bispectral power on the third. The spectral peaks generated by the primary oscillators lie on the diagonal of the frequency–frequency plane, and spectral peaks originating in nonlinear interactions are arrayed symmetrically off the diagonal. The results are shown in Fig. 7, upper panel. There are six significant off-diagonal peaks at the frequency pairs (0.023, 0.070), (0.023, 0.094), (0.023, 0.117), (0.023, 0.141), (0.023, 0.168), and (0.0234, 0.188), indicating the presence of nonlinear interactions. The bicoherence of these peaks is in excess of 0.99. The power spectrum in Fig. 6 shows peaks at 0.023, 0.070, 0.094, 0.117, 0.141, 0.168, and 0.188 Hz. All have counterparts in off-diagonal peaks in the bispectrum. This result shows that there are significant nonlinear interactions and that the larger peaks at frequencies above 0.07 Hz in Fig. 6 can be attributed, at least in part, to them. This result is consistent with the analysis of data from whole-kidney blood flow measurements reported by Chon et al. (7).

View larger version (51K): [in this window] [in a new window]   Fig. 7. A: bispectrum calculated from data of Fig. 4. B: bispectrum calculated from a simulation with TGF fixed. Arterial pressure: 100 mmHg.

Thus the frequency peak at 0.0703 is partially generated by the addition of the primary peak 0.0234 Hz and the harmonic peak of 0.0469. The spectral peak at 0.0937 Hz is the second largest peak in the power spectrum, and it likely represents another spontaneously excited primary mode. Although 0.0937 Hz could be a candidate for the third harmonic of 0.0234 Hz, its peak is much larger than that of the first harmonic (0.0469 Hz). Higher harmonics generally decrease in magnitude. Thus the bispectrum peak represents a nonlinear interaction between TGF (0.0234 Hz) and the myogenic mechanism (0.0937 Hz). Both 0.117 Hz and 0.0703 Hz are partially generated by the sum and difference of 0.0234 and 0.0937 Hz. The spectral peak at 0.141 Hz might arise as a second harmonic of 0.0703 Hz, but the bispectral peak at 0.0703 Hz arises from a nonlinear interaction, and the amplitude of the peak at 0.141 Hz is not much different from that at 0.0703 Hz. Thus 0.141 Hz is more likely another primary mode, and the bispectral peak at (0.0234, 0.141) represents another phase coupling between two primary modes. Two frequencies of 0.169 Hz and 0.117 Hz are consequently generated, in part by the sum and difference of the bispectrum peak at (0.0234, 0.141). We conclude that there are three spontaneously excited primary modes 0.0234, 0.0937, and 0.141 Hz, and these peaks interact to generate other peaks observed in the spectrum, as well as the harmonic frequency at 0.0469 Hz. Wavelet analyses of the data in Fig. 4 are presented in the companion paper (28). The result shows a TGF peak at 0.0234 Hz and frequency modulation of the myogenic oscillation, which switches between 0.0937 Hz and 0.141 Hz in rhythm with the TGF oscillation. The results of the two different types of analysis are therefore consistent with each other.

As a further test of the suggestion that the complexity of the power spectrum originates in nonlinear interactions, we performed a simulation with TGF input clamped to eliminate such interactions. The power spectrum is shown in Fig. 6, bottom. For a more convenient comparison, the power spectrum shown in Fig. 6, top is redrawn to show the absolute value of the power, rather than the logarithm of its ratio, as in the Fig. 6, top. This redrawn spectrum occupies Fig. 6, middle. The power of the myogenic oscillation is comparable in the two cases represented in the middle and lower panels. With TGF clamped, however, the power is confined to a single peak at 0.11 Hz and a harmonic at double the frequency. This finding is consistent with nonlinear interactions as the source of the additional spectral peaks in Fig. 6, top.

We also calculated the bispectrum on the results with TGF clamped. The result is shown in Fig. 7, bottom, which is drawn to the same scale as Fig. 7, top. In the absence of TGF oscillations there are no off diagonal peaks, which is also consistent with the suggestion that nonlinear interactions cause the additional peaks in the spectra.

The structure of the model does not permit the myogenic mechanism to be disabled while TGF is left intact to cause oscillations in blood flow. We therefore exercised an earlier model, which couples TGF to the afferent arteriole linearly, with no representation of a myogenic mechanism (17). The spectrum of the predicted renal plasma flow rate showed a single peak, with harmonics at twice and three times the frequency of the primary peak, and power at –29 and –34 dB, respectively, compared with the primary peak. These values are small compared with the peaks of Fig. 6, top, in the frequency range 0.07–0.17 Hz. The bispectrum showed no off-diagonal peaks except for one at the primary frequency, which indicates an interaction between the primary mode and the first harmonic, as detailed above. Thus, when either TGF or the myogenic mechanism operates alone, the power spectrum becomes simpler, and many of the peaks, which we have interpreted to arise from interactions, disappear.

The ratio of period lengths, as can be seen in Figs. 3–5 is 1:5. This ratio is the same as in experimental results from single nephron blood flow measurements (35).

Parameter Sensitivity

The model presented in this paper is the concatenation of three other models: the tubular model of Holstein-Rathlou and Marsh (15), the GFR model of Deen et al. (10), and the arteriolar model of Gonzalez-Fernandez and Ermentrout (12). Each model has undergone testing independent of this study, including tests of parameter sensitivity. Rather than repeat that work, we have elected to concentrate on the new parameter added to merge the models, and those parameters that required change, which are exclusively in the arteriolar model.

The new parameter is , which couples tubuloglomerular feedback to the myogenic mechanism. Fig. 8 shows the results of 108 separate simulations, in which arterial pressure and were varied. The results depicted are of the power spectral calculations over the frequency band 0.015–0.05 Hz. The spectral peak of the slow oscillation is typically spread over three or more frequencies in a 2,048 point power spectrum. Representing the result of a simulation by a single number therefore requires a summation of the spectral power over a frequency band. The results are presented logarithmically.

View larger version (33K): [in this window] [in a new window]   Fig. 8. Log spectral power of plasma flow rate in the frequency band 0.015–0.05 Hz, as a function of the blood pressure, and of the coupling coefficient, .

With set to 0.150, there is no significant spectral power at the frequencies of the TGF oscillation at any blood pressure from 80 to 115 mmHg. As is increased, oscillations first become prominent at the higher blood pressures. With at 0.275, oscillations are present at all blood pressures from 95 to 135 mmHg. Oscillations are found experimentally over this blood pressure range. Still higher values of induce TGF oscillations at 90 mmHg, and increase the amplitude of the oscillations at the higher blood pressures.

With = 0.25 and the arterial blood pressure at 95 mmHg, the amplitude of the predicted oscillation in tubular pressure is around 0.8 mmHg, or 7% of the mean value. Because of the various sources of noise in the tubular pressure signal, it is possible that an oscillation of this magnitude would not have been detected in records from experiments. A value of 0.275 for generates tubular pressure oscillations of 1.2 mmHg, or 10% of the mean value, with arterial pressure at 95 mmHg. Values of greater than 0.275 begin to produce larger oscillations of tubular hydrostatic pressure and macula densa NaCl concentration than are found experimentally at normal blood pressures (16), although there is no abrupt transition. The value of 0.275 for the coupling parameter provides oscillations of GFR, renal plasma flow, tubular hydrostatic pressure, and macula densa NaCl concentration, and allows the model to predict appropriate amplitudes of these oscillations. The increases in spectral power with increasing values of are gradual so that the choice of the value of is not critical. We have chosen to use 0.275 because it is the lowest value that is consistent with the body of experimental measures of renal blood flow dynamics.

Autoregulation of GFR and renal blood flow is an important whole system function provided by TGF and the myogenic mechanism. We tested the ability of the model to achieve autoregulation, which is usually measured while attempting to hold blood pressure constant. Therefore, we took the average value of GFR and renal plasma flow over 400 s. Fig. 9 shows the GFR as predicted by the model over a range of blood pressures from 80 to 135 mmHg. The result is shown for several values of . At the values of used for the simulations in Figs. 1–5, the GFR is invariant over this range of blood pressures. The same is true at higher values of . Autoregulation is less efficient at lower values of the coupling parameter, where TGF is weaker.

View larger version (17K): [in this window] [in a new window]   Fig. 9. Mean single nephron GFR as a function of arterial blood pressure, for five different values of , the coefficient coupling TGF to the myogenic mechanism.

Table 1 shows the results of variation in k and blood pressure on the spectral power of the TGF oscillation. The magnitude varies directly with the parameter value and the blood pressure. The value 0.011 was chosen because it yielded oscillations of tubular pressure and of macula densa NaCl concentration that were closest to measured values at normal arterial blood pressure.

View this table: [in this window] [in a new window]   Table 1. Model results: effect of variation in k on spectral power of the TGF oscillation at different blood pressures.

View larger version (20K): [in this window] [in a new window]   Fig. 10. Myogenic spectral power and nephron plasma flow rate as a function of arterial blood pressure. A: with TGF input to the afferent arteriole fixed, disabling TGF. B: with TGF intact, but using a low value of , the coupling coefficient, to reduce the modulating effects of TGF on the myogenic mechanism. The myogenic frequency in the left panel was 0.11 Hz at all blood pressures; in the right panel, it was 0.11 Hz at the lower blood pressures and increased to 0.12 Hz at 105 mmHg and higher.

View this table: [in this window] [in a new window]   Table 2. Model results: effect of variation in k on the fast oscillation, with = 150, and the blood pressure at 100 mmHg

The results presented to this point are from studies with constant arterial pressure. Although this is a necessary starting point for the presentation of the model, it is far from the reality confronting the renal circulation in live animals, where there are significant fluctuations with a 1/f pattern in blood pressure (19). We tested the robustness of the model by using natural blood pressure records from our previous study (19).

The results of a single simulation are shown in Fig. 11. The blood pressure is from a conscious unrestrained rat, outfitted with telemetering equipment. The blood pressure record was passed through a low-pass filter to remove beat-to-beat fluctuations, and was then sampled at 3 Hz. Compared with the result shown in Fig. 4, the predicted plasma flow rate shows random fluctuations, which may be attributed to the variation in the arterial pressure input, but it maintains a stable blood flow despite long-term variation in blood pressure. The blood pressure in this experiment varied over the range 95–110 mmHg, with occasional transients to more extreme values. As can be seen in Fig. 11, nephron plasma flow rate remained relatively constant despite this variation in the input, the main response to the blood pressure variation being a change in the amplitude of the slow oscillations. We interpret the variation in the amplitude of the oscillation to reflect increased activation of TGF at higher blood pressures. A similar pattern may be seen in the tubular pressure record. The power spectrum from the simulation of Fig. 11 is shown in Fig. 12, top. The power spectrum is much the same as in Fig. 6, with a single slow oscillation and additional power distributed over the frequency band 0.09–0.2 Hz.

View larger version (31K): [in this window] [in a new window]   Fig. 11. Measured arterial blood pressure as the input to the model. A: mean arterial blood pressure. B: predicted renal plasma flow rate. C: predicted tubular hydrostatic pressure.

View larger version (20K): [in this window] [in a new window]   Fig. 12. Power spectra from using measured arterial pressure as the input to the model. A: result of a single simulation using the time series of Fig. 10 (B). C: average result of nine simulations. Dashed lines show ± 1 SE.

	DISCUSSION

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
We presented a model of the regulation of renal blood flow that successfully predicts several known characteristics of the renal circulation:

Autoregulation of GFR and renal blood flow over the pressure range 80–135 mmHg. The model predicts little variation of these variables over this range. Reducing the value of the coefficient coupling TGF to the myogenic mechanism makes autoregulation less effective. This result supports a role for TGF in autoregulation.

The model contains no membrane-sensing site for pressure in the myogenic mechanism. The contractile mechanism has been given length dependence, a well-known property of vascular smooth muscle, and this length dependence adjusts the force of contraction to the length of the contractile mechanism. Thus the force of contraction is related ultimately to the arterial pressure. The length dependence therefore confers an intrinsic sensitivity to hydrostatic pressure.

Mean values of GFR, proximal tubule hydrostatic pressure, and macula densa NaCl concentration comparable to those measured experimentally (16). The measurements of single-nephron plasma flow rate, made with laser-Doppler technology, provide only fractional changes relative to an undefined mean value. The model's prediction of GFR, when divided into a standard value for whole kidney GFR, correctly predicts the number of nephrons in the rat kidney. The model also predicts a value for the filtration fraction that conforms to experimental results (16). Thus the predicted value of nephron plasma flow rate should correspond to the flow rate in rats.

Oscillations of GFR and tubular flow, renal blood flow, tubular pressure, and NaCl concentration at the macula densa. The frequency of these oscillations is in the range of those that have been measured experimentally. The magnitude of the GFR oscillation is comparable to those of the oscillation in early proximal tubular flow that has been measured. There are no measurements of the absolute value of the instantaneous single nephron blood flow for comparison. The tubule pressure oscillation is approximately ±10% of the mean value, which is the same as measured values.

The tubular segment between the macula densa and the early distal tubule, where measurements of NaCl concentration are first possible, reabsorbs NaCl and is water-impermeable. These fluxes will alter the amplitude of the oscillation in NaCl concentration, and direct comparisons between the predicted NaCl concentration and experimental result are therefore not possible. The phase relationships among GFR, tubular hydrostatic pressure, and early distal NaCl concentration are the same as have been measured (16). The phase relationship between nephron blood flow and tubular pressure oscillations is also the same as has been measured experimentally (43).

Oscillations representing spontaneous vasomotion of the afferent arterioles. The predicted oscillations are the same magnitude relative to the TGF oscillation, and the ratio of the frequencies is the same as in rats (43). Because both TGF and the myogenic mechanism carry out their actions via a single contractile mechanism, an interaction between them is inevitable and is present in the output of the model. This conclusion is reinforced by the results of bispectral analysis, which gives results similar to those of Volterra–Wiener analysis applied to whole kidney blood flow measurements (7). This higher frequency oscillation is propagated through the proximal tubule, but is absent from the predicted flow at the macula densa, a result that is consistent with the function of the thick ascending limb as a low-pass filter, as has been described (17, 22, 23, 34).

Modulation of the frequency of the myogenic oscillation by the TGF oscillation. Not previously known in the renal or other circulations, this phenomenon has now been identified in single-nephron blood flow measurements by the application of wavelet analysis (9, 28). The modulation occurs if the coupling of TGF to the myogenic mechanism is via a Ca channel in the cell membrane, because these voltage-gated mechanisms are essential elements in the generation of the oscillations of spontaneous vasomotion.

The amplitude but not the frequency of the TGF oscillation varied directly with arterial pressure. The same result was observed in those studies that used measured arterial blood pressures as an input to the model. Yip et al. measured nephron blood flow at two different arterial pressures and found that the power in the TGF oscillation increased at a higher blood pressure (43). The model result is therefore consistent with experimental observation. We interpret the increase in the amplitude of the TGF oscillation to reflect increased activation of the TGF mechanism by higher pressures.

The model predicts no TGF oscillations at arterial pressures lower than 95 mmHg, which is consistent with experimental observation. Barfred et al. (2) and Layton et al. (21) have used simple models to show the presence of a Hopf bifurcation in TGF, using TGF gain as the bifurcation parameter. Our results are consistent with this result, since the coupling parameter combines with the feedback gain in equations A22 and A23. Although the blood pressure is not, strictly speaking, a parameter of the model, appropriate values are needed to permit the bifurcation, as shown in this paper and in experiments.

The response of the myogenic mechanism to increased arterial pressure is more complex because TGF modulates both the amplitude and the frequency of the myogenic oscillation (28). The modulation and the nonlinear interaction it reflects causes the myogenic oscillation to generate several smaller peaks in the power spectrum in the frequency range 0.070–0.188 Hz, as seen in Fig. 6, rather than a single one. This redistribution of spectral power makes it difficult to estimate either the frequency or the amplitude of the myogenic oscillation. At arterial pressures in the range of 80–90 mmHg, there is no TGF oscillation, and the amplitude of the myogenic oscillation varies directly with arterial pressure. Modulation effects appear at higher pressures because the TGF oscillation emerges, and the interaction becomes more pronounced.

The two oscillations are different in the sense that the rapid one is driven by a discrete clock consisting of the K, Ca, and leak currents and their interaction. Deliberate changes in the parameters governing the K and Ca currents cause predictable changes in the frequency of spontaneous vasomotion. The TGF oscillation has no clock but arises autonomously in the operation of a negative feedback loop with a long delay, provided mainly by the reabsorption of NaCl in the thick ascending limb of Henle's loop (17).

The model merges three other models. This merger required one new parameter, , coupling TGF to the myogenic mechanism. We selected a single value as an optimum choice, but the results of Fig. 8 indicate that there is a range of values that could provide a satisfactory fit of the output to experimental measurements.

Blood pressure fluctuates (19), and autoregulation protects the distal tubule and its hormone-dependent epithelial transporters from the resulting load variations. Recorded blood pressure records used as an input provide some insight into the system's operation. We applied this input; blood flow rate was maintained at a constant mean value despite changes of blood pressure of 20–25 mmHg. The amplitude of the TGF oscillations varied with the mean blood pressure, but the mean flow rate did not change. This result shows that the model reproduces autoregulation over a frequency range of blood pressures where the mechanisms are known to act in animals. The result also illustrates the adaptations of TGF to the blood pressure fluctuations, but the response of the myogenic mechanism is obscured by the 1/f noise of the blood pressure, and by frequency modulation. In simulations at fixed blood pressure and low values of , the amplitude of the myogenic oscillation varied directly with arterial pressure, as it does in experiments (43).

In the complete absence of TGF the model predicts that autoregulation is less effective, as may be seen in Fig. 10, left. This result may be at variance with some experimental results. For example, Moore (29) showed that interrupting tubular flow to the macula densa in rat kidneys attenuated but did not eliminate autoregulation of single nephron GFR. Moore and Casellas (30) found a similar result in the perfused juxtamedullary nephron preparation using furosemide perfusion to inhibit TGF. The residual autoregulation might be due to the myogenic mechanism, but it could also be due to vascular communication from nearby nephrons with intact TGF (6, 20, 39, 42). The model presented in this paper is of a single nephron, so we cannot use it to resolve this question. Pressure-dependent vasoconstriction is found in isolated perfused afferent arterioles, which are freed from any TGF signal (44). Therefore, it is likely that the myogenic mechanism is capable of providing some autoregulation, which the model does not predict, even though the model with intact TGF predicts perfect autoregulation, as shown in Fig. 9.

Conclusion. We developed a simulation of the renal circulation based on models of tubuloglomerular feedback and the myogenic mechanism. The model successfully simulates the known dynamical behavior of that microcirculation, and predicts that TGF modulates the frequency of the myogenic oscillation. Frequency modulation was previously unsuspected in this or other circulations. Reanalysis of records of nephron blood flow confirm the existence of the phenomenon. In the accompanying paper we present the results of that analysis and use the model to illustrate the mechanism of modulation (28).

	APPENDIX A. THE MODEL

TOP ABSTRACT MODELS NUMERICAL METHODS RESULTS DISCUSSION APPENDIX A. THE MODEL APPENDIX B. THE BISPECTRUM REFERENCES

 
Glossary

Subscripts Ca_i(Ca_i, m,v,t

A

afferent arteriolar

E

efferent arteriolar

GC

glomerular capillary

I

interstitial

S

NaCl

T

renal tubule

Independent Variables Ca_i(Ca_i, m,v,t)

t

time, s

z

tubule position, cm

z_GC

glomerular capillary position, cm

Dependent Tubular Variables Ca_i(Ca_i, m,v,t)

C_S(z,t)

NaCl concentration, mM

P(z,t)

pressure, mmHg

Q(z,t)

flow, cm³ s^–1

Dependent Glomerular Capillary Variable Ca_i(Ca_i, m,v,t)

C(z_GC)

glomerular capillary protein concentratin, g/l

Auxilliary Variables Ca_i(Ca_i, m,v,t)

J(z,t)

solute (mmol/min) or volume (nl/min) flux

R(t)

vascular resistance, dyn·s·cm^–5

r(z,t)

tubular radius, cm

(C,z_GC)

plasma oncotic pressure, dyn/cm²

Dependent Arteriolar Variables Ca_i(Ca_i, m,v,t)

Ca_i(Ca_i, m,v,t)

intracellular Ca concentration, mM

n(v,Ca_i,t)

fraction of open K channels

v(Ca_i,n,t)

membrane electrical potential difference, mV

x(x,P_A,t)

length of parallel eleastic element, cm

y(y,,,t)

length of the contractile element, cm

(,,t)

fraction of actual: possible actin-myosin cross bridges

Auxilliary Arteriolar Variables Ca_i(Ca_i, m,v,t)

m(v)

equilibirum distribution of Ca open channel states

n(v, Ca_i)

equilibrium distribution of K open channel states

u(x,y)

length of series elastic component, cm

(Ca_i)

fraction of phosphorylated myosin light chain

The TGF–Tubule Model
Tubular model. The equations describing the variation in pressure, P_T, and flow, Q_T, for a noncompressible Newtonian fluid in a compliant, reabsorbing tubule at low Reynolds number are

(A1)

and

(A2)

where is the fluid density and is the viscosity (Table A1).

(A3)

The proximal tubule and the descending and the ascending limbs of the loop of Henle are each represented separately.

Proximal tubular fluid reabsorption, J_v is

(A4)

For the descending limb of Henle's loop, J_v is a function of the transtubular gradient of solute:

(A5)

Please note that to reduce the number of constants, Lp is defined here to include the universal gas constant and the absolute temperature and is defined as moles per second.

In the ascending limb of the Henle's loop, volume reabsorption is always set equal to zero. The tubular radius is expressed as a linear function of transtubular pressure:

(A6)

NaCl is treated as a single solute, s. For a solute s in a reabsorbing tubule with an axial flow, Q_T, mass balance requires that

(A7)

where A is the tubule cross section area and C_S is the solute concentration.

Because the proximal tubule is an isosmotic reabsorbing epithelium for which NaCl is the principal osmotically active component, C_s is assumed to be constant in that segment and Equation A7 is solved only for the loop of Henle. For the initial condition, C_s is set equal to the interstitial concentration, 150 mM, at the beginning of the descending limb. J_s is specified as a passive transport term in the descending limb of Henle's loop and as a sum of a passive transport term and an active term following Michaelis-Menten kinetics in the ascending limb of Henle's loop.

Thus

(A8)

where V_max =0 for the descending limb, since there is no active transport in that segment.

The interstitial NaCl concentration, C_i, varies linearly from 150 mmol/l at the cortico-medullary boundary to 300 mmol/l at the bend of the loop of Henle. For the cortical thick ascending limb of the loop of Henle, C_i was set to 150 mmol/l.

Tubuloglomerular Feedback–Afferent Arteriolar Model Following convention, the action of tubuloglomerular feedback (Table A2) is described by a logistic equation:

(A9)

where _max is the maximum obtainable response, is the dynamic range, C_S(md) is the NaCl concentration at the macula densa, k is the feedback gain, and C_1/2 is the NaCl concentration that gives the half-maximum response.

(A10)

The plasma flow, Q_A, to the glomerular capillaries is given by the conventional flow equation:

(A11)

Glomerular capillary hydrostatic pressure, P_GC, is given by

(A12)

which is obtained by assuming that afferent blood flow less the glomerular filtrate rate, Q_T(0), passes through the efferent arteriole and that thereafter the tubular reabsorbate is added to the efferent blood flow, which then passes through a distal resistance, R_v, to the venous compartment, in which the hydrostatic pressure is assumed to be zero.

(A13)

where z_GC is the fractional position along the glomerular capillary, K_f is the filtration coefficient, and L is the length of an idealized glomerular capillary. The initial condition, P_GC(0), is derived from the arterial pressure and the calculation of afferent arteriolar vascular resistance, as detailed below.The plasma oncotic pressure, (C), is found from the Landis-Pappenheimer equation:

(A14)

according to Deen et al. (10).

Myogenic Model The preglomerular vascular bed is modeled as three segments in series. The first segment represents the larger renal vessels which do not contribute to autoregulation, and whose resistance R_A1 is therefore assumed to be constant. The second and third segmen