点击显示 收起
【关键词】 myogenic
1Department of Biomedical Engineering, State University of New York at Stony Brook, Stony Brook, New York
2Department of Physiology and Biophysics, University South Florida, Tampa, Florida
3Department of Molecular Pharmacology, Physiology
Biotechnology, Brown University, Providence, Rhode Island
ABSTRACT
We previously showed that nonlinear interactions between the two renal autoregulatory mechanics (tubuloglomerular feedback and the myogenic mechanism) were observed in the stop flow pressure (SFP) and whole kidney blood flow data from Sprague-Dawley rats (SDR) using time-invariant bispectrum analysis (3, 4). No such nonlinear interactions were observed in either SFP or whole kidney blood flow data obtained from spontaneously hypertensive rats (SHR). We speculated that the failure to detect nonlinear interactions in the SHR data may be related to our observation that these interactions were not continuous and therefore had time-varying characteristics. Thus the absence of such nonlinear interactions may be due to an inappropriate time-invariant method being applied to data that are especially time varying in nature. We examine this possibility in this paper by using a time-varying bispectrum approach, which we developed for this purpose. Indeed, we found significant nonlinear interactions in SHR (n = 18 for SFP; n = 12 for whole kidney blood flow). Moreover, the duration of nonlinear coupling is found statistically to be longer (P = 0.001) in SFP data from either SDR or SHR than it is in whole kidney data from either type of rat. We conclude that nonlinear coupling is present at both the single nephron as well as the whole kidney level for SDR and SHR. In addition, SHR data at the whole kidney level exhibit the most transient nonlinear coupling phenomena.
time-varying bispectrum; nonlinear interactions; tubuloglomerular feedback; phase coupling
PHYSIOLOGICAL SYSTEMS OFTEN exhibit complex time-dependent behavior. The beating of the heart, respiration, the organized activity of the nervous system, the 24-h rhythm, the menstrual cycle, and oscillations in the system regulating renal blood flow (9, 16) are examples. Physiological systems consist of components that usually exhibit nonlinear behavior; a list of examples may include simple enzymes following Michaelis-Menten kinetics, membrane transporters, intracellular signal cascades and regulatory networks, neurons, networks of neurons, and gene expression. A question at issue in each case is whether nonlinear interactions occur among such components. 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 (e.g., the Poincare oscillator), or from the nonlinear coupling of two linear systems (11). The reason for the interest in the detection of nonlinear interactions is that they can give rise to a number of system properties, including chaos, synchronization, and frequency modulation (14), which may be physiologically important, and which do not occur in linear systems.
We presented evidence for nonlinear interactions between tubuloglomerular feedback (TGF) and the myogenic mechanism. These regulatory systems combine to produce autoregulation of renal blood flow and glomerular filtration rate (2, 4, 6, 19, 20, 2225, 27). Most of this work has used methods that assume the properties of the signal being analyzed do not vary with time. For a signal with an oscillation, the condition of time invariance (also termed stationarity) means that the frequency of the oscillation does not change with time. We have shown that time invariance is not satisfied in the kidney (4, 19, 20, 22, 27). In this paper, we offer a new method for the analysis of time-varying nonlinear interactions in the renal circulation and apply it to experimental data. The method was developed in the course of our work on the kidney, but the development is quite general, and the method could be applied to other systems.
Because time series can be approximated as the sum of sine and cosine terms, they can be characterized in terms of their frequency and phase. Quadratic frequency coupling and quadratic phase coupling are manifestations of nonlinear interactions (21). Frequency and amplitude modulations are the common examples corresponding to the quadratic frequency coupling and quadratic phase coupling phenomena, respectively. When there is quadratic frequency coupling between the characteristic and the modulation frequency, the power spectrum contains a third peak at a frequency that is either the difference or sum of these two frequencies. With phase coupling, the phase associated with the additional peak is the sum of the two phases associated with the two frequencies. In addition, the self-coupling peak, while physically meaningless, can also occur with both quadratic frequency and phase coupling (21). The power spectrum, being a linear approach, is unable to determine whether a third peak that is either the difference or sum of the two frequencies is due to quadratic frequency, phase coupling, or both. Therefore, nonlinear methods are needed to determine the presence of quadratic frequency and phase coupling. The definition of quadratic phase coupling implies the presence of frequency coupling. However, quadratic frequency coupling can occur without associated quadratic phase coupling. Quadratic frequency and phase coupling are phenomena that arise from the transformation of time series into the frequency domain. They are not descriptions of specific kinds of interactions, which require mechanistic models to identify.
Analytic methods such as Volterra analysis and the bispectrum can be used to test whether harmonically related peaks are the result of nonlinear interactions. Volterra analysis can only be used to determine quadratic frequency coupling, whereas the bispectrum provides evidence of coupling if the conditions of quadratic frequency and phase coupling are both satisfied. For example, for the time-invariant bispectrum, each of the harmonically related peaks will be characterized by statistically independent random phases. The procedure involved in the computation of the bispectrum, which involves statistical averaging of many segments of the bispectrum, will result in a zero value as random phases will cancel out. This phenomenon of a zero bispectrum value is also true for quadratic frequency coupling alone using the time-invariant bispectrum. However, if the third peak generated by the sum or difference of the two peaks is generated by true nonlinear interactions, phase coherence exists and the statistical averaging will not result in a zero bispectrum value (21). Therefore, for a peak to appear in the time-invariant bispectrum, both quadratic frequency and phase coupling must occur. The time-varying bispectrum, because it involves estimation of the bispectrum for each time instant, can only detect quadratic frequency coupling and not phase coupling. As is the case with the time-invariant bispectrum, only the two frequencies responsible for generation of a peak due to the sum or difference of their frequencies (but without an additional phase coupling requirement) will appear in the time-varying bispectrum. A special case occurs in which a self-coupling peak may also appear in a time-varying bispectrum. Summarizing, quadratic phase coupling detection using the time-invariant bispectrum or quadratic frequency coupling in the time-varying bispectrum are both manifestations of nonlinear interactions.
Using Volterra analysis, we previously showed nonlinear interactions between the myogenic mechanism and TGF at the level of the whole kidney blood flow in normotensive rats using randomly forced arterial pressure fluctuations (3). Most recently, using time-invariant bispectral analysis, we detected nonlinear interactions between the two autoregulatory mechanisms in the stop flow tubular pressure recordings obtained from normotensive rats, but not in spontaneously hypertensive rats (SHR) (4). The Volterra analysis method was not used in either that study or the present one because it requires a broadband signal to ensure fidelity of the results, and data from single nephron measurements do not necessarily satisfy this requirement.
With the time-invariant method of bispectral analysis, we were able to detect nonlinear interactions between TGF and the myogenic mechanism in normal Sprague-Dawley rats (SDR), but not in SHR, which have a genetic form of hypertension (4). There is no evidence we know of that either TGF or the myogenic mechanism is inactive in SHR, and our inability to detect nonlinear interactions is more likely due to a shortcoming of the analytic method. Blood flow and renal tubular pressure records from SHR contain irregular, chaotic-like interactions in the TGF fluctuations (24, 25), and it seems likely that these irregular fluctuations increase the time variance of the interactions. Another reason for our inability to detect nonlinear interactions may be that the magnitude of the myogenic oscillation is considerably smaller in SHR than in normotensive SDR (4).
Two studies using time-varying spectral analysis on renal blood flow data revealed nonstationary dynamics of the renal autoregulatory mechanisms (4, 27) with SHR exhibiting greater nonstationarity than SDR. Both the bispectrum and Volterra analyses used to detect nonlinear interactions in the whole kidney blood flow and tubular pressure data are time-invariant methods. We were limited to these time-invariant methods because no reliable time-varying methods were available. It is likely, then, that the time-invariant methods used in the past to detect nonlinear interactions may not be appropriate to detect possible transient interactions, especially in SHR. Thus the goals of the present study were 1) to develop a high- resolution and high-fidelity time-varying bispectral method, 2) to use this new time-varying bispectral method to determine whether there are nonlinear interactions between the two renal autoregulatory mechanisms in SHR at the level of whole kidney blood flow as well as single-nephron tubular pressure data, and 3) to verify that the presence of nonlinear interactions previously observed in SDR data using time-invariant methods is still found when using the new time-varying method.
METHODS
Experimental Methods
All experiments were performed under protocols approved by The Institutional Animal Care and Use Committee at Brown University and The University of South Florida, in accordance with The National Institutes of Health policy on the care and use of laboratory animals.
All experimental data were from previous studies (2, 10). To obtain measurements of proximal tubule stop flow pressure (SFP) in SDR (240300 g) and SHR (12 wk old), the subjects were anesthetized with halothane administered in an oxygen-nitrogen mixture and were artificially ventilated after the administration of a muscle relaxant. The left kidney was exposed. Tubular flow was interrupted with bone wax in a selected proximal tubule. Intratubular hydraulic pressure proximal to the wax block was measured with a micropipette (1- to 3-μm diameter) attached to a servo-nulling pressure circuit. For measurement of whole kidney blood flow, an electromagnetic flow probe was placed around the left renal artery. Arterial pressure was measured in the superior mesenteric artery while broadband fluctuations were induced into the distal aorta (10). Arterial pressure, proximal tubular pressure, and renal blood flow were recorded simultaneously on a TEAC R-61 4-channel cassette data recorder for off-line analysis. The interval of recording ranged from 15 to 20 min. The recorded data were replayed through an electronic low-pass filter with a roll-off frequency of 1.5 Hz, and sampled digitally at 4 Hz for tubular pressure, and at 1 Hz for renal blood flow.
Data Analysis
We further limited the data bandwidth to 0.5 Hz from the experimental sampling rate of 4 Hz for both tubular pressure and the whole kidney blood flow (by downsampling by a factor of 4 after low-pass filtering at 0.5 Hz) to focus more on the frequency bands of interest to autoregulation. We analyzed 15 SFP data records for SDR and 18 for SHR. In addition, 12 whole kidney blood flow recordings were used for each type of rat, SDR and SHR. Each time series, consisting of 7502,400 data points, was subjected to linear trend removal (fitting a first-order least squares-error polynomial to the data and subtracting it from the data record), demeaned (by subtracting out the mean value), and normalized to unit variance (18).
Spectral analysis. Time-invariant power spectra for tubular pressure data were obtained using the modified covariance autoregressive (MODCOVAR) spectrum (18), defined as:
(1)
where w, a constant value, is an estimate of the driving noise variance, T is the sampling interval, and a(n) are the autoregressive (AR) parameter estimates. The model order P, which determines how many of the a(n) terms are necessary to obtain the most essential periodic components in the signal, was determined using a model order criterion recently developed in our laboratory, known as the optimal parameter search (OPS) (17). The MODCOVAR spectrum involves extraction of a few a(n) parameters that characterize essential periodic components in the data, followed by the computation of Eq. 1 to obtain the power spectrum. The MODCOVAR spectrum was used instead of the classic power spectrum because it is less data-length dependent and provides higher frequency resolution (18).
To determine whether nonlinear interactions occur between the myogenic mechanism and TGF for the whole kidney blood flow data and single nephron tubular pressure data, time-varying bispectral analyses were performed. Unlike the time-invariant and time-varying power spectrum, a time-varying bispectrum is a useful tool for detecting and quantifying the presence of time-varying quadratic frequency coupling (QFC), a phenomenon in which two components interact and generate a third component with a frequency equal to the sum or difference of the two components. The goal of this work is to overcome this limitation and to determine whether there are time-varying QFC phenomena in the blood flow and tubular pressure data.
The time-varying bispectrum can be estimated using either Fourier Transform (FT) approaches or model-based approaches such as the AR process. For our analysis, we used a time-varying AR bispectrum as AR approaches have been shown to provide higher-resolution QFC detection than do FT-based methods (21). Details regarding the implementation of the time-varying AR bispectrum are provided in the APPENDIX. It should be noted that the time-invariant bispectrum detects only the phase-coupled frequencies, due to averaging of the third-order cumulant sequences as provided by Eq. A5 in the APPENDIX. With the time-invariant bispectrum, due to this averaging process, any frequency coupling without phase coupling will not appear in the bispectrum plot. However, with the time-varying bispectrum, because of the nature of its design, only frequency coupling can be detected (as it is a time-varying algorithm, averaging cannot be done, thus only frequency coupling can be detected). An additional difference between time-invariant and time-varying bispectrum is that the former is a three-dimensional (f1, f2, and magnitude) plot, whereas the latter involves a four-dimensional (t, f1, f2, and magnitude) plot. The presence of either frequency or phase coupling indicates that the system is nonlinear, and the presence of phase coupling provides additional information about the time nature of the system. The quadratic phase coupling or frequency coupling can be seen by passing the signal
(2)
where A1 and A2 are constants, through the following quadratic zero-memory nonlinear system
(3)
in which a and b are also constants. The output of the quadratic system [y(t) = ax(t) + bx2(t)] contains cosine terms with the following frequencies and phases: (f1, 1), (f2, 2), (2f1, 21), (2f2, 22), (f1 + f2, 1 + 2), and (f1 f2, 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 quadratic phase and frequency coupling. If only one sinusoidal term, (f1, 1), is passed through the quadratic nonlinear system without the linear term in Eq. 3, only the component (2f1, 21) will result at the output. This case exhibits only self-phase and self-frequency coupling. This indicates the presence of a zero-memory nonlinear system.
RESULTS
Simulation Example
With the use of a simulation example, we demonstrate the concept of phase coupling and frequency coupling via the time-invariant bispectrum and time-varying bispectrum, respectively. Consider a signal with 1,024 data points in which the first 512 s of the signal are expressed as:
(5)
and the second 512 s of the signal are expressed as:
(6)
where f1 = 0.1 Hz, f2 = 0.04 Hz, f3 = 0.27 Hz, and f4 = 0.2 Hz and 1, 2, 3, 4, 5, 6, and 7 were randomly selected from a uniform distribution on [0, 2]. Note that while the first half of the signal contains one frequency pair that is both frequency and phase coupled [(f1 + f2) in Eq. 5] and one pair that is only frequency coupled [(f2 + f3) in Eq. 5], the second half only contains one frequency pair, which is both frequency and phase coupled [(f2 + f4) in Eq. 6].
The above expressions as described by Eqs. 5 and 6 represent a signal in which frequencies are time varying, as the frequencies in Eq. 5 are present only during the first 512 s of the data and the frequencies in Eq. 6 are present only in the last 512 s of the data. To illustrate the time-varying characteristics of Eqs. 5 and 6, we computed the smoothed pseudo-Wigner-Ville (SPWV) time-frequency spectrum, shown in Fig. 1B, along with its time-invariant power spectrum counterpart, shown in Fig. 1A. It can be seen that the difference between the power spectrum and the SPWV is that the former provides only the frequencies present in the signal, while the latter method provides concomitant time and frequency information. Although the SPWV provides time-varying characteristics of the frequencies present in the signal, neither the SPWV nor the power spectrum provides either the frequency or phase coupling present in these frequencies.
To detect phase and frequency coupling, the time-invariant bispectrum and time-varying bispectrum were used, and the results are shown in Figs. 2, A and B, respectively. For the time-invariant bispectrum, the result is based on using a model order of 20 in the OPS (16) as described in the APPENDIX. As discussed previously, only the phase-coupling peaks (0.1, 0.04 and 0.2, 0.04 Hz) including self-phase coupling (0.1, 0.1 Hz) should be present in the time-invariant bispectrum plot. Indeed, as shown in Fig. 2A, two prominent peaks due to phase coupling, at 0.1, 0.04 Hz and 0.2, 0.04 Hz, and a smaller magnitude peak at 0.1, 0.1 Hz, due to self-coupling, are identified. Note that the frequency coupling peak at 0.27, 0.04 Hz is not present in the time-invariant bispectrum. This is because the frequency pair is only frequency coupled and not phase coupled. Note that detection of phase coupling with the time-invariant bispectrum requires both frequency and phase coupling. Although the time-invariant bispectrum shows the presence of the correct phase-coupling peaks at 0.1, 0.04 Hz and 0.2, 0.04 Hz, these peaks are inappropriately shown to occur over the entire duration of the simulation. As the simulation example was designed, the phase-coupled peak at (0.1, 0.04) Hz and (0.2, 0.04) Hz should only occurs for the first 512 s and the second 512 s, respectively.
As with the missing time information with the power spectrum, the time-invariant bispectrum also does not show the time during which phase coupling occurs; the (0.1, 0.04) Hz pair exists only in the first half of the signal, whereas the (0.2, 0.04) Hz pair exists in the second half of the signal. The information about when these frequency couplings occur is provided in Fig. 2B using the time-varying bispectrum approach. As noted above, the time-varying bispectrum is presented as a four-dimensional plot in which time is the additional variable. Figure 2B shows when the frequency coupling occurs. For the simulation example, the presence of frequency coupling at (0.1, 0.04) Hz and (0.27, 0.04) Hz for the first 512 s and only the frequency coupling at (0.2, 0.04) Hz for every time instant of the second 512 s is expected. Indeed, Fig. 2B illustrates that we have two frequency couplings, at (0.1, 0.04) Hz and (0.27, 0.04) Hz, shown in "x" and "o" lines, respectively, occurring during the first 512 s. In the last 512 s, we have a single frequency coupling occurring at (0.2, 0.04) Hz, shown in "+" lines. Note that a peak at (0.1, 0.04) Hz is detected even though it is a phase coupling, as the definition of phase coupling refers to the condition where both frequency and phase couplings are satisfied. One important difference between the time-invariant and time-varying bispectra is the inability of the time-varying bispectrum to distinguish peaks that are both frequency and phase coupled; it can only detect frequency coupling. On the other hand, as the time-invariant bispectrum involves averaging of the triple correlation over a number of data records, the peaks not phase coupled are suppressed in the resulting time-invariant bispectrum. However, the time-varying method calculates the instantaneous bispectrum and as such does not involve any averaging over time. Thus frequency coupling alone is required for a peak to appear in the time-varying bispectrum. However, this is not a serious handicap in our analysis of renal autoregulation because phase coupling provides information on which signal component leads or lags. Without phase information, we may still reliably capture the dynamics of nonlinear interactions, albeit without knowledge of the lead/lag between the signal components.
The time-varying bispectrum estimation technique just presented is based on fitting an AR model to the signal. An AR model-based time-varying bispectrum is an extension of the time-invariant AR bispectrum (21), as the coefficients in Eq. A16 are functions of time. As is the case with the time-invariant bispectrum, an AR model-based time-varying bispectrum provides better frequency resolution as well as the capability to be more accurate in the face of noise corruption than non-model-based methods (21). To demonstrate this, we compare our time-varying bispectrum with a non-model-based time-varying bispectrum estimation technique, which relies on a moving-window scheme (11). This method assumes that within a small window, the signal is stationary, and the bispectrum is computed within this defined window. The final bispectrum is then the average of all windowed bispectra.
The simulated signal is based on Eqs. 5 and 6, and we also consider a case in which white noise is added to the clean signal so that the signal-to-noise ratio is 10 dB. All results are averaged based on 100 realizations of Monte Carlo simulations. The signal-to-noise ratio of 10 dB represents a signal in which the variance of the signal is approximately three times greater than the variance of the noise. Figure 3, left, top and bottom, shows the comparison of the AR time-varying bispectrum vs. the windowed time-varying bispectrum for a clean signal at time instant t = 1,024 s. As per our simulation, at t = 1,024 s we should observe a clear peak in the bispectrum at (0.2, 0.04) Hz. This peak is clearly seen in both the AR (our method) and windowed (method described in Ref. 11) bispectrum estimates, but the resolution is much greater and there are fewer spurious peaks with the AR time-varying bispectrum. The same simulation but with 10 dB of white noise corruption are shown in Fig. 3, right, top and bottom. Clearly, the AR method provides a higher resolution bispectrum estimate and is less affected by the noise corruption than the windowed bispectrum estimate. Although not shown, the AR bispectrum method consistently provides accurate results up to 0 dB of additive white noise and better results at these noise levels than does the windowed bispectrum approach.
To better explain the concept of a second-order nonlinear (quadratic) frequency coupling, we designed the following simulation in which signal y is generated as follows:
(7)
where s = sin(2f1) + sin(2f2) in which f1 = 0.1 Hz and f2 = 0.22 Hz. Using simple trigonometric identities, we may expand Eq. 7 to:
Therefore, the signal y contains components which oscillate at frequencies f1, f2, 2f1, (f1 f2), and (f1 + f2). Indeed, the time-frequency spectrum of y in Fig. 4A depicts these frequencies precisely and that these frequencies occur over all time, hence correctly suggesting that the signal is time invariant. Note that signal y was generated by quadratic frequency coupling between signal components oscillating at frequencies f1 and f2. Neither the time-frequency spectrum plot nor power spectrum (not shown) provides information on the existence of this quadratic frequency coupling. On the other hand, the time-varying bispectrum at time instant t = 97 s clearly shows the existence of this quadratic frequency coupling, shown in Fig. 4B. This peak is a direct result of the existence of quadratic frequency coupling at the pair of frequencies at f1 = 0.1 Hz and f2 = 0.22 Hz.
Application to Experimental Data
We applied the time-varying bispectrum to detect QFC in the SFP data recordings from 15 SDR and 18 SHR as well as 12 each of SDR and SHR whole kidney blood flow recordings. Figure 5A, top and bottom, shows a typical tracing of SDR SFP data and its power spectrum, respectively. The power spectrum was computed using the MODCOVAR technique, as described in Eq. 1 (18). For this SFP record, shown in Fig. 5A, top, we observe dominant spectral peaks at f1 = 20, f2 = 80, f3 = 100, and f4 = 120 mHz, shown in Fig. 5A, bottom. We consider peaks at 20 and 100 mHz to be the characteristic frequencies of the TGF and myogenic mechanisms as these are in agreement with the previously reported frequency ranges of the TGF (2050 mHz) and myogenic (100200 mHz) mechanisms (810). The time-varying characteristics of these peaks are illustrated in Fig. 5B. To investigate whether peaks observed in the power spectral plot are spontaneously excited normal modes, coupled modes, or due to harmonic generation, and to determine whether these scenarios are time dependent, the time-varying bispectrum was estimated. A time-varying bispectral analysis is well suited to test whether the frequency modes observed in the power spectrum are spontaneously excited normal modes of the system or whether they are coupled modes, and whether they are time dependent. Furthermore, we can use the time-varying bispectrum to determine whether a frequency mode was due to harmonic generation (e.g., fx= 2fy) that is time dependent. Figure 5, C and D, shows the time-varying bispectral magnitude and its contour plots, respectively, obtained at time instant t = 322 s. The time-varying bispectral plots (Fig. 5, C and D) exhibit a dominant peak at the frequency pair (100, 20) mHz. Thus we conclude that the primary modes, f1 = 20 mHz and f3 = 100 mHz, interact and generate modes at f4 = 120 mHz and f2 = 80 mHz (f4 = f1 + f3 and f2 = f3 f1). Note that these coupled modes correspond to nonlinear interactions between the myogenic and TGF mechanisms. Furthermore, a self-coupled peak is also observed at the TGF frequency.
Figure 6A provides typical SFP data from a SHR (top) and the resultant power spectrum (bottom). In addition, the corresponding time-varying spectrum is provided in Fig. 6B. There are greater time-varying characteristics than for the SDR, shown in Fig. 5B. We observe in Fig. 6A dominant spectral peaks at f1 = 50, f2 = 70, f3 = 100, f4 = 120, f5 = 140, and f6 = 170 mHz. We consider the peaks at 50 and 100 mHz to be the characteristic frequencies of the TGF and myogenic mechanisms as these are in agreement with the previously reported values (810). Time-varying bispectral analysis, justified by the time-varying characteristics of these peaks, was conducted to determine whether the peaks are due to nonlinear interactions between the TGF and myogenic mechanisms. The time-varying bispectral plots (Fig. 6, C and D) are shown for the time instant t = 1,387 s. The plots reveal a dominant peak at the frequency pair (120, 50) mHz. Thus we can conclude that the primary modes f1 = 50 mHz and f4 = 120 mHz interact and generate modes at f2 = f4 f1 = 70 mHz and f6 = f4 + f1 = 170 mHz. The modes at f3 = 100 mHz and f5 = 140 mHz are harmonics of f1 = 50 mHz and f2 = 70 mHz, respectively (f3 = 2f1 and f5 = 2f2).
For each time series, the time-varying bispectrum was estimated at each sampled data point. The model order, P, as denoted in Eq. A10 in the APPENDIX, was estimated using the OPS, which we previously demonstrated to be one of the most accurate methods for finding model orders to date (16). The time-varying bispectrum, estimated via Eq. A17, was then analyzed to detect those frequency pairs where QFC was statistically significant. To test whether the nonlinear phase coupling in the SFP time series is significant, the time-varying bispectral power value from Eq. A17 was calculated at every frequency pair and was used as an index to measure the degree of coupling. To ensure that the bispectral power does indeed provide statistical evidence of frequency coupling, we employed the white noise technique described in Ref. 5.
A statistical test was applied to the null hypothesis that there is no coupling. The threshold for this test was computed as follows: 100 realizations of a Gaussian white noise process with zero mean and unit variance were generated. An AR model (with the same model order as that of the SFP and the whole kidney blood flow data) was fit to this process. A Kalman filter (12) was used to estimate the mean time-varying AR coefficients. The time-varying bispectrum was estimated at each sampled data point, and at every instant the mean bispectral power was found by averaging over all frequencies. At a statistical significance level (0.05 in our analysis), the (1 )th percentile of the bispectral power was taken to be the threshold. We thus obtained a time-varying bispectral threshold. This was compared with the bispectrum obtained from the stop pressure and whole kidney flow data to determine whether coupling was significant.
Because it is not possible to document the bispectral plot and QFC details for every time instant, we present our results in a summarized format. For each of the four data categories (SFP SDR, SFP SHR, whole kidney SDR, and whole kidney SHR), we provide summarized values of the following: 1) fraction of time when no significant QFC occurs, which is a measure of how often coupling does not occur; and 2) length of the longest coupling duration and mean length of continuous QFC (defined below).
We define duration of the QFC as those consecutive time instants where QFC is present. In addition, we impose an additional constraint on the consecutive duration of the QFC such that the difference between the myogenic frequencies at the current time instant and the previous time instant does not exceed the TGF frequency at the current time instant. This constraint is imposed to ensure that variations in the coupling frequency are not caused by changes in TGF alone and that they indeed reflect actual presence of the QFC. Figure 7 depicts the detection of coupling for a typical data set from each data category. The value 1 represents time instants where coupling is present while the value 0 represents the absence of coupling. It can be seen that QFC was detected a majority of the time in single nephron data for both SDR and SHR. Compared with SFP data, whole kidney data from both SDR and SHR exhibited less QFC.
To quantify the extent of the QFC in a time series, we computed the amount of time that QFC is absent. This number, expressed as a fraction of the total time duration, is a measure of the extent of the QFC in the signal. Histograms showing the distribution of this fraction for all experimental subjects for each type of data set are provided in Fig. 8. We can see that in single-nephron SDR data (SFP SDR), this fraction was 0 in 12 out of 15 data sets. This means that in these 12 data sets, QFC was detected at each sampled data point. In SFP SHR, this fraction was 0 in 9 of 18 data sets. There is no statistical difference in coupling duration between SFP SHR vs. SFP SDR, and WK SDR vs. WK SHR. However, whole kidney SDR data (WK SDR) show less coupling duration than SFP SDR and SFP SHR (chi-square, P < 0.02). Likewise, WK SHR show less coupling duration than SFP SDR and SFP SHR (2, P < 0.03). Figure 9, A and B, plots the maximum length and mean length of the chains of continuous QFC, respectively, over different categories. It can be seen that maximum and mean chain lengths are significantly longer in SFP data than in whole kidney data for both SDR and SHR. These results indicate that the length of the chains in SFP data obtained from SDR was statistically the greatest among all data categories tested. Furthermore, even though chain length in SFP data from SHR was less than in SFP data from SDR, it was still significantly greater than chain lengths in whole kidney data from either SDR or SHR. There was no significant difference between WK-SDR and WK-SHR for either maximum or mean chain lengths.
DISCUSSION
In this study, we developed and implemented a novel approach to detecting nonlinear interactions between TGF and the myogenic mechanism without the implicit assumption of stationarity of the signal. This was achieved by sequential calculations of Eqs. A16 and A17, which define the time-varying bispectrum. The validity of this approach was confirmed by using a simulated data set.
The main findings are that nonlinear couplings between TGF and the myogenic mechanism are found in both SFP and whole kidney blood flow data obtained from SHR and SDR using the time-varying bispectrum method, which we have developed for this application. No such coupling was obtained with this data using a time-invariant bispectrum (4) and time-invariant Volterra kernel approaches (3). Moreover, it was found that SFP data from both SDR and SHR exhibited more persistent coupling than did the whole kidney blood flow data obtained from either SDR or SHR. For example, mean coupling length was significantly longer for the SFP than the whole kidney for both SDR and SHR data.
Before the current study, we were limited to utilizing time-invariant bispectral methods because no accurate time-varying bispectral methods were available. In our recent study, we raised the possibility that the absence of nonlinear interactions in the SFP data from SHR may be due to time-varying characteristics, which the time-invariant bispectrum is not designed to capture (4). This motivated us to develop a time-varying version of the bispectrum to test this conjecture. The newly developed time-varying bispectrum is based on an AR model (with parameters estimated via Kalman filtering), which provides one of the highest time and frequency resolutions. This property ensures that even subtle time-varying dynamics (nonlinear coupling) can be detected using this time-varying bispectrum approach.
The time-invariant bispectrum is formulated based on the assumption that nonlinear couplings are stationary over all time segments. Thus it essentially provides the mean nonlinear couplings over the whole duration of the data. However, most physiological systems, including renal autoregulation, involve feedback loops, which are ultimately time dependent. In addition, under certain conditions physiological systems may only exhibit short-duration time-varying characteristics, which are then subsequently masked by time-invariant approaches.
It should be noted that the OPS (17), although an optimal model order estimation method, could produce some ambiguity in the choice of the optimum model order. In certain cases, an over- or underdetermined model order may lead to introduction of nonessential additional peaks or absence of essential peaks in the bispectrum plot, respectively. To ascertain whether results are highly dependent on the choice of the selected model order, we also performed bispectral analysis with model orders that are both lower and greater than the chosen optimal model order (e.g., optimal ± 3). The time-varying bispectrum results did not vary within these confined bounds of the optimal model order, suggesting that the optimal choice of the model order is not critical to the outcome.
Spectral methods (including the bispectrum) of data analysis are especially well suited for renal tubular data because both TGF and the myogenic mechanism exhibit self-sustained oscillations(3, 810, 16, 26). The oscillations of the TGF system are the larger and slower of the two and can easily be detected as an oscillation in tubular pressure. The oscillations of the myogenic mechanism are smaller and more frequent and are not usually seen in tubular pressure recordings because the tubular compliance, acting as low-pass mechanical filter, removes them. Special analytic techniques are required to detect the myogenic oscillations in records of tubular pressure (22, 27).
The SFP technique was designed to study TGF in the open loop mode. By blocking the tubular outflow path, the technique causes tubular pressure to rise about threefold until it approximately equals the sum of glomerular capillary pressure and plasma colloid oncotic pressure. Tubular obstruction increases the magnitude of the oscillations in tubular pressure caused by the myogenic mechanism because the inhibitory effect of TGF is removed (19, 20), because the rise in tubular pressure distends the tubule (D. J. Marsh, personal observation), which should reduce the attenuation of pressure fluctuations caused by tubular compliance, and because a short, stiff, obstructed tubule has different resonance and mechanical filter properties than a long, relaxed, nonobstructed tubule.
Despite tubular blockade, a strong TGF signal persists in the tubular pressure. We described the phenomenon of nephron-to-nephron coupling via signal conduction over the vascular wall (2, 6, 12, 23, 25); gap junctions formed by connexins provide the structural basis for this phenomenon (1). The coupling is stronger in SHR than in SDR, and synchronization of the chaotic tubular pressure signals persists in SHR (2, 24, 25). Thus vascular signal conduction from nearby unobstructed nephrons is the most likely cause of the persistence of TGF signals in tubular pressure records from SFP experiments.
The results presented in this paper show more significant nonlinear interactions between TGF and the myogenic mechanism in measurements made from single nephrons than from whole kidney blood flow. Nephron-to-nephron communication promotes synchronization over nephrons supplied by a single cortical radial artery, but not among nephrons supplied by different cortical radial arteries (2, 6, 23, 25). As shown in this and other papers, there is frequency coupling between TGF and the myogenic mechanism in single nephrons (19, 20, 22). In a single kidney, because there are nephrons whose TGF oscillations occur over a range of frequencies, so too will there be a range of frequencies over which the myogenic oscillation occurs, and the result would be reduced opportunities to observe the interaction between the two mechanisms in measurements made in whole kidney blood flow.
The results of this study show the persistence of nonlinear interactions between TGF and the myogenic mechanism in SHR, when the tubular pressure signal from TGF becomes irregular (7, 24, 25); a previous study failed to detect them (4). There is no evidence that the myogenic mechanism is inhibited in hypertension, and one would expect that these interactions should persist, as this study confirms. Our present results show the importance of the new analytic methodology we presented and support additional efforts to explain the interaction of the two regulatory mechanisms in hypertension.
In summary, a novel time-varying bispectrum approach was developed to specifically account for nonstationarity, which in the past may have precluded the detection of nonlinear interactions between the TGF and myogenic mechanisms in tubular pressure recordings as well as whole kidney data from SHR. Indeed, by accounting for nonstationarity in the data via the time-varying bispectrum, we found nonlinear interactions between the TGF and myogenic mechanisms in tubular pressure data from SDR and SHR as well as whole kidney data also from both SDR and SHR.
APPENDIX
To estimate the AR bispectrum, let {y(1),y(2),...,y(N)} be the given tubular pressure data set. The AR process of the given tubular pressure data can be expressed as:
(A1)
where e(n) are independent and identically distributed with E[e(n) = 0, E[e3(n)] = 0 and y(m) is independent of e(n) for m (A2) where R(j,k) is the third-order correlation function, and (j,k) is the two-dimensional impulse response function. The impulse response function is the system's transfer function or admittance gain represented in the time domain. Equation A2 can be written in matrix form as Ra = b where, (A3) Since the true third-order correlation functions as defined above are not known, they are estimated by the following procedure: 1) form the third-order moment estimate, i.e., segment the data into K records of M samples each, i.e., N = KM. Let {yi(l); l = 1, 2,...,M} be the data samples of the i-th record; 2) form the third-order moment estimate: (A4) where i = 1, 2,..., K; p1 = max(0, -m, -n); and p2 = max(M-1, M-1-m, M-1-n). The average ri(m,n) over all segments; (A5) Solve the following matrix equation (A6) with the same matrix and vector notation shown in Eq. A3. The desired AR parameters, a, in Eq. A6 are obtained using the least-squares method. Once the AR parameters, a, are estimated using Eq. A6, the next goal is to compute the bispectrum estimate as (A7) where (A8) It should be noted that the following symmetry exists for the bispectrum: Furthermore, B(1,2) is periodic in 1 and 2 with period 2. Thus knowledge of the bispectrum in the region 2 0, 1 > 2, 1 + 2 is sufficient for a complete description of the bispectrum (21). Estimation of Time-Varying Coefficients Using the Kalman Filter The method explained above can be used to obtain the time-invariant bispectrum of stationary signals. However, to estimate the time-varying bispectrum of nonstationary signals, we need to estimate the time-varying coefficients of the AR process. Consider a nonstationary AR process of order p (A9) where ai(k) are the linear time-varying coefficients to be estimated, y(k) is the measured signal, and v(k) is additive Gaussian white noise. The coefficients of this system can be estimated using a Kalman filter (13). The Kalman filter is used to estimate the state x of a discrete time process that is governed by the linear stochastic equation: (A10) with a measurement y that is given by (A11) where w(k) and v(k) represent the process and measurement noise, respectively. In our application, the state x is the time-varying system coefficient we wish to estimate and the measurement y(k) is the available signal. The matrix B relates the state at the previous time step k 1 to the state at the current step k, while the matrix C in the measurement equation relates the state to the measurement y(k). In our application, the matrix B was defined as B = 1 and the matrix C contains the data lag terms such that C = [y(k 1) y(k 2)... y(k p)] where the model order p was defined in Eq. A9. In relating Eqs. A10 and A11 to Eq. A9, these Eqs. A10 and A11 can be written as (A12) where a(k) now directly refers to the linear AR coefficients in Eq. A9. The task is then to use a Kalman filter to estimate these a(k) coefficients. The Kalman filter involves finding the linear mean square estimator given the measurements y(k). The noise processes are assumed to be independent and white with normal probability distributions (A13) where the first and last arguments in the right side of Eq. A13 represent zero-mean and variances of Q and R, respectively. The Kalman filter recursively uses a set of "prediction" and "correction" equations to obtain the state estimates. The Kalman filter we utilize is based on the minimum mean-squared error principle, known as the innovation approach (13). The prediction equations use the value of the state at the previous time step to project onto the current time step and obtain an a priori estimate of the current state. The prediction equations are: (A14) where P is the prediction error covariance, the superscript "" indicates a priori estimates and the subscript denotes the time instant. The prediction error is defined as and the prediction error covariance as Once the state and error covariance have been projected onto the current time step, the correction equations are used to obtain a more accurate estimate of the current state by utilizing the value of the measurement at the current instant. The correction equations are: (A15) where Kk is known as the Kalman gain and is defined as that gain factor that minimizes the a posteriori error covariance, and R is the covariance matrix of the measurement noise v(k). The last two equations of (A15) are known as the measurement update (13). By recursively using the prediction and correction equations, we can estimate the AR coefficients, ai(k) at each sampled data point. The instantaneous system transfer function is then given by (A16) The time-varying bispectrum can be obtained by (A17) By using the Kalman filter, the linear time-varying coefficients ai(k) and a linear filter H(k,f) in Eq. A16 can be obtained at each sampled data point. A time-varying bispectrum can be estimated at each sampled data point by the triple product of linear filters shown in Eq. A17. The time-varying bispectrum in Eq. A17 is nonlinear because it is the result of this triple product of linear time-varying filters. Thus it should be noted that the extended Kalman filter, which is primarily used for nonlinear systems, does not need to be employed as the time-varying coefficients in Eq. A16 are linear and a standard linear Kalman filter can be used to estimate these coefficients. For further details on the implementation of the Kalman filter, see reference (13). GRANTS This work was supported in part by National Institutes of Health Grants HL-69629 and EB-3508. 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. REFERENCES Arensbak B, Mikkelsen HB, Gustafsson F, Christensen T, and Holstein-Rathlou NH. Expression of connexin 37, 40, and 43 mRNA and protein in renal preglomerular arterioles. Histochem Cell Biol 115: 479487, 2001. Chen YM, Yip KP, Marsh DJ, and Holstein-Rathlou NH. Magnitude of TGF-initiated nephron-nephron interactions is increased in SHR. Am J Physiol Renal Fluid Electrolyte Physiol 269: F198F204, 1995. Chon KH, Chen YM, Marmarelis VZ, Marsh DJ, and Holstein-Rathlou NH. Detection of interactions between myogenic and TGF mechanisms using nonlinear analysis. Am J Physiol Renal Fluid Electrolyte Physiol 267: F160F173, 1994. Chon KH, Raghavan R, Chen YM, Marsh DJ, and Yip KP. Interactions of TGF-dependent and myogenic oscillations in tubular pressure. Am J Physiol Renal Physiol 288: F298F307, 2005. Faes L, Nollo G, and Antolini R. Experimental approach for testing the uncoupling between cardiovascular variability series. Med Biol Eng Comput 40: 565570, 2002. Holstein-Rathlou NH. Synchronization of proximal intratubular pressure oscillations: evidence for interaction between nephrons. Pflügers Arch 408: 438443, 1987. Holstein-Rathlou NH and Leyssac PP. Differences in tubuloglomerular feedbackoscillatory activity between spontaneously hypertensive and Wistar-Kyoto rats. J Hypertens Suppl 3, Suppl 3: S343S345, 1985. Holstein-Rathlou NH and Leyssac PP. TGF-mediated oscillations in the proximal intratubular pressure: differences between spontaneously hypertensive rats and Wistar-Kyoto rats. Acta Physiol Scand 126: 333339, 1986. Holstein-Rathlou NH and Marsh DJ. Oscillations of tubular pressure, flow, and distal chloride concentration in rats. Am J Physiol Renal Fluid Electrolyte Physiol 256: F1007F1014, 1989. Holstein-Rathlou NH, Wagner AJ, and Marsh DJ. Dynamics of renal blood flow autoregulation in rats. Kidney Int Suppl 32: S98S101, 1991. Jamsek J, Stefanovska A, Mcclintock PVE, and Khovanov IA. Time-phase bispectral analysis. Physiol Rev 68: 016201, 2003. Kallskog O and Marsh DJ. TGF-initiated vascular interactions between adjacent nephrons in the rat kidney. Am J Physiol Renal Fluid Electrolyte Physiol 259: F60F64, 1990. Kalman RE. A new approach to linear filtering and prediction problems. Asme J Basic Eng 82: 3545, 1960. Kantz H and Schreiber T. Nonlinear Time Series Analysis. Cambridge, UK: Cambridge Univ. Press, 1997. Kim YC and Powers EJ. Digital bispectral analysis of self-excited fluctuation spectra. Phys Fluids 21: 14521453, 1978. Leyssac PP and Baumbach L. An oscillating intratubular pressure response to alterations in Henle loop flow in the rat kidney. Acta Physiol Scand 117: 415419, 1983. Lu S, Ju K, and Chon KH. A new algorithm for linear and nonlinear arma model parameter estimation using affine geometry. IEEE Trans Biomed Eng 48: 11161124, 2001. Marple SL. Digital Spectral Analysis With Applications. Englewood Cliffs, NJ: Prentice-Hall, 1987. Marsh DJ, Sosnovtseva OV, Chon KH, and Holstein-Rathlou NH. Nonlinear interactions in renal blood flow regulation. Am J Physiol Regul Integr Comp Physiol 288: R1143R1159, 2005. Marsh DJ, Sosnovtseva OV, Pavlov AN, Yip KP, and Holstein-Rathlou NH. Frequency encoding in renal blood flow regulation. Am J Physiol Regul Integr Comp Physiol 288: R1160R1167, 2005. Nikias CP and Petropulu AP. Higher-Order Spectra Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1993. Sosnovtseva OV, Pavlov AN, Mosekilde E, and Holstein-Rathlou NH. Bimodal oscillations in nephron autoregulation. Phys Rev E Stat Nonlin Soft Matter Phys 66: 061909, 2002. Wagner AJ, Holstein-Rathlou NH, and Marsh DJ. Internephron coupling by conducted vasomotor responses in normotensive and spontaneously hypertensive rats. Am J Physiol Renal Physiol 272: F372F379, 1997. Yip KP, Holstein-Rathlou NH, and Marsh DJ. Chaos in blood flow control in genetic and renovascular hypertensive rats. Am J Physiol Renal Fluid Electrolyte Physiol 261: F400F408, 1991. Yip KP, Holstein-Rathlou NH, and Marsh DJ. Dynamics of TGF-initiated nephron-nephron interactions in normotensive rats and SHR. Am J Physiol Renal Fluid Electrolyte Physiol 262: F980F988, 1992. Yip KP, Holstein-Rathlou NH, and Marsh DJ. Mechanisms of temporal variation in single-nephron blood flow in rats. Am J Physiol Renal Fluid Electrolyte Physiol 264: F427F434, 1993. Zou R, Cupples WA, Yip KP, Holstein-Rathlou NH, and Chon KH. Time-varying properties of renal autoregulatory mechanisms. IEEE Trans Biomed Eng 49: 11121120, 2002.