Abstract
The functional units of the kidney, called nephrons, utilize mechanisms that allow the individual nephron to regulate the incoming blood flow in response to fluctuations in the arterial pressure. This regulation tends to be unstable and to generate selfsustained oscillations, perioddoubling bifurcations, modelocking and other nonlinear dynamic phenomena in the tubular pressures and flows. Using a simplified nephron model, the paper examines how the regulatory mechanisms react to an external periodic variation in arterial pressure near a region of resonance with one of the internally generated modelocked cycles. We show how the stable and unstable resonance cycles generated in this response undergo interconnected cascades of perioddoubling bifurcations and how each period doubling leads to the formation of a new pair of saddlenode bifurcation curves along the edges of the resonance zone. We also show how period doubling of the resonance cycles is accompanied by a torusdoubling process in the quasiperiodic regime that exists outside of the resonance zone.
1. Introduction
Living organisms operate under farfromequilibrium conditions [1]. This implies that many of the feedback regulations that control the biological processes at different time and space scales are unstable and produce selfsustained oscillations and other complex nonlinear dynamic phenomena.
The regulation of the incoming blood flow to the individual functional unit of the kidney, for instance, involves at least two different mechanisms that tend to produce oscillatory dynamics [2,3]. The tubuloglomerular feedback (TGF) mechanism produces selfsustained oscillations in the 40 mHz regime because of a delay associated with the flow of fluid through the loop of Henle, and the myogenic mechanism produces vasomotoric oscillations in the 200 mHz regime in connection with the synchronized activation of the smooth muscle cells in the arteriolar wall in response to increasing transmural pressures. Detailed statistical analyses [4,5] have clearly shown how the oscillatory mechanisms interact to produce different forms of synchronization between the two modes. Similar analyses have revealed episodes of perioddoubling dynamics in nearly 50 per cent of the available time series for normotensive rats. In rats with 2kidney, 1clip Goldblatt hypertension and in spontaneously hypertensive rats, one observes chaotic dynamics caused, presumably, by the combination of an increased sensitivity in the myogenic response and a stronger feedback gain for the TGF mechanism.
We have recently examined the bifurcation structure associated with the Ctype perioddoubling transition to chaos in a periodically forced Rössler oscillator [6,7]. This transition is characteristic of multidimensional systems where the perioddoubling bifurcations take place in the presence of transitions between states of synchronized and nonsynchronized dynamics. We have demonstrated how the stable (node) and unstable (saddle) cycles in the forced Rössler system undergo interconnected cascades of perioddoubling bifurcations. Near the edge of the synchronization zone, node and saddle cycles with the same periodicity approach one another and at the very edge, the perioddoubling bifurcations occur simultaneously. Away from the edge of synchronization, the node solution in the forced Rössler system tends to bifurcate first, i.e. for lower values of the parameter that controls the nonlinearity of the system. More importantly, however, these period doublings take place in a direction transverse to the unstable manifold of the saddle cycle that connects the resonance modes, and perioddoubling transitions of the periodic modes are accompanied by a torusdoubling process in the quasiperiodic region outside the resonance zone.
Each pair of perioddoubling bifurcations gives rise to a new set of saddlenode bifurcation curves along the sides of the resonance zone. As the perioddoubling cascade unfolds, the edges of the resonance zone thus accumulate more and more of almost parallel saddlenode bifurcation curves, each defining the boundaries of the resonance zone for a particular level of the interconnected perioddoubling cascades. The new saddlenode bifurcation curves are supported by the corresponding perioddoubling curves, but are formed a little away from the existing zone edge. This produces a gap between the saddlenode bifurcation curves, and additional local and global bifurcations are required to complete the zone boundary [7]. We have also demonstrated how the saddlenode bifurcations along the edge of the resonance zone accumulate and, in an alternating manner, approach a welldefined threshold curve for the transition between phasesynchronized chaos and nonsynchronous chaos.
The purpose of the present paper is to demonstrate how a similar structure arises in a model of the physiological mechanisms by which the individual functional unit (nephron) of the kidney regulates the incoming blood flow to compensate for variations in arterial pressure arising, for instance, from changing levels of physical activity. As mentioned above, the individual nephron generates two different and clearly recognizable internal modes of oscillation. Interaction between these modes leads to a variety of additional nonlinear dynamic phenomena, including modelocking, multistability, period doubling and transitions to deterministic chaos. Over the years, these phenomena have been extensively studied both in experiments on anaesthetized rats [2,8] and through model simulations [9,10]. In the frequency regime of interest for the present discussion, the nephron autoregulation functions as a mechanical highpass filter: rapid fluctuations in arterial pressure are allowed to propagate into the tubular system, but the delicate regulatory processes that take place in the distal tubule are protected from more lasting variations.
After a short introduction to the considered nephron model [10,11] and a discussion of a few of the underlying experiments, the paper first presents a twodimensional bifurcation diagram in order to clarify the organization of the internal modes in the physiologically relevant regime. We then apply an external periodic forcing and determine the resulting bifurcation structure as a function of the afferent arteriolar resistance and the frequency of the applied forcing. Although the details of the organization are somewhat different, we rediscover many of the characteristic features associated with the Ctype criticality in the forced Rössler system: the interconnected cascades of perioddoubling bifurcations for the resonant node and saddle cycles, the accumulation of saddlenode bifurcation curves along the edges of the resonance tongue, the gaps that open up between the saddlenode bifurcation curves and the torusdoubling process that takes place outside of the resonance zone.
2. Nephron pressure and flow regulation
It is well known that the kidneys function as a filter to remove metabolic end products and regulate the excretion of water and salts so as to maintain a proper volume and osmolality of the blood and appropriate conditions for the cells in general. The kidneys also play a part in the regulation of the blood pressure, both directly by controlling the extracellular fluid volume and through the production of hormones that, together with hormones from other organs, regulate the peripheral flow resistance of the vascular system. At the same time, to protect their own function, the individual nephrons in the kidney utilize two different mechanisms that can compensate for variations in the arterial blood pressure and maintain a relatively constant filtration rate.
The myogenic mechanism depends on the inherent propensity of the smooth muscle cells in the afferent arteriolar wall to contract in response to an increasing pressure difference across the arteriolar wall [3,12]. This contraction causes the flow resistance to increase and, hence, leads to lower glomerular pressure and reduced filtration. The TGF mechanism, on the other hand, depends on a signal from specialized cells (the macula densa cells) near the end of the loop of Henle. These cells respond to changes in the salt concentration in the fluid leaving the loop of Henle. A high rate of glomerular filtration will lead to a faster flow through the loop of Henle, to incomplete reabsorption of salt from the tubular fluid, and to increasing concentrations of salt at the macula densa. This produces a signal to the smooth muscle cells in the arteriolar wall to contract, thus causing the filtration to go down.
As mentioned above, experiments on rats [2,8] have demonstrated that the TGF mechanism is unstable and produces clearly detectable selfsustained oscillations in the tubular pressures and flows with periods in the 30–40 s range. This relatively long periodicity is directly related to the time of approximately 15 s that it takes for the tubular fluid to flow through the loop of Henle. The myogenic (or vasomotoric) mechanism is also capable of producing selfsustained oscillations through synchronization of the cytoplasmic Ca^{2+} dynamics among the smooth muscle cells in the arteriolar wall as the muscular activation increases [3,13]. As the entire signalling pathway for the myogenic mechanism is present with the individual smooth muscle cell, its response time is faster than that of the TGF mechanism. As a consequence, the corresponding myogenic oscillations have a period of only 5–10 s [14,15]. The two oscillatory modes both work through activation of the smooth muscle cells in the arteriolar wall. This allows the oscillations to interact and to produce modelocking and other nonlinear dynamic phenomena [4,5]. The frequency ratios most commonly observed between the two modes are 1 :4, 1 : 5 and 1 : 6 [16].
The autoregulatory system involves a number of nonlinear relations. A main nonlinearity acting to constrain the amplitudes of the oscillatory modes is associated with the limited dynamic range for the arteriolar contraction. This range, which differs between normotensive and spontaneously hypertensive rats [14,15], can be determined through openloop experiments where the rate of glomerular filtration is measured as a function of the rate of infusion of artificial tubular fluid into the loop of Henle, while preventing flow through the proximal tubule by means of a wax seal [14,17]. The slope of this relation determines the feedback gain factor α for the TGF mode. This gain factor is typically 30–50% larger for spontaneously hypertensive rats than for rats with normal blood pressure [14,15]. Interaction between the two regulatory mechanisms is also found to be significantly stronger in hypertensive than in normotensive rats [5,18].
Figure 1 shows a typical example of the pressure oscillations observed in the proximal tubule of a normotensive rat. The mean tubular pressure in this experiment is about 7 mmHg. However, this pressure is modulated by relatively slow TGFmediated oscillations with amplitudes of the order of 1–2 mmHg. The faster myogenic oscillations reveal themselves as a ripple on top of the TGF oscillations. Wavelet analysis clearly confirms their existence and can also detect how their amplitudes and frequencies are modified by the amplitude of the TGF oscillations [5]. Closer inspection of the time trace in figure 1a reveals an episode of period2 dynamics, extending approximately from time 200 s to 500 s, where the tubular pressure alternates between shallow and deep minima. This part of the time trace is magnified in figure 1b. The spectral distribution shown in figure 1c also provides clear evidence for the presence of subharmonic components corresponding to half the frequency (or twice the period) of the TGF signal. Approximately 50 per cent of our time series for normotensive rats show evidence of period doubling.
The tubular pressure variations in hypertensive rats are much more irregular, even though the amplitudes are often smaller. This is illustrated for a spontaneously hypertensive rat in figure 2. The differences between the two types of dynamics can probably not be related to a single factor. We have already indicated that both the feedback gain and the interaction between the two regulatory modes are higher for hypertensive than for normotensive rats. It is also likely that the higher arterial pressure has shifted the working point on the TGF feedback curve and enhanced the sensitivity of the myogenic response. The main conclusions that can be drawn from the time series for hypertensive rats are that the myogenic component plays a stronger role, that synchronization between the two modes is less common, and that the experimentally observed pressure variations for most practical purposes can be considered as chaotic [18,19]. The spectral distribution in figure 2c demonstrates the presence of three to four subharmonic components to the TGFmediated oscillation.
Over the years, a range of different experiments have been performed to determine the precise mechanisms underlying the two modes of oscillation [8,15,17] and to explicitly measure the different parameters and nonlinear relations. The proximal and distal tubule of superficial nephrons can be detected on the surface of the kidney and the tubular pressures can be measured by means of a small glass pipette [2]. It is also possible to examine the effect of a perturbation to the system by infusing artificial tubular fluid into the loop of Henle at rates of 5–10 nl min^{−1}. Besides openloop experiments to measure the feedback relation, experiments have been performed to determine the tubular flow resistance, the damping of the pressure oscillations along the loop of Henle and the delay in the fluid flow [20]. Experiments have also been made to determine the amplitude and phase of the NaCl concentration oscillations near the macula densa and to determine the degree of crosstalk between neighbouring nephrons [21].
The ability of the renal autoregulation to handle external pressure fluctuations can be tested by applying a forcing signal to the arterial pressure while simultaneously recording the spectral distributions of the variations observed in this pressure and in the renal blood flow. In practice, the experiment can be performed [15] by connecting a computeroperated pump that generates broadband fluctuations at the distal end of the abdominal aorta through a bloodfilled cannula.
Figure 3 shows an example of the frequency characteristics observed in such an experiment and covering the range from 1 mHz to 1 Hz. This characteristic clearly shows the damping of the pressure oscillations in the nephron blood flow at frequencies lower than 30 mHz. Also revealed are the resonances characteristic of the autoregulatory system: the relatively slow TGF mode at about 40 mHz and the faster myogenic mode around 200 mHz. We conclude that the nephron autoregulation functions as a highpass filter that protects the nephron against more longterm variations in the arterial pressure. Fluctuations above the range of the myogenic oscillations are likely to be damped out by the dissipative processes associated with the fluid flow through the loop of Henle. The TGF mechanism is unique to the nephrons. By virtue of the enormous blood flow that the kidneys have to handle, the TGF mechanism is required as a reinforcement of the myogenic mechanism. The frequency response of the TGF mechanism is restricted, though, by the delay in the Henle flow and, while reduced in amplitude by a factor of the order of two, the TGF oscillations are still present in the distal tubular pressure and salt concentration. By allowing the feedback to be oscillatory, the system presumably achieves the fastest possible reaction with the given delay.
3. Simplified single nephron model
Together with the experimental work we have constructed a number of different models of the autoregulatory system [10,14,19,22]. These models each emphasize particular aspects of the problem such as, for instance, the absorption processes in the loop of Henle [2] or the mechanisms involved in the feedback from the macula densa cells to the smooth muscle cells in the arteriolar wall [22]. In the present paper, we shall use the simplified model described by Mosekilde et al. [11]. This model is particularly useful for more detailed bifurcation studies. The model integrates the most essential aspects of nephron autoregulation into a consistent picture and, over the years, it has been able to predict several phenomena that have subsequently been detected in the experimental data. The model has also been used to examine different types of nephron–nephron interaction [23–25].
The first component in the model is a simple conservation equation that relates changes in the proximal tubular pressure to the rate of glomerular filtration, the absorption in the proximal tubule and the flow into the loop of Henle. This is combined with a number of algebraic equations that determine the rate of glomerular filtration in terms of the arterial pressure and the afferent arteriolar resistance and determine the rate of flow into the loop of Henle in terms of the proximal tubular pressure, the distal pressure and the tubular flow resistance. Calculation of the rate of glomerular filtration also involves an equation that relates the protein osmotic pressure to the protein concentration in the arteriolar blood. The TGFmediated variation in the afferent arteriolar resistance is determined from the measured openloop feedback relation, but also includes a smooth delay to represent the time it takes for the fluid to pass the loop of Henle. Finally, we have applied a couple of firstorder differential equations to simulate the myogenic mechanism that generates the fast oscillatory component. In the form we use the model here, these equations describe a damped oscillator, but the myogenic oscillations are continuously excited through the variation in the arteriolar resistance caused by the TGF mechanism.
Figure 4 shows a twodimensional bifurcation diagram for the unforced nephron model in a parameter plane spanned by the delay time T in the loop of Henle and the gain factor α for the TGF mechanism. Normal values for the feedback delay are about T = 15 s and, as defined in the model, the feedback gain factor is typically α = 10–15 for normotensive rats and some 30–50% larger for spontaneously hypertensive rats. However, as explained above, other parameters also differ between the two strains of rats and it is, therefore, of interest to consider higher values of α where excitation of the myogenic oscillations becomes stronger.
The horizontal curve H at α ≈ 10 represents the Hopf bifurcation at which the onset of the TGFmediated oscillations takes place. Below this curve, the nephron autoregulation displays a stable equilibrium point, and TGFmediated oscillations do not occur. Experimentally, nephrons in normotensive rats are found to operate above, but relatively close to this threshold in most cases. As previously mentioned, the myogenic oscillations are assumed to be damped, but to be continuously excited by the TGFmediated oscillations. The regions delineated by the dotted curves and denoted as 1 : 1, 1 : 2, 1 : 3, 1 : 4 and 1 : 5 represent the different regions of synchronization between the myogenic and the TGFmediated oscillations. For small values of α, the myogenic oscillations lock directly into a 1 : 1 mode as one expects for a periodically driven oscillator at low amplitudes. As the gain factor increases, the excitation of the myogenic oscillator becomes stronger, and a shift to higher and higher excitation frequencies occurs. Doublewavelet analysis of the experimental data has confirmed that the most common frequency ratios under realistic conditions are 1 : 4 and 1 : 5 [4], i.e. the myogenic mode completes four or five oscillations each time the slower TGFmediated mode completes one full period.
The fully drawn curves represent perioddoubling bifurcations. In the middle of the figure, for instance, we observe two cascades of perioddoubling bifurcations denoted PD_{1}^{1:5}, PD_{2a}^{1:5}, PD_{2b}^{1:5}, etc., for the 1 : 5 modelocked solution. In a perioddoubling bifurcation, an existing stable periodic motion loses its stability and is replaced by a stable motion that oscillates alternatingly between two different amplitudes (see, e.g. the experimental pressure trace in figure 1). At the curve PD_{1}^{1:5}, for instance, the synchronized 1 : 5 mode is replaced by a 2 : 10 mode in which the oscillation repeats itself only after two full TGF oscillations and 10 myogenic oscillations. This type of transition is common for nonlinear dynamic systems and often leads to a complete cascade of perioddoubling bifurcations to deterministic chaos. We notice that modelocking is maintained between the two oscillatory components during the perioddoubling process.
A little to the left, the 1 : 4 modelocked solution is found to undergo similar cascades of perioddoubling bifurcations. The dashed curve denoted SN_{1}^{4–5} is a saddlenode (or fold) bifurcation curve that delineates the region of coexistence between the 1 : 4 and 1 : 5 modelocked solutions. In this region, depending on the initial conditions, the model will approach a dynamics with either four myogenic oscillations per TGF oscillation or five myogenic oscillations per full cycle of the TGF oscillation. A saddlenode bifurcation in general denotes a transition in which a stable periodic motion (the node) appears or disappears through collision with an unstable cycle (the saddle) of the same periodicity. In the present case, the region of 1 : 5 dynamics extends all the way to the lefthand branch of SN_{1}^{4–5}, where the stable 1 : 5 cycle collides with an unstable 1 : 5 cycle. Similarly, the region of 1 : 4 dynamics extends all the way to the righthand branch of SN_{1}^{4–5}.
4. Periodically forced nephron
Let us now start to examine the response of the single nephron model to an externally applied periodic variation in the arterial pressure of the form 4.1where P_{a,0} = 13.3 kPa is the average arterial pressure. A = 0.0075 is the relative amplitude of the sinusoidal pressure variation and ω the angular frequency. This frequency is chosen to be close to the TGF frequency of the unforced system.
To be specific, we have chosen a point of operation O corresponding to a delay T = 16 s for the Henle flow and a gain factor α = 24 in the TGF feedback (figure 4). With these parameters, the TGFmediated oscillation has a period of about 40 s, corresponding to an angular frequency ω_{0} = 0.155 s^{−1}. Moreover, the model operates in a regime where the 1 : 5 modelocked solution for the two internally generated oscillations coexists with the 2 : 8 solution (i.e. with the 1:4 solution after its first period doubling). The forcing amplitudes that we can apply must, obviously, be fairly small both to avoid breaking up the modelocking between the internally generated modes and to reduce the effect of mixing between the coexisting modes.
Figure 5 provides an overview of the twodimensional bifurcation diagram for the forced nephron model using the angular forcing frequency ω and the equilibrium value of the afferent arteriolar resistance R_{a,0} as control parameters. We recall that the instantaneous value of the afferent resistance is the variable through which both the myogenic and the TGFmediated regulations are effectuated.
The first thing to notice when inspecting figure 5 is the different shading we have applied to distinguish among different types of dynamics. White denotes synchronized periodic dynamics, i.e. in the white regions the internally generated modelocked solutions synchronize with the externally applied forcing signal. This synchronization involves either the internally generated 1 : 5 periodic solution (in the middle of the figure) or the internally generated 1 : 4 solution (in the lower right corner of the figure). Hatched light grey areas are regions with ergodic twomode dynamics (quasiperiodicity), while the light grey regions represent phasesynchronized chaos. In the region of ergodic twomode dynamics (quasiperiodicity), the internally generated 1 : 4 (or 1 : 5) mode fails to synchronize with the external forcing signal. The result is a type of motion that never repeats itself but gradually fills out a twodimensional torus surface in phase space. In the region of phasesynchronized chaos, the system displays dense sets of unstable periodic solutions, but in such a manner that all of these solutions are synchronized with the external forcing [11,26,27]. Finally, the dark grey areas to both sides of the 1 : 5 synchronization zone represent nonsynchronous chaos. The curves denoted T are torus bifurcation curves where the periodic solution loses stability as two complex conjugated eigenvalues cross out of the unit circle, and quasiperiodic dynamics results.
The following discussion will focus on the bifurcations observed for the internally synchronized 1 : 5 mode. Hence, we have dropped the superscript 1 : 5 on the bifurcation curves that relate to this mode and retained only the subscripts 1, 2, 4, etc., to denote the periodicity of the bifurcating mode relative to the forcing signal. For the perioddoubling curves, superscript S denotes a period doubling of a stable cycle (node) and superscript U denotes period doubling of an unstable cycle (saddle).
Before we leave the 1 : 4 mode, let us note that the region in figure 5 where this mode synchronizes with the external forcing continues to exist in the part of the region of resonance with the 1 : 5 mode. We also notice that the torus arising in the torus bifurcation T_{1}^{1:4} destroys almost immediately and gives rise to synchronized chaotic dynamics. Besides, by this torus bifurcation, the visible part of the 1 : 4 synchronization regime is primarily delineated by the saddlenode bifurcation curve SN_{1}^{1:4}. At this saddlenode bifurcation, a pair of 1 : 4 saddle and stable node cycles is born. The saddle cycle undergoes a perioddoubling bifurcation at PD_{1}^{U,1:4}, and the node doubles its period when crossing the curve PD_{1}^{S,1:4} in the downward direction. Together, these two curves form a closed perioddoubling curve that is tangent to the edge of the resonance zone on both sides. Below PD_{1}^{S,1:4}, period doubling of the 1 : 4 stable node and saddle solution continues (not shown) in a manner similar to the development described below for the 1 : 5 modes.
To better understand the bifurcation structure for the region in parameter space where the internally generated 1 : 5 modelocked solution synchronizes with the external forcing, we have drawn the simplified bifurcation diagrams reproduced as figure 6a,b. Along its left, respectively, righthand side, the 1 : 5 synchronization regime is primarily delineated by the two branches SN_{1}^{L} and SN_{1}^{R} of the saddlenode bifurcation curve SN_{1}. In the projection shown in the figure, this curve displays an extra fold structure, causing some added complexity to the variation of the stability of the 1 : 5 synchronized mode.
At the boundary of the resonance zone, a pair of saddle and stable node 1 : 5 modelocked solutions arises. If we follow the development along the main diagonal starting in the upper right corner, the stable 1 : 5 solution undergoes its first period doubling at the curve PD_{1}^{S}, and the saddle solution doubles its period at PD_{1}^{U}. The two branches of the perioddoubling curve meet at the points a_{1} and b_{1}, where the curve is tangent to the saddlenode bifurcation curve SN_{1}^{L}. Close to the saddlenode bifurcation curve, the 1 : 5 node and saddle cycles approach one another, and at the very edge of the synchronization regime, the perioddoubling bifurcations of the node and saddle solutions occur simultaneously. The points denoted a_{1} and b_{1} (open circles) are thus socalled codimension2 points in which two Floquet multipliers simultaneously cross the unit circle, one through +1 and one through −1. (The Floquet multipliers measure the rates at which neighbouring trajectories approach or move away from a given periodic orbit. A Floquet multipler that exceeds 1 in numerical value signals instability. A Floquet multipler of +1 indicates a saddlenode bifurcation and a multiplier of −1 indicates a perioddoubling bifurcation.)
To the left of the perioddoubling curves PD_{1}^{S} and PD_{1}^{U}, the system, besides a saddle and a doubly unstable 1 : 5 node solution, displays a pair of saddle and stable node 2 : 10 solutions. The region of existence for the 1 : 5 cycles continues to be delineated by SN_{1}^{L} and SN_{1}^{R}. To bound the range of existence for the two 2 : 10 cycles, a new saddlenode bifurcation curve SN_{2}^{L} is born. This saddlenode bifurcation curve is supported by the perioddoubling curves PD_{2}^{S} and PD_{2}^{U} in points (full circles) close to, but not coinciding with a_{1} and b_{1}. We shall return to the detailed bifurcation structure around these points in §5.
Closer inspection of figure 6a shows that the perioddoubling curves PD_{1}^{S} and PD_{1}^{U} together have four additional points of tangency with the saddlenode bifurcation curve SN_{1}, denoted c_{1}, d_{1}, e_{1} and f_{1}. The presence of these points gives rise to further complexity in the bifurcation structure. However, this complexity is not directly in the focus of the present study. Hence, we shall not discuss the bifurcations here, but only mention that the full analysis includes various torus birth bifurcations in which the number of unstable directions for the involved cycles is adjusted up or down by two. These processes are similar in several respects to the processes we shall discuss in §5.
As illustrated in figure 6b, the stable 2 : 10 cycle undergoes a new period doubling at PD_{2}^{S}, and the 2 : 10 saddle solution period doubles at PD_{2}^{U}. The two branches of the second perioddoubling curve meet at the points a_{2} and b_{2} on the saddlenode bifurcation curve SN_{2}^{L}, where the node and saddle 2 : 10 cycles simultaneously perioddouble, merge and disappear. We can hereafter follow the repetition of the same process through a complete perioddoubling cascade and via a region of phasesynchronized chaos (figure 5) to the large period3 window (denoted 3 : 15) in the middle of the figure. Each time a pair of saddle and stable node cycles perioddouble, a new saddlenode bifurcation curve arises to delineate the range of existence for the perioddoubled cycles.
To illustrate the results of the above discussion, figure 7a,b shows a couple of onedimensional bifurcation diagrams obtained by scanning across the 1 : 5 synchronization regime along the lines A and B in figure 5. Following the transitions observed in figure 7a from right to left, we first observe the formation of a pair of saddle (dashed curve) and stable node (solid curve) 1 : 5 cycles from the ergodic torus that exists to the right of the saddlenode bifurcation curve SN_{1}^{R}. At PD_{1}^{S}, the stable node cycle undergoes a perioddoubling bifurcation, and at PD_{1}^{U}, the 1 : 5 saddle cycle suffers a similar bifurcation. At SN_{1}^{L}, the two 1 : 5 modes (now a saddle and a doubly unstable node) merge and disappear. The 2 : 10 saddle and stable node cycles continue to the saddlenode bifurcation SN_{2}^{L} where they merge and disappear. It is interesting to note that the ergodic torus formed in this process is a perioddoubled torus [28,29]. As discussed below, we generally expect the merger and disappearance of a pair of period2^{n} (n = 0,1,2 …) saddle and stable node cycles to produce an ergodic torus that has the same periodicity. In many cases, the range of existence for these perioddoubled tori is very narrow. They tend to disintegrate into nonsynchronous chaos a little away from the resonance zone.
The processes observed along scan line B in figure 7b initially follow the same pattern as those along scan line A. After period doubling of the stable 1 : 5 node at PD_{1}^{S}, there is a small interval with ergodic torus dynamics (quasiperiodicity) before the stable period 2 : 10 node proceeds through a complete perioddoubling cascade to phasesynchronized chaos. The resonant 1 : 5 saddle cycle undergoes its first period doubling at PD_{1}^{U} and thereafter completes a perioddoubling cascade leading to a dense set of doubly unstable nodes. In the mean time, the doubly unstable 1 : 5 node and the 1 : 5 saddle cycle produced in the perioddoubling processes PD_{1}^{U} and PD_{1}^{S} merge and disappear in the saddlenode bifurcation SN_{1}^{L}, and the doubly unstable 2 : 10 node and the 2 : 10 saddle cycle produced in PD_{2}^{U} and PD_{2}^{S} have merged and disappeared in the saddlenode bifurcation SN_{2}^{L}. Similarly, as we approach the left edge of the synchronization regime, many of the other solutions produced in the interconnected perioddoubling cascades gradually merge and disappear, leading the system into a state of nonsynchronous chaos.
Figure 8 illustrates the phenomenon of torus doubling [28,29]. Here we have plotted simultaneous values of the Henle flow and the proximal tubular pressure as determined from a sequence of Poincaré sections of the ergodic torus that exists along the outside edge of the resonance zone a little to the left of the point b_{1} where the first perioddoubling curve PD_{1} is tangent to the saddlenode bifurcation curve SN_{1}^{L}. The equilibrium value of the afferent arteriolar flow resistance R_{a,0} is used as a bifurcation parameter, and the angular frequency of the forcing signal is kept constant at ω = 0.1505 s^{−1}. To the left in the figure, i.e. for R_{a,0} < 2.44 kPa (nl s^{−1})^{−1}, we observe the crosssection of a normal ergodic torus with points distributed along a closed invariant curve. As we increase the arterial resistance R_{a,0}, the return points of the trajectory alternatingly follow two different loops of a perioddoubled invariant curve.
5. Details of the bifurcation structure
The bifurcation diagram in figure 5 includes a considerable number of additional details, some of which are generic for the Ctype perioddoubling transition while others are specific for the considered model. Among the generic aspect is the fact that the new saddlenode bifurcation curves that emerge after a perioddoubling bifurcation of a pair of saddle and stable node cycles emanate in points of the perioddoubling curves close to, but not in the points of tangency with the previous saddlenode bifurcation curves. To illustrate this, figure 9 shows an enlargement of a part of the bifurcation diagram around the point b_{2} in figure 6b. We notice how the saddlenode bifurcation curve SN_{4}^{L} emerges from a point Q on the unstable branch of the perioddoubling curve PD_{2} different from the point of tangency b_{2} with SN_{2}^{L}. A little above Q we observe another point (closed circle). This is the point in which the torus bifurcation curve T_{4} that bridges between SN_{2}^{L} and SN_{4}^{L} is supported by SN_{4}^{L}.
As described above, the saddlenode bifurcation curve SN_{2}^{L} delineates the range of existence for cycles with 2:10 modelocking. Below b_{2}, a pair of 2:10 saddle and stable node cycles are born as the system crosses through SN_{2}^{L} from left to right. Above b_{2}, the pair of 2:10 saddle and doubly unstable node cycles arising from the perioddoubling bifurcations PD_{2}^{S} and PD_{2}^{U} merge and disappear as the system passes through SN_{2}^{L} from right to left. SN_{4}^{L} similarly delineates the range of existence for solutions with 4:20 modelocking. The gap between the two saddlenode bifurcation curves would allow 4:20 modelocked solutions to escape from the synchronization regime had it not been for the presence of the supercritical torus bifurcation curve T_{4}. At this curve, the stable 4 : 20 cycle (now of focus type) looses its stability as two complex conjugated eigenvalues cross out of the unit circle and a period4 ergodic torus is formed. The doubly unstable 4 : 20 cycle (now again a node) subsequently undergoes a reverse perioddoubling bifurcation at PD_{2}^{U}, and the doubly unstable 2 : 10 node produced disappears in the saddlenode bifurcation SN_{2}^{L}.
Similar phenomena take place in the neighbourhood of other codimension2 points. In some cases, the torus bifurcation is subcritical and accompanied by a torus fold bifurcation. The saddlenode bifurcation curve SN_{4} that delineates the region of stable period4 dynamics arises from a point on the stable branch of PD_{2}. The saddlenode bifurcation curve SN_{8} hereafter emerges from a point on the unstable branch of PD_{4}, and the saddlenode bifurcation curve SN_{16} emerges from a point on the stable branch of PD_{8}. In general, we find that the point of emergence for a new saddlenode bifurcation curve tends to alternate through the perioddoubling cascade between the stable and unstable branches of the perioddoubling curves. Moreover, if the saddlenode bifurcation curve in one side of the perioddoubling curves emanates from the stable branch, the saddlenode bifurcation curve at the other side emanates from the unstable branch. As the cascade of perioddoubling bifurcations unfold, more and more saddlenode bifurcation curves accumulate along the edge of the resonance zone. Except for the first few saddlenode bifurcation curves, this accumulation process takes place in an alternating manner and leads to a limiting curve that delineates the existence of phasesynchronized chaos.
6. Discussion
The formulation of dynamic and mechanismbased descriptions of living systems represent an enormous challenge, not only to the biological sciences, but also to mathematics, physics and computer science. Part of this challenge arises from the fantastic number of different processes and feedback mechanisms involved in the regulation of the biological functions and the extraordinary ranges of time and space over which this regulation is maintained. Another contribution derives from the heterogeneity of biological tissues and structures. In the years to come, the biological sciences are likely to become an important test bed for the development of new advanced concepts and approaches in mathematics as well as in physics.
However, as emphasized in §1, an additional source of complexity derives from the fact that living systems operate under farfromequilibrium conditions. This implies that many feedback regulations are unstable and produce selfsustained oscillatory dynamics. Nonlinearity in the feedback and interaction between the oscillatory processes give rise to even more complicated dynamics in the form, for instance, of perioddoubling transitions or transitions between different forms of resonance dynamics. The cells and organs make use of these complex phenomena both to regulate their internal processes and as a means to communicate with other cells and organs.
Nephron autoregulation can be viewed as a set of mechanisms that serve to protect the delicate processes in the distal tubule from fluctuations in the arterial pressure. Experimental studies applying broadband random perturbations to the arterial pressure have shown [15] that the regulation works as a highpass filter that provides damping for arterial pressure variations slower than the response time of the TGF regulation. However, experimental studies also show that the two mechanisms involved in nephron autoregulation—the myogenic and the TGF mechanism—both produce oscillatory dynamics and that coupling between these oscillations and nonlinearity in the systems produce both modelocking and perioddoubling transitions.
Recent theoretical analyses have shown that the perioddoubling transition to chaos along the edge of a resonance zone displays an unusual organization and scaling behaviour, denoted as cyclic or Ctype criticality [30]. It is also known that forced perioddoubling systems may be associated with the appearance of perioddoubled ergodic tori along the outside edge of the resonance zone. This phenomenon was first observed by Arnéodo et al. [28] and by Kaneko [29] who found torus doubling both for multidimensional maps and for timecontinuous systems and described, for instance, how a finite number of torusdoubling bifurcations can lead to chaotic dynamics. More recently, Sekikawa et al. [31] have demonstrated the transition to chaos through a series of subsequent torusdoubling bifurcations in an electronic circuit, and Sekikawa et al. [32] have illustrated the occurrence of torusdoubling bifurcations in coupled logistic and sinecircle maps.
Kuznetsov and coworkers [30,33,34] determined the scaling relations for the perioddoubling process along the edge of an Arnold tongue, and Kuznetsov [35] examined a number of related problems, including the effect of noise on the dynamics near the torusdoubling terminal point for a quadratic map subjected to quasiperiodic forcing. To our knowledge, however, no realistic example of this type of behaviour has so far been described, and the associated bifurcation structure still remains largely unexplored.
Using a simplified nephron model, we have performed a detailed analysis of the response of the regulatory system to variations of the arterial pressure with periods close to the period of the TGFmediated oscillations. To our knowledge, this analysis also represents the first study of a Ctype perioddoubling transition in a realistic system. This transition, which is generic for forced perioddoubling systems [30], is characterized both by a particular bifurcation structure and by specific scaling relations.
The bifurcation structure involves a cascade of simultaneously occurring saddlenode and perioddoubling bifurcations on the edge of the synchronization zone. At the same time, one can observe the process of torus doubling along the outside boundary of this zone. Each pair of perioddoubling bifurcations of the resonant saddle and stable node solutions leads to the formation of a new pair of saddlenode bifurcations along the zone edge. These saddlenode bifurcation curves arise from points on the perioddoubling curves away from the points of tangency with the previously formed saddlenode bifurcation curves. This leaves a gap in the zone edge, a gap that is closed by either a sub or a supercritical torus bifurcation curve. In the case of a subcritical torus bifurcation, this will in general be accompanied by a global (torus fold) bifurcation. In this way, the stable resonance node and the corresponding saddle cycle always communicate with an ergodic torus of the same periodicity across the edge of the resonance zone.
As the perioddoubling cascade unfolds, more and more saddlenode bifurcations accumulate along the zone boundary. In accordance with the scaling theory [30], our detailed simulations indicate that this accumulation takes place in an alternating manner and in the limit leads to a welldefined curve that delineates the region of phasesynchronized chaos from that of nonsynchronous chaos.
The presence of noise in the biological system will obviously wash out some of the finer details in the bifurcation structure. However, doublewavelet analysis and other effective statistical methods have allowed us to clearly demonstrate the existence both of episodes of perioddoubling behaviour and of intranephron synchronization in the temporal variation of the proximal tubular pressure [5]. Transitions between different states of dynamics may represent elements in the normal physiological regulation or they may signal the development of a disease. Some types of hypertensive rats, for instance, are clearly distinct in their tubular pressure dynamics as well as in their intra and internephron synchronization behaviours.
Acknowledgments
We thank A. N. Pavlov for performing the wavelet analyses of the tubular pressure variations. The project was sponsored by the European Union through the Network of Excellence BioSim (contract no. LSHBCT2004005137).
Footnotes

One contribution to a Theme Issue ‘Advancing systems medicine and therapeutics through biosimulation’.
 Received September 3, 2010.
 Accepted November 10, 2010.
 This Journal is © 2010 The Royal Society