## Abstract

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

glomerular filtration rate and renal blood flow are regulated at the level of the single nephron by an intrinsic action of afferent arterioles, known as the myogenic mechanism and by a signal conducted to afferent arterioles in response to changes in tubular fluid electrolyte concentrations, chiefly Cl^{−}, at the macula densa. The second mechanism is tubuloglomerular feedback (TGF). The two mechanisms interact to provide effective autoregulation over a range of blood pressures from ∼75 mmHg to at least 125 mmHg. An oscillation of renal blood flow has been identified with each mechanism. The myogenic oscillation has a period of 6–10 s (4, 8, 26, 35, 43); the TGF oscillation is 20–50 s long (16, 17, 24, 25). The TGF oscillation is the larger of the two(43). Each oscillation reflects the operation of a nonlinear system. Activation of TGF causes retrograde propagation of vasoconstriction (4–6, 14, 20, 30, 42), and the two oscillations show coupling (35). Retrograde propagation indicates either hemodynamic or some other form of interaction, and the observation of phase coupling between the two oscillations suggests significant intracellular interaction between the two mechanisms. Arteriolar cells have a single contractile mechanism, and the convergence of two signaling mechanisms requires that they interact. We hypothesize that the two mechanisms are coupled at voltage-gated Ca channels in the plasma membrane, and we develop a computer model to predict the consequences of this hypothesis on the dynamical behavior of the system regulating renal blood flow.

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

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

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 × 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 × 10^{−4} cm. All simulations, including the one using the standard 2.5 × 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

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.

Figure 1 also shows periodic variation in the magnitude of the ionic currents, the magnitude of the action potentials, and the frequency of the action potentials. In the model, TGF acts on voltage-gated Ca channels and affects the Ca current. Since the Ca current is an essential element in the mechanism that generates the myogenic oscillation, TGF should affect the interval between action potentials. The amplitude-modulating effects of TGF are more pronounced in the distal afferent arteriolar segment. This difference is expected because the coupling was made stronger in that segment, to reflect the spatial distribution of the TGF action. The frequency modulation is the same in the two segments because they are electrically coupled. Experimental evidence for frequency modulation of the myogenic mechanism by TGF is presented in the accompanying paper (28).

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.

Figures 3–5 show the predicted behavior of glomerular filtration rate, renal plasma flow, and proximal tubular hydrostatic pressure. Cl^{−} concentration in tubular fluid at the macula densa, the system variable sensed by the feedback mechanism, is presented in Fig. 2. An arterial blood pressure of 100 mmHg was used in the simulation. The mean values were: single-nephron GFR, 33.6 nl/min; single-nephron plasma flow rate, 136.9 nl/min; proximal tubule hydrostatic pressure, 10.9 mmHg; and NaCl concentration in tubular fluid at the macula densa, 36.1 mM. Whole-kidney GFR in the rat is ∼1 ml/min. Dividing this latter figure by the predicted nephron GFR yields a value of 29,762 nephrons; 30,000 is the generally accepted value. The predicted filtration fraction is 0.247; the usual fraction is ∼0.25. We reported mean proximal tubule hydrostatic pressure of 12.8 ± 0.3 mmHg, and Cl^{−} concentration in the early distal tubule of 39.8 ± 2.2 mM (16). The agreement between the model's predictions and experimentally measured values is good.

Figs. 2–5 show that each of the variables oscillates at 0.023 Hz, which is consistent with results of a number of experimental studies (14, 16, 18, 24, 25, 41–43). Proximal tubule pressure leads Cl^{−} in the early distal tubule, as it does in rats (16). The myogenic oscillation is present as a second, lower-amplitude oscillation at frequencies above 0.1 Hz. It is prominent in the single-nephron GFR and plasma flow rate results. The myogenic oscillation is also visible in the tubular pressure result, although it is considerably damped, compared with the GFR and plasma flow rate results. The faster oscillation is not readily apparent in the macula densa NaCl concentration results. We assume that tubule compliance, as expressed in equations A2, A3, and A6, acts as a low-pass filter and removes the faster oscillation. Two oscillations have been found in experimental measurements of nephron blood flow, using standard spectral methods for analysis (43), and in tubular pressure records, using wavelet transforms (35). Frequency modulation causes the frequency of the myogenic oscillation to vary periodically with TGF. The average frequency is 0.12 Hz.

The top panel of Fig. 6 shows the power spectrum calculated from the results of Fig. 4. Dominant spectral peaks are at *f*_{1} = 0.023 Hz, *f*_{2} = 0.070 Hz; *f*_{3} = 0.094 Hz, *f*_{4} = 0.117 Hz, *f*_{5} = 0.141 Hz, and *f*_{6} = 0.168 Hz. The largest peak, *f*_{1}, 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).

