| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
1Department of Mathematics, Duke University, Durham, North Carolina; and 2Department of Physiology and Biophysics, State University of New York, Stony Brook, New York
Submitted 3 February 2005 ; accepted in final form 3 October 2005
 |
ABSTRACT |
|---|
 TOP ABSTRACT MATHEMATICAL MODELRESULTSDISCUSSIONAPPENDIXGRANTSREFERENCES |
|---|
renal hemodynamics; hypertension; mathematical model; nonlinear dynamical system
30 s. Figure 1, A1 and A2, illustrates their findings. Experiments have shown that these regular oscillations are TGF mediated (14, 17, 31). TGF, a negative feedback loop with time delays, has a resonant frequency: mathematical models have indicated that if the feedback is sufficiently strong (i.e., the feedback has sufficiently large gain magnitude), then the stable behavior of the TGF system is a self-sustaining regular oscillation and not a time-independent steady state (SS1) (18, 26, 29). Such a regular oscillation is called a limit-cycle oscillation (LCO). When nephron parameters change, and thereby change TGF loop delays and gain, a bifurcation may occur; i.e., there may be a parameter-dependent change from one type of stable solution (e.g., a SS) to another type of stable solution (e.g., a LCO).
|
We have recently developed a systematic method for investigating, in the context of our model of the TGF system, the impact of NaCl backleak from the interstitium into the TAL. In the absence of backleak, NaCl concentration in the model TAL has been shown to depend on TAL transit time only (28, 39), but the addition of backleak introduces a significant degree of spatial dependence (26, 27). Our investigation, described herein, revealed unexpectedly rich dynamics in a parameter regime that appears to correspond to that reported in SHR. This finding motivated us to reexamine the following question: to what extent can the irregular TGF oscillations in SHR be explained by the intrinsic dynamic characteristics of the TGF system in a single nephron or in a small number of coupled nephrons? Our strategy is to identify unusual behaviors in SHR and key features of the complex power spectra that are computed from the irregular oscillations. Then, we will use our mathematical model to identify factors that may contribute to that spectral complexity.
Many of the key features of TGF-mediated oscillations in SHR, which we seek to explain, are illustrated in Fig. 1B2. This figure, a power spectrum, shows the frequencies, and relative strengths, of the sinusoidal components that, when summed together, make up the long-time record of the waveshape shown in Fig. 1B1. In the spectrum, the numbers above the peaks show the frequencies of the sinusoidal components (in mHz). Careful analysis of this spectrum suggests that elements of order reside amidst the appearance of chaos. First, several strong peaks are found in the range (25–50 mHz) that is typical of the fundamental frequency of TGF-mediated oscillations in normotensive rats (14, 16, 17). Second, some of the peaks appear to be closely spaced doublets or triplets; candidates are seen at 13.5 and 16 mHz; 18.5 and 22 mHz; 27 and 30 mHz; 30 and 33 mHz; 30, 33, and 36.5 mHz; 40 and 44 mHz; 49, 51.5, and 54 mHz; and 72 and 76 mHz. Third, some peak frequencies may be harmonics of lower fundamental frequencies (i.e., higher frequencies that are integer multiples of a lower frequency). Examples of frequencies that appear to be in harmonic relationship are 13.5, 27, 40, 54, and 67 mHz; 16, 33, 49, and 65 mHz; 18.5, 36.5, and 54 mHz; 22, 44, and 65/67 mHz; 33 and 65/67 mHz; and 36.5 and 72 mHz. Another published SHR power spectrum (Fig. 1E in Ref. 52) shows two major peaks that appear to be in harmonic relationship: one apparent fundamental oscillation at 23 mHz (somewhat lower than typical in control rats) and a second oscillation at 47 mHz (somewhat higher than typical in control rats). Fourth, peaks in the SHR power spectrum are often broad compared with those in spectra from normotensive rats. Examples in Fig. 1B2 include peaks in the following ranges: 11.5–17.5 mHz; 17–26 mHz; 28–38 mHz; 38–45 mHz; and 45–56 mHz. Finally, more spectral power is found at very low frequencies (less than 20 mHz) in SHR than in the Wistar-Kyoto (WKY) rat controls (compare Fig. 1, A2 and B2; see also Fig. 1E in Ref. 52 and 
Fig. 8 in Ref. 54). Fluctuations in this frequency range are likely to arise from sources extrinsic to the TGF system, such as shifts in the neurohumoral systems that regulate renal vascular and tubular function.
|
|
|
|
|
|
|
In this study, we examine four factors that may contribute to spectral complexity. The first is the intrinsic complex dynamics of our model TGF system in a parameter regime of high gain consistent with measurements in SHR (11, 32, 46). In this region, we will show that TGF exhibits multistability, in that perturbations can produce LCO at one of two frequencies, a frequency of f1 and a second frequency of f2 that is approximately equal to 2f1 (however, f2 is not a harmonic of f1). Moreover, parameter changes can produce a bifurcation directly from a SS to the LCO of frequency f2. The second factor is alterations in TGF system time delays, which are determined by nephron flow rates, nephron dimensions, and the time required for propagation of the TGF response through the juxtaglomerular apparatus (JGA). A previous analysis of our model predicted that TGF time delays are important determinants of TGF dynamics (26). In this study, we present a more extensive bifurcation analysis in which steady-state TAL flow rate (or, alternatively, TAL volume) is a bifurcation parameter. The results suggest that elongated time delays, and especially TAL transit time, can introduce substantial spectral peaks below the 30-mHz TGF fundamental seen in normotensive rats (see below, Fig. 5B).
The third factor is vascular coupling between nephrons. Experimental studies have shown that nephrons sharing a common feed artery are functionally coupled, in that TGF responses in one nephron will influence the vascular resistance of closely coupled nephrons, leading to synchronous oscillations (7, 14, 22, 53). Finally, our simulations show that gradual shifts in system time delays within a bifurcation region lead to broadened spectral peaks, whereas abrupt changes can produce closely spaced doublets, multiple peaks, and frequency-shifted harmonic series.
Taken together, the results of our study provide a possible explanation for irregular TGF-mediated oscillations in SHR that does not require the postulation of deterministic chaos. Rather, the irregular oscillations in SHR and the associated complex power spectra may be consequences of intrinsic multistability of the TGF system in the parameter regime of the SHR, switching between oscillatory modes, strong nephron coupling between a few nephrons, and temporal variation in key parameters.
|
MATHEMATICAL MODEL |
|---|
|
TOPABSTRACTMATHEMATICAL MODELRESULTSDISCUSSIONAPPENDIXGRANTSREFERENCES |
|---|
Independent Variables
Specified Functions
(t)
Dependent Variables
(t)
Parameters
r2 (µm2)
r (µm)
Q
Q/2Qop (dimensionless)
i
c
i
ipi + /2, effective delay interval at the JGA (s)
pi
ij
Formulation
The mathematical model, which is similar to that previously described (26, 29), represents the TGF system in a short-looped rat nephron. The TAL is modeled as a rigid flow tube, and the actions of the afferent arteriole (AA), glomerulus, proximal tubule, and descending limb are represented by phenomenological relations. The only solute represented in the model TAL is the chloride ion (Cl–), which is believed to be the principal signaling agent for TGF activation along the macula densa (MD) (41). The model has been expanded to represent coupling, through the preglomerular vasculature, with the TGF systems of other model nephrons (each labeled by an index i). Also, the model now includes a parameter i, for each model instance, that scales TAL cross-sectional area. Model variables and parameters, with their dimensional units as commonly reported, are given in the GLOSSARY. A schematic diagram for the ith model TGF system is given by Fig. 2. The spatial variable x, which designates position along the TAL, is oriented so that it extends from the TAL entrance (x = 0), through the outer medulla, and through the cortex to the site of the MD (x = L).
In the ith model the TGF system, TAL tubular fluid chloride concentration Ci(x,t), TAL tubular fluid flow 
, coupling with other nephrons, and a TGF signal delay at the JGA are represented by the following equations
 | (1) |
 | (2) |
 | (3) |
 | (4) |
 | (5) |