There are several processes that could be responsible for the spectral peaks *f*_{2} − *f*_{6}, including the myogenic oscillation, nonlinear interactions among the various peaks of the myogenic oscillation and the TGF oscillation, and harmonics. Using whole-kidney blood flow measurements made in experiments with white noise forcing of arterial blood pressure, Chon et al. (7) detected nonlinear interactions between the slow and fast oscillations, and these would be expected to produce peaks in the standard power spectrum.

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*_{3} = *f*_{1} + *f*_{2}). 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*_{1}, *f*_{2}), then this suggests that any frequency modes present in the power spectrum at *f*_{1}, *f*_{2}, and *f*_{1} + *f*_{2}, 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).

The bispectrum peak at (0.0234, 0.0234) reflects a high degree of phase-coupling between the peak at 0.0234 and the peak at 0.0234 + 0.0234 = 0.0469, indicating that the frequency mode at 0.0469 is the second harmonic of 0.0234 and not a spontaneously excited mode. The bispectrum peak at (0.0234, 0.0469) arises because of the relationship 0.0234 + 0.0469 = 0.0703.

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.

In Fig. 8, values of TGF power greater than −0.91 represent significant oscillations. All simulations with lower values were assigned a value of −1.0 to reduce clutter in the figure.

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.

Gonzalez-Fernandez and Ermentrout formulated their model to simulate the behavior of the rat cerebral artery, which has a myogenic oscillation at 1 Hz (12). The afferent arteriole oscillates at 0.10–0.15 Hz (43). We reduced the frequency by changing φ_{n}, a parameter in equation A20 that determines the fraction of open K channels. The need to reduce the frequency of the myogenic oscillation changed the dynamics of Ca_{i}, and therefore of both oscillations. We found it necessary to reduce the value of the Ca conductance, *g*_{Ca}, by 16% and to increase the value of the rate coefficient for Ca removal, *k*_{Ca}, by 25%. The largest change was a reduction in the value of *k*_{ξ} from 3.3 to .011. The parameter *k*_{ξ} governs the rate of cross bridge formation between actin and myosin in equation A27.

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.

To estimate the effect of variation of *k*_{ξ} on the fast oscillation, it is necessary to clamp the macula densa NaCl concentration to provide a time invariant TGF signal. This is most easily accomplished by setting the value of ζ, the TGF signal, to an arbitrary fixed level, which eliminates the TGF oscillation. The effects of clamping TGF on the spectrum and bispectrum are described above. This modification permits us to study the myogenic mechanism as it varies in response to blood pressure, without an accompanying change in the TGF signal. Figure 10 shows the spectral power of the myogenic oscillation and the nephron plasma flow rate with the TGF signal fixed. The TGF signal to the afferent arterioles varies during the TGF oscillation between values of −0.4 and 0.4. A value of −0.25 was used in the simulations shown in Fig. 10, *left*. The power of the myogenic oscillation varies inversely with blood pressure under these circumstances, as has also been reported in the mesenteric circulation (37), and plasma flow rate increases with blood pressure, indicating less effective autoregulation. Increased blood pressure has been shown experimentally to increase the power of the renal myogenic oscillation (43). The model can reproduce this result when TGF is active, even when there are no TGF oscillations. This result can be seen in Fig. 10, *right*, which shows the peak spectral power of the myogenic oscillation when ζ, the TGF–myogenic coupling coefficient, is assigned a value of 0.15. This is a relatively low value and does not elicit TGF oscillations except at blood pressures in excess of 120 mmHg. Using this low value of ζ reduces the nonlinear interactions that generate the complex power spectrum shown in Fig. 6. Figure 10 shows that the power of the myogenic oscillation increases with arterial pressure, provided that tubular flow adjusts macula densa NaCl concentration and sends the corresponding signal to the afferent arteriole.

Table 2 presents the effects of parameter variation on the frequency and spectral power of the fast oscillation with ζ = 0.15, and at an arterial pressure of 100 mmHg. Varying the parameter had no effect on the frequency with the greatest spectral power, but increasing the parameter increased power. Thus the effect of increasing *k*_{ξ} is to increase the power of both oscillations.

### Model Robustness

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.

We performed nine simulations using three different 600-s segments of the recorded blood pressure of three different animals. The average power spectrum calculated from these time series is shown in Fig. 12, *bottom*. The average power spectrum has the same features as the spectrum from the single simulation. This series of results indicates that the model is capable of responding to variation in the arterial pressure as it occurs naturally. As the system responds, the magnitude of the TGF oscillation adjusts to the blood pressure, but the mean plasma flow rate remains stable around its mean value.