Equation 2 gives i(t) as the sum of a nephron's own TGF-mediated TAL tubular flow response Fi(t) and a term from the other model TGF systems (as in Ref. 38); each summand in that term is a weighted difference between another system's TGF response Fj(t) and the time-independent steady-state flow Q. The weight, or coupling coefficient, ij is a nonnegative constant parameter that scales the influence of the jth nephron on the ith nephron.
Equation 3, the feedback response for an individual nephron, is an empirical relationship that has been well established by steady-state experiments (41). The constant Cop is the time-independent steady-state TAL tubular fluid Cl– concentration at the MD when 
i(t) assumes a steady-state value of Q;C(t) represents the effective TAL tubular fluid Cl– concentration after a discrete and distributed time delay described by Eq. 4. The final term in Eq. 3, 
i(t), represents transient perturbations that were introduced to test the stability of solutions.
Dynamic experimental results (5) indicate that a change in MD concentration does not significantly affect AA muscle tension until after a discrete (or "pure") delay time 
and that the subsequent effect is temporally distributed so that a full response requires additional time, with the greatest weight in a time interval [t––, t–], where is a second delay parameter. The convolution integral in Eq. 4 describes the distributed delay signal received by the AA at time t; CMDi is an "effective" chloride concentration at the MD that serves to transmit that delay to the TGF response represented in Eq. 3. Equation 5 provides a weighting for a sigmoidal distribution of the distributed delay; a sigmoidal distribution is consistent with extant experimental data (Fig. 6 in Ref. 5).
For every choice of parameters within the physiological regime, the model equations have a time-independent steady-state solution; however, this solution may be stable (hence, a SS solution), or unstable, in which case an oscillatory solution is stable. The steady-state Cl– profile along the TAL, denoted S(x), can be found by setting 
to the steady-state TGF operating point flow (Fop, which equals Q) and solving the time-independent form of Eq. 1 obtained when the term containing the time derivative is set to zero; the resulting equation is
 | (6) |
i, neither does S(x); therefore S(x) is the same for each model TGF system i.
The dimensionless parameter i, which scales TAL cross-sectional area, has an important alternative interpretation: one can consider i to be a simultaneous rescaling of the steady-state flow rate and of the transport parameters Vmax and P. Thus one can rewrite Eq. 1 as
 | (7) |
(t)/i, Vmax/i, and P/i, respectively. Notwithstanding the reinterpretation of i, Eqs. 7 and 1 have exactly the same solution.The equation for the steady-state solution corresponding to Eq. 7 is
 | (8) |
i's in Eq. 8 cancel, Eqs. 8 and 6 are exactly the same equation and have exactly the same solution S(x).
Thus one can interpret i in two ways. One can consider i to be a TAL volume scaling (relative to the base case of i = 1) that represents a uniform inflation (i > 1) or deflation (i< 1) of the TAL lumen all along the ith TAL, relative to a reference cross-sectional area Ao; if i = 1, the cross-sectional area is considered to be Ao. The intratubular perimeter length corresponding to the reference area Ao is denoted c; that perimeter length is considered to be fixed, so that transport capacity per unit TAL length remains constant when i is varied. Thus although the cross section may vary from simulation to simulation (or vary as a function of time), total TAL transport capacity is unchanged. This interpretation of i would be appropriate to describe the effects of a short-term change in intrarenal pressure that changes tubular volume without changing TAL transport capacity.
Alternatively, one can consider i to be a scaling of the TAL flow rate in which the transport characteristics of the TAL have been also scaled in such a way that the same steady-state concentration profile S(x) is preserved. This situation might arise from a long-term adaptation of the TAL to a change in flow rate, as, for example, from a long-term change in either SNGFR or fractional absorption from the proximal tubule.
Both interpretations of i can be unified under the concept of the steady-state transit time interval for the TAL; this interval is the time required for a water molecule to travel up the TAL, from the TAL entrance to the MD, at a constant speed, assuming plug flow. For our steady-state base-case parameters, that transit time, denoted to, equals the TAL volume, AoL, divided by the steady-state flow Fop. However, for a change in i, under either of the interpretations advanced above, the steady-state transit time becomes ito and thus i may be considered to be a dimensionless scaling of the base-case steady-state TAL transit time to.
Parameters
Base-case values for model parameters are given in Table 1. Extratubular concentration is specified by Ce(x) = Co [A1 exp(–A3 x) + A2], where A1 = (1 – Ce(L)/Co)/(1 – exp(–A3)), A2 = 1 – A1, and A3 = 2, and where Ce(L) corresponds to a cortical interstitial concentration of 150 mM. Graphs for Ce(x) and the steady-state profile S(x) were given in Fig. 1 of Ref. 26. The steady-state concentration Cop corresponds to S(L).
|
of 2 s to be followed by a transition interval of 3 s, providing an overall delay i = + /2 of 3.5 s, in approximate agreement with experiments (5). The overall delay i was varied by varying , whereas was considered to be a fixed parameter in all cases (see supplementary material, item 1; all supplementary material in this paper can be found at http://ajprenal.physiology.org/cgi/content/full/00048.2005.DC1). Analysis of Model Behavior
We used two methods to investigate model behavior: 1) analysis of the model's characteristic equation and 2) direct computation of numerical solutions of the model equations. The characteristic equation contains information that allows systematic predictions of generic model behavior as a function of model parameters. By "generic behavior," we mean, e.g., whether the model's long-time behavior will approximate a time-independent, stable steady state (i.e., a SS) or a sustained oscillation. A numerical solution of the model equations allows one to determine model behavior for a specific set of parameters. These two methods are complementary, in that the first provides helpful, but not always certain, predictions for many parameter combinations, whereas the second is a more nearly certain predictor for behavior, but only for a particular set of parameters. Moreover, some features of the predictions by both methods should be entirely consistent, so each method provides a validation of the other. The methods that were used to compute numerical solutions from the model equations are summarized in the APPENDIX; the APPENDIX also contains a synopsis of the procedures that were used to confirm or identify bifurcation boundaries.
For this study, we used only the characteristic equation for a single, uncoupled nephron (i,j = 0 for all j in Eq. 2). We have previously described how such a characteristic equation is obtained (26). Briefly, the nonlinear model equations (Eqs. 1–5) are combined into a single nonlinear equation and linearized to obtain a linear partial differential equation (PDE). An expression representing potential steady-state and oscillatory solutions is substituted into that linear PDE, and from the resulting equation, the conditions for the potential solutions that satisfy the linearized equation may be derived. Those conditions are embodied in the characteristic equation, which is given by
 | (9) |