## DISCUSSION

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

### 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
^{3}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
^{2}

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

The boundary conditions for equations A1 and A2 are the tubular inflow rate, Q_{T}(0), as calculated from the equations of the glomerular model, and the tubular outlet pressure, which is estimated from experimental measurements (25): (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.

#### Glomerular Model

Because the glomerular capillaries are impermeable to protein, glomerular filtration rate (Table A3), Q_{T}(0), is given by (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.

The filtration process that causes the change in protein concentration is proportional to the sum of local hydrostatic and oncotic pressure differences: (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 segments, which are downstream from the first, each have a variable resistance under the control of both the myogenic mechanism and TGF. The resistance in these segments is given as (A15) where Λ_{j} combines the segment's length and the blood viscosity.

Each of the two arteriolar segments is modeled separately as a set of six ordinary differential equations and associated constituitive relationships. The model draws heavily from the work of Gonzalez-Fernandez and Ermentrout (12).

##### Transport of ions and membrane potential.

The equilibrium distribution of Ca open-channel states, *m*_{∞,j}, as a function of the membrane voltage, *v*_{j}, of the arteriolar segment *j* is given by (A16) where *v*_{1} is the voltage at which half the channels are open, *v*_{2} is a measure of the spread of the distribution, and *j* is the second and third arteriolar segments. This equation is intended to represent a mixture of L- and T- type Ca channels.

For K channels, the distribution is given by (A17) where (A18) Equation 18 provides a Ca-dependent shift in the distribution of K open states with respect to membrane voltage. The time course of the fraction of open K channels, *n*, is given by (A19) where (A20) The time rate of change of the membrane potential is related to membrane currents by (A21) where *C*_{A} is membrane capacitance, and *I*_{L,j}, *I*_{K,j}, and *I*_{Ca,j} are the leak, potassium, and calcium currents of the *j*th arteriolar segment, respectively, and *I*_{IC} is intercellular current. Assuming an ohmic voltage–current relationship: (A22) where *g*_{L}, *g*_{K}, *g*_{Ca}, and *g*_{IC} are the whole-cell membrane conductances for the leak, K, Ca, and intercellular currents, respectively; *v*_{L}, *v*_{K}, and *v*_{Ca} are the corresponding Nernst potentials (Tables A4 and A5), and *g*_{Ca,j,c} = (1 + ζϑ_{j})*g*_{Ca}. ϑ*j* is the input from TGF to the *j*th arteriolar segment, and ζ is the TGF myogenic coupling parameter. The coefficient *g*_{Ca,j,c} therefore represents the membrane Ca conductance, coupled to TGF. The TGF input, ϑ , varies approximately symmetrically about 0.

The rate of change of total Ca concentration in the cytosol is given by (A23) where α_{A} = 1/2β_{A}V_{cell}*F*, β_{A} is the fraction of the total cell volume occupied by the cytosol, V_{cell} is the total cell volume, *F* is Faraday's constant, and *k*_{Ca} is the first-order rate constant for Ca removal from the cytosol. By assuming that cystosolic-free Ca is in equilibrium with various Ca buffers whose total concentration is *B*_{T} and that can be represented by a single equilibrium constant, *K*_{d}, one may arrive at (A24) where (A25)

##### Myosin light-chain phosphorylation.

Activation of smooth muscle is assumed to occur in response to Ca-dependent phosphorylation of myosin light chain. One assumes that the kinetics of phosphorylation are rapid compared with the time scale of the vascular events, and defines ξ as the fraction of myosin light chain sites that are phosphorylated: (A26) Gonzalez-Fernandez and Ermentrout (12) described the relationship between phosphorylated light-chain myosin and cross-bridge formation phenomenologically in terms of a binding distribution. Letting ω_{j} represent the fraction formed of the total possible cross bridges in the *j*th arteriolar segment, they arrived at an expression that we modified for use in the two segment arteriolar model: (A27)

##### Wall stress and contraction velocity.

The contractile mechanism has length *y*, the series elastic component length *u*, and the parallel elastic component length *x* = *u* + *y* (Table A6). The stresses α normal to the surface of a longitudinal slice, *S*_{A}, through the vessel wall and associated with each element are (A28) (A29) and (A30) where *x** = *x*/*x*_{0}, *x*_{0} is a reference length, *s*(*y*) = [*y*_{1}/(*y* + *y*_{2})]^{y4}, and ω_{ref} = ω(Ca_{i,ref})/[ξ_{m} + ξ(Ca_{i,ref})]. The *x* _{1−9}, *u* _{1−3}, and *y*_{1−3} are coefficients used to fit the various expressions to data from the literature.

The force velocity relationship for the contractile mechanism is given by or (A31) where υ_{ref}, *a*_{A}, *b*_{A}, *c*_{A}, and *d*_{A} are constants (Table A7). The hoop forces *f* acting on a longitudinal section *S*_{A} of the blood vessel wall are given by (A32) where Δ*P* is the transmural pressure difference between the vascular wall and the interstitial space, and *A* is the cross-section area of the vascular smooth muscle cells. (A33) (A34) *w*_{e} and *w*_{m} are weighting factors.

The rate of change of the circumference *x* is given by (A35) subject to the geometric compatibility condition (Table A8) (A36)

## APPENDIX B. THE BISPECTRUM

The bispectrum can be estimated using the Fourier transform of the third-moment sequence. To estimate the autoregressive bispectrum of the data set {*y*(1), *y*(2),…, *y*(N)}, the data are first segmented into *K* records of *M* samples each, i.e., *N* = *KM*. Then {*y*_{i}(l);l = 1, 2,…, *M*} are the data samples of the *i*th record.

The third-order moment estimate is (B1) where *i* = 1, 2,…, *K*; *p*_{1} = max(0, −*m*, −n); and *p*_{2} = max(*M* − 1, *M* − 1 − *m*, *M* − 1 − *n*). The *r*^{i}(*m*,*n*) are averaged over all segments: (B2) The bispectrum is estimated from (B3) where *L* < *m* − 1 and *W*(*m*,*n*) is a two-dimensional window. Bispectrum analysis as defined here detects nonlinear interactions, but it provides no information about the coupling strength. To satisfy this need, the bicoherence can be estimated: (B4) where *B*(ω_{1}, ω_{2})is the bispectrum estimate as defined and P(ω) is the power spectrum. A bicoherence value of 1 indicates a strong interaction, a value of 0.5 indicates a borderline interaction, and 0 indicates an absence of interaction.

To understand why the bispectrum detects nonlinear interactions, consider the process: (B5) where λ_{3} = λ_{1} + λ_{2}, λ_{6} = λ_{4} + λ_{5}, φ_{1}, φ_{2}…, φ_{5}, are all independent, uniformly distributed random variables over (0, 2π), and φ_{6} = φ_{4} + φ_{5}. Note that while (λ_{1}, λ_{2}, λ_{3}) and (λ_{4}, λ_{5}, λ_{6}) are at harmonically related positions, only λ_{6} is a result of phase coupling between those at λ_{4} and λ_{5}, and λ_{3} is due to an independent harmonic component. By definition, the bispectrum is obtained by taking a two-dimensional Fourier transform of the third-order cumulant sequence. Thus, if we take the third-order cumulant sequence of the process *x*(*k*), we obtain (B6) Note that in the above equation, only the phase coupled components (λ_{4}, λ_{5}, λ_{6}) appear. Thus the bispectrum is a useful technique for detecting quadratic phase-coupled components.

The quadratic phase coupling or nonlinear interactions can be seen by passing the signal (B7) where A_{1} and A_{2} are constants, through the following quadratic zero-memory nonlinear system (B8) in which *a* and *b* are also constants. The output of the quadratic system {*y*(*t*) = *ax*(*t*) + *bx*^{2}(*t*)} contains consinusoidal terms with the following frequencies and phases: (*f*_{1},φ_{1}), (*f*_{2},φ_{2}), (2*f*_{1}, φ_{2}), (*f*_{1} + *f*_{2},φ_{1} + φ_{2}), and (*f*1 − *f*2, φ_{1} − φ_{2}).

Such a phenomenon that gives rise to the sum or difference of phases with the same relationship to that of the frequencies is called phase and frequency coupling. If only the single sinusoidal term (*f*_{1},φ_{1}) is passed through the quadratic nonlinear system, only the component (2*f*_{1}, 2φ_{1}) will result at the output. This case exhibits only self-phase and self-frequency coupling. Therefore, the presence of (self) frequency and (self) phase coupling indicate the presence of a zero-memory nonlinear system. The presence of frequency coupling indicates that the system is nonlinear, and the presence of phase coupling provides additional information about the time nature of the system.

## Acknowledgments

This work was supported by NIH Grants DK19568 to D. J. Marsh and HL69629 to K. H. Chon, a grant from the EU commission (BioSim, Project 005137) to N.-H. Holstein-Rathlou, and by a grant from the Lundbeck Foundation to O.V. Sosnovtseva. The authors are indebted to Professor G. Bard Ermentrout for his helpful discussions and advice on the implementation of his model of vascular smooth muscle.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

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

- Copyright © 2005 the American Physiological Society