p and have been rendered dimensionless through division by the base-case steady-state TAL transit time to. However, , already dimensionless, is unchanged. The additional parameter appearing in Eq. 9 is the feedback-loop gain magnitude, a dimensionless quantity that is the product of two factors (26, 27): 1) the dimensionless slope of the TGF response curve at the steady-state operating point (Fop, Cop), which is given in nondimensional form by K1K2 (see APPENDIX); and 2) the magnitude of the nondimensional slope of the TAL tubular fluid Cl– concentration profile alongside the MD [i.e., –S'(L), but expressed in nondimensional form]. The two factors that make up depend on the model parameters. In this study, K1 and –S'(L) are assumed to have fixed values; to obtain particular values of , K2 is varied.
A solution of the characteristic equation is a complex number 
that depends on the values assigned to the model parameters. Because is a complex number, it can be written in the form
 |
is the real part of and Im is the imaginary part of . The real and imaginary parts of are indicators of the strength and angular frequency, respectively, of a solution to the linearized model equations. A solution to the characteristic equation is not unique; rather it is a number appearing in an infinite series of solutions 1, 2, 3, .... However, as shown in RESULTS, only 1 and 2 appear likely to have a tangible impact on TGF system behavior. Because a solution of the characteristic equation may reasonably be considered to be an eigenvalue of a particular differential operator (34), we will refer to such solutions as eigenvalues (see supplementary material, item 2).
When backleak permeability P is set to zero, the characteristic equation given in Eq. 9 has a simpler form, given by
 | (10) |
The characteristic equation makes predictions about the stable solutions of the linearized model equations, as a function of parameter values (26). If, for a particular set of parameters, Re n < 0 for all indices n, then the only stable solution is an SS. If, on the other hand, Re n > 0 for one, and only one index n, the only stable solution is a regular oscillation having frequency, for that index n, of (Im n)/(2to). For cases where Re n > 0 for multiple values of the index n, there may be multiple stable oscillatory solutions. In such case, the nature of the solutions must be investigated directly by solving the model equations to determine potential solution behaviors, behaviors that may vary for differing initial conditions. The value of the gain that corresponds to the transition from a case where Re n < 0, for all indices n, to a case where, for at least one index n, Re n > 0, is a called a critical gain, denoted c (thus for that n, Re n = 0 at the point of transition). For the base-case parameters, c = 3.24 where Re 1 = 0. (For additional information regarding the applicability and limitations of the characteristic equation, see supplementary material, item 3).
|
RESULTS |
|---|
|
TOPABSTRACTMATHEMATICAL MODELRESULTSDISCUSSIONAPPENDIXGRANTSREFERENCES |
|---|
As we now demonstrate, for the case of nonzero TAL backleak permeability, the qualitative behavior of our TGF system model is very different from that previously inferred by us from the case of no backleak permeability (26).
The diagrams in Fig. 3, A1 and B1, for a parameter region in the the - plane, are based entirely on our analysis of the characteristic equations Eqs. 10 and 9, for no backleak (P = 0) and for nonzero backleak (P = 1.5 x 10–5 cm/s), respectively. In both cases, steady-state TAL transit time was fixed at its base-case value, thus = 1. Throughout RESULTS, (or i), which characterizes the TGF signal delay at the MD, designates the sum of the nondimensional values of p (or ) and (which correspond to the durations of pure delay p and a distributed delay , as reported in Ref. 5). The gain magnitude (or i) characterizes the strength of the TGF response (see MATHEMATICAL MODEL), as in experimental studies (15, 47). Solid curves in these diagrams mark loci where the real part of an eigenvalue n is zero and the sign of the real part changes from negative to positive, as indicated by the inequalities in the figures. Dashed curves mark loci where the real part of an eigenvalue overtakes the real part of another eigenvalue. Thus, for example, for values of and that fall above the curves marked "Re 2 = Re 1," Re 2 exceeds Re 1.
Although, as will be apparent below, the diagrams in Fig. 3, A1 and B1, provide a good guide to the generic solutions of our model equations, they are nonetheless based on the linearized forms of those equations, and detailed model behavior must be investigated and confirmed by means of simulations in which numerical solutions to the unaltered model equations are computed. Figure 3, A2 and B2, which portray the same region of the - plane as do Fig. 3, A1 and B1, indicate the behaviors of model solutions based on numerical solutions, which were computed as described in the APPENDIX. Figure 3, A2 and B2, which are bifurcation diagrams, show the true model behaviors, which are substantially consistent with the indications of Fig. 3, A1 and B1.
Figure 3A2, for the no-backleak case, corresponds to Fig. 3A1. For this case, only two generic stable behaviors are observed: if the values of and fall below the curve labeled "Re 1 = 0," then, for any initial conditions, or for any transient perturbation of a steady-state solution, the model solutions always more and more closely approximate the SS solution. However, for values of and above the curve Re 1 = 0, the only stable solution that was observed was a regular oscillation that converged to a fixed, periodic trajectory; the period and amplitude of that oscillation vary as a function of and . This oscillation, which is an LCO, has a frequency f1 that approximates the frequencies [viz., (Im 1)/(2to)] that can be found by solving the characteristic equation along the locus where Re 1 = 0. Because the LCO can be characterized by the frequency f1, we call it an "f1-LCO"; this LCO may be considered to be the eigenmode associated with the eigenvalue 1. The f1-LCO solutions for SNGFR and for MD Cl– concentration, at the point labeled "A3" in Fig. 3A2 ( = 0.135, = 8.00), are shown in Fig. 3A3. A large number of numerical simulations, using sundry perturbations (which included the transient addition of oscillations of various frequencies to TAL flow), revealed no other stable solution for the region above the loci of Re 1 = 0.
However, when TAL backleak permeability P was increased from 0 to P = 1.5 x 10–5 cm/s, the diagram obtained by solving the characteristic equation (Eq. 9), and the observed model behavior, both changed dramatically, as shown in Fig. 3, B1 and B2, which are analogous to Fig. 3, A1 and A2. In Fig. 3B1 the curves for the conditions Re 2 = 0 and Re 3 = 0 cross the curve for Re 1 = 0. Our numerical simulations revealed that such crossings can have important consequences for the nature of the stable solutions that can occur when the real part of at least one eigenvalue is positive. As in Fig. 3A1, a bifurcation locus Re 1 = 0 marks the transition from an SS to an f1-LCO. However, this bifurcation locus is now restricted to the portion of the diagram that lies to the right of the crossing of Re 1 = 0 and Re 2 = 0. Moreover, numerical solutions revealed a curve, labeled "Empirical Boundary," that lies above the curve for Re 2 = Re 1 and that marks the boundary between the f1-LCO and a different type of dynamic behavior: in the parameter region between the Empirical Boundary and the curve for Re 2 = 0, there is a potential stable LCO solution that corresponds to 2; we call that solution an "f2-LCO." The numerical solutions showed that an f2-LCO has a frequency approximately twice that of an f1-LCO. In the top left where only Re 2 > 0, which is labeled by "f2," the only stable solution is an f2-LCO. Moreover, above the curves for the Empirical Boundary and for Re 1 = 0, two stable LCO, of differing frequencies, can be elicited, depending on initial conditions. One of these solutions is an f1-LCO, and the other is an f2-LCO; we call this behavior an "f1, f2-LCO." Thus this parameter region, where the real parts of both 1 and 2 are both positive, and where Re 2 exceeds Re 1, exhibits bistability with respect to the f1 and f2 solutions. We did not observe any f3-LCO. The region in which Re 3 > Re 2 is completely contained within the region where Re 2 > 0, and thus the relationship of the real parts of these eigenvalues is analogous to that for the real parts of 2 and 1 in the case of P = 0 shown in Fig. 3A1. The corresponding solution behaviors for these analogous cases are consistent in that f2-LCO were not observed in the parameter region of Fig. 3A2, and f3-LCO were not observed in Fig. 3B2.
The gray disk, which has its center at values of 
0.2228 and 3.24 (and that thus lies on the curve for Re 1 = 0), represents an approximate region for typical parameters in normotensive rats, which, based on our previous studies, appears to be near the bifurcation boundary from the SS to f1-LCO (26, 30).
LCO in SNGFR and tubular fluid Cl– concentration at the MD that were computed for points A3 and B3 are shown in Fig. 3, A2 and B2, respectively. The delay and gain for both points were set to (, ) = (0.135, 8.00). In the zero backleak case, SNGFR and Cl– concentration, shown in Fig. 3A3, both oscillate with an f1-LCO frequency of 50.1 mHz. However, for the nonzero backleak case, the analogous point B3 lies in the region for f1, f2-LCO; thus two different stable solutions are possible. In Fig. 3B3, we exhibit the f2-LCO, which oscillates with a frequency of 101 mHz.
Three considerations may clarify these findings. First, in the regions where LCO are the stable solutions, the frequencies of those solutions vary as a function of , , and . Moreover, within a region where a particular type of LCO is observed, its frequency is much more sensitive to variations in time delays, and thus to changes in and 1/, than to changes in the gain . Thus, for example, an f1-LCO has a frequency along a vertical line (i.e., a fixed value of ) in Fig. 3, A2 or B2, that closely approximates the frequency predicted by the characteristic equation at the point where that line crosses the bifurcation curve for Re 1 = 0. An analogous result was found for f2-LCO.
Second, although an f2-LCO has a frequency about twice that of an f1-LCO (provided that the values of and in the f2 case do not differ or differ little from the corresponding values in the f1 case), these frequencies are not in precise harmonic relationship. Thus even though the frequencies f1 and f2 correspond to the first two eigenmodes for resonate standing waves within the TAL, they are differentially influenced by their respective bifurcations and by nonlinearities inherent in the full model equations.
Finally, for ease of comparison among our TGF investigations, we have used a uniform set of base-case model parameters (26–30, 35). These parameters, which were based on experimental evidence (26), yield an f1-LCO at 45 mHz (29), which is in the upper frequency range of experimentally measured LCO oscillations in rats (19). However, our model framework is generic, in the sense that the results depend on dimensionless lumped parameters (e.g., , , and ), and differing combinations of those parameters predict oscillations of varying frequencies. Indeed, as shown in Fig. 5 (gray dashed lines) and exemplified in numerical results of Figs. 6 and 7, f1-LCO frequencies can change substantially with changes in TAL transit time, which is characterized by , and these frequencies encompass the physiological range.
TAL Transit Time May Be an Important Bifurcation Parameter
We varied TAL transit time by varying , which may be considered to scale either the TAL cross-sectional area (and hence TAL volume) or the base-case steady-state TAL flow Fop (see MATHEMATICAL MODEL). To investigate the impact of in the context of the previously identified bifurcation parameters and , we extended the planar parameter region in Fig. 3 to three dimensions, with axes oriented as shown in the inset in Fig. 4A1. Figure 4, A1 and A2, shows diagrams corresponding to a portion of the (1/)- parameter plane for fixed feedback delay = 0.2228 (dimensional value, 3.5 s), and for P = 1.5 x 10–5 cm/s. These diagrams are analogous to, and share a qualitative similarity with, the - diagrams given by Fig. 3, B1 and B2; the interpretation of Fig. 4, A1 and A2, is therefore similar to the interpretation of Fig. 3, B1 and B2.
Below the curves in Fig. 4, A1 and A2, the only stable solution is a time-independent steady state (labeled "Re n< 0" and "SS" in Fig. 4, A1 and A2, respectively). In Fig. 4A2, the parameter regions labeled "f1", "f2," and "f1 or f2" have the same model solution properties as described for the analogously labeled regions in Fig. 3B2. In addition, however, Fig. 4A2 contains a distinct region, labeled by "f3," where Re 3 is sufficiently large that an oscillation having a frequency of 3f1 is stable (note that in Fig. 4A1, the loci for Re 3 and Re 2 cross, so there is a region in which only f3-LCO are stable, analogous to the region in Fig. 3B2 where only f2-LCO are stable).
These results indicate that is also a bifurcation parameter, and that 1/ has effects similar to those of (compare Figs. 3B2 and 4A2). Both and can be considered delay parameters, because characterizes the delay in TGF signal transmission at the JGA whereas is proportional to TAL transit time. However, we have previously shown, in the simpler model context of zero TAL backleak, that the bifurcation behavior depends on the ratio of the delay at the JGA to the TAL transit time (26); in the present context, that ratio corresponds to the ratio /. This explains the strong similarities between the effects of changes in and of changes in 1/.
Multiple Stable Solutions and Multistability May Introduce Spectral Complexity
Figure 4, B1 and B2, shows simulated oscillations of SNGFR and the associated power spectrum corresponding to point B in Fig. 4A2, where (1/, ) = (2/3, 9/2). (The power spectra reported in this study were produced by the same method as described in Ref. 29; however, no broadband perturbations were used in the present study.) Because point B lies below the curves, where the real parts of all eigenvalues are negative, the stable solution is a SS: a solution that is transiently perturbed will always return more and more closely to that steady state, provided that no further perturbations are applied. Such a return is an evolving process of successive approximation, like that exemplified in Fig. 4B1. The model's steady-state SNGFR is 30 nl/min, and the power spectrum for that steady state shows that all frequency components are zero (all power spectra presented herein are based on the deviation from the steady-state SNGFR). Point C, where (1/, ) = (2/3, 11/2), corresponds to a single stable solution, the f1-LCO. Simulated SNGFR oscillations and the associated power spectrum for point C are shown in Fig. 4, C1 and C2. The oscillation frequency is 35 mHz, and because the oscillation is nonsinusoidal it has harmonics (29), the first of which is at 70 mHz. Point D, at (1/, ) = (2/3, 13/2), lies in the bistable region of f1, f2-LCO. Simulated SNGFR oscillations (for the f2-LCO) and the associated power spectrum are shown in Fig. 4, D1 and D2. The frequency of this oscillation is 66 mHz, which is nearly 2 x 35 mHz. It is noteworthy that three very different stable model solutions exist for values of gain that span about 2/9 of the reported experimental range of 1.5 to 9.9 (15). Indeed, examination of Fig. 4A2 (near points B, C, and D) shows that the three differing solution types (viz., SS, f1-LCO, and f2-LCO) can be obtained for an even smaller range of than was used in the example presented here. Hence, in this parameter region, even small stochastic variations in model parameters can result in switching between oscillatory modes.
It is noteworthy that point D in Fig. 4A2 is in a bistable "f1 or f2" region having two stable solutions, viz., an f1-LCO and an f2-LCO (an analogous region was shown in Fig. 3B2). The f1-LCO at point D is very similar to the f1-LCO found at point C (see Fig. 4C1); the f2-LCO at point D was shown in Fig. 4D1. The two stable solutions at point D are an instance of multistability. Whether the f1-LCO or the f2-LCO emerges as the stable mode depends on initial conditions (see text describing Fig. 7 below); moreover, perturbations can cause switching from one mode to another. Because perturbations and parameter variations can occur within short time intervals both multistability and multiple stable solutions occurring within a limited parameter regime may result in multiple strong peaks in power spectra.
Parameter Nonstationarity Can Introduce Spectral Complexity
Figure 5, A and B, shows -(1/) bifurcation diagrams for = 5 and = 7, respectively; in both cases, P = 1.5 x 10–5 cm/s. These diagrams correspond to a portion of the -(1/) plane, which is perpendicular to the -axis (see inset in Fig. 4A1).
In Fig. 5A, the region above the curve labeled by "Re 1 = 0" corresponds only to stable solutions exhibiting an f1-LCO; the region below the curve corresponds only to the time-independent steady-state solution, SS. The solid dot, at (o, 1/o) = (0.2228, 0.7282), lies just within the region having an f1-LCO solution. The three gray dashed curves mark the -(1/) loci for which the frequencies of f1-LCO are predicted by the characteristic equation to be 20, 35, and 50 mHz.
For fixed gain ( = 5), time-dependent variations were introduced in the values of 1/ and . In the " case," was fixed at o, but was varied linearly, over an interval of 4,096 s (68 min), at a uniform rate of change, from 80 to 120% of o. Subsequently, in the " case," was fixed at o, but was varied linearly, over the same time interval, at a uniform rate of change, from 120 to 80% of . In Fig. 6, A1 and A2, for linear and logarithmic ordinate scales, respectively, power spectra for these simulations (solid curves for and dashed curves for ) and for the stationary parameter case (gray curves) of (, ) = (o, o), are shown. In the simulations with variation in or , peaks in the power spectra are broadened; i.e., SNGFR oscillations contain a wider range of frequencies, relative to the stationary case. Peaks in the power spectra of the case are slightly to the right of stationary case, because the higher frequencies correspond to larger 1/'s, which in turn correspond to parameter values farther away from the SS region; thus larger 1/ corresponds to oscillations of larger amplitude. In contrast, the peaks in the power spectra of the case are slightly to the left of the base case. In this case, the lower frequencies correspond to larger 's and thus to oscillations of larger amplitude. (Note the good agreement between the frequencies in the power spectra, which were computed from the full model, and the frequencies predicted by the characteristic equation, as shown by the gray dashed curves in Fig. 5A.)
In another simulation, also for = 5, was fixed, but was set to be 80% of o for approximately half of the time interval; then, in two nondimensional time units (2to = 31.42 s), was increased linearly to 120% of o; and for the remainder of the time interval, = 120% of o. The analogous simulation was also conducted for . Power spectra obtained for the case and for the stationary case (, ) = (o, o) are shown in Fig. 6, A3 and A4, in linear and logarithmic ordinate scales, respectively. Two large peaks, rather than one broadened peak, were obtained for the case; these peaks, at 43 and 37 mHz, correspond to = 80% of 0 and 120% of o, respectively. The 43-mHz peak is larger because oscillations have larger amplitude for smaller 's. In contrast to the arrow, which lies approximately perpendicular to the "35 mHz" curve, the arrow lies approximately parallel to the 35 mHz curve. Thus the two peaks obtained by varying (not shown) lie much closer to each other, compared with the case.
The experiments described above were repeated for = 7; the bifurcation diagram for this case is shown in Fig. 5B, and corresponding results are in Fig. 7. Because the base parameters (dot, Fig. 5B) lie within the region of f1-LCO or f2-LCO, one of these two stable solutions had to be chosen; we chose the f2-LCO of 70 mHz (gray curves, Fig. 7). However, because the f1, f2-LCO region is bistable, the mode corresponding to an f1-LCO was also excited, as can be seen in the logarithmic power spectrum (gray curve, Fig. 7B2). As in the previous experiments, and were varied linearly, at a constant rate of change, corresponding to paths along the arrows indicated in Fig. 5B, over a interval of 68 min. The tails of both arrows lie within the f1 region, whereas the heads lie within the bistable region. As shown in Fig. 7, B1 and B2, f1-LCO were obtained for both cases. As in the = 5 case, the peaks (37 mHz) were broadened, because the oscillations involve a range of frequencies, with the peak wider than the peak. It is noteworthy that if and were varied along different arrows such that the parameters began in the f2 region and moved into the bistable region, or if an f1-LCO in the bistable region was perturbed using a sinusoidal wave having a frequency approximating the frequency of the f2-LCO, the resulting solutions were f2-LCO (results not shown).
In the next set of simulations, was set to be 120% of o for approximately half of the time interval; then, in two nondimensional time units, was decreased at a constant rate to 80% of o; and for the remainder of the interval, = 80% of o. The power spectrum associated with the resulting oscillations has three prominent peaks (Fig. 7, B3 and B4). The 35-mHz peak corresponds to the f1-LCO obtained for = 120% of o, which lies within the f1 region; the 70-mHz peak is its first harmonic. The 67-mHz peak, which is also excited, corresponds to the f2-LCO obtained for = 80% of o, which lies within the bistable region. A goal of this numerical experiment was to exhibit distinct peaks at approximately the fundamental frequency and the first harmonic, similar to the relationships seen in some power spectra from physiological experiments (see DISCUSSION). Thus it is desirable in the simulation for the fundamental frequency to remain relatively constant as the parameter is varied. Because, as we have seen, the same variations in (as a percentage of the base case) introduce smaller changes in the fundamental frequency than do variations in , results are shown in Fig. 7, B3 and B4, for the case rather than for the case.
Nephron-to-Nephron Coupling May Introduce Spectral Complexity
We investigated the effects of coupling a model nephron having parameters in the bistable ("f1 or f2") region with nephrons in the f1 region. The bifurcation conditions for a system of coupled nephrons differ somewhat from those for a single nephron (38). However, for simplicity, the bifurcation conditions for a single nephron were used in this study as a guide for predicting coupled behavior. First, two nephrons, designated 1 and 2, were coupled. The parameters for the nephrons are given at the top of Table 2; the coupling parameter ij was set to 0.2 for both nephrons. Figure 8, A1–A3, show oscillations in SNGFR for nephron 1 and corresponding power spectra, in linear and logarithmic scales. In this nephron, the frequencies of 27 and 60 mHz were most excited. These frequencies correspond to the fundamental frequency of nephron 2, which is in the f1 region, and the f2-LCO of nephron 1, which is in the bistable region. The fundamental frequency of nephron 1 (33 mHz) and the first harmonic of nephron 2 (54 mHz) were also excited.
|
ij was set to 0.2 in all couplings among the nephrons. With these parameters, nephron 1 is in the bistable region, whereas nephrons 2 and 3 are in the f1-LCO region. Figures 8B1-8B3 show oscillations in SNGFR for nephron 1 and corresponding power spectra. In this case, the most excited frequencies are 27 mHz (f1-LCO of nephron 3), 42 mHz (f1-LCO of nephron 2), and 62 mHz (f2-LCO of nephron 1). The f1 frequency of nephron 1 (34 mHz) and the first harmonics of nephrons 2 and 3 (84 and 55 mHz) were also excited. The f1-LCO's in nephrons 2 and 3 do not synchronize because their delays are sufficiently different (23, 38). These results predict that irregular oscillations in SNGFR and complexities in corresponding power spectra can be obtained via internephron coupling. Combined Effects of Bifurcations, Coupling, and Parameter Lability
The resul
