The Effect of Different Forms for the Delay in a Model of the Nephron

We investigate how the dynamics of a mathematical model of a nephron depend on the precise form of the delay in the tubuloglomerular feedback loop. Although qualitative behavioral similarities emerge for different orders of delay, we find that significant quantitative differences occur. Without more knowledge of the form of the delay, this places restrictions on how reasonable it is to expect close quantitative agreement between the mathematical model and experimental data. 1. Introduction. Delay arises naturally in many mathematical models of biological systems, including population dynamics [7], epidemiology [2], and physiology [6]. In some cases the precise form of the delay is well defined, for example, if it represents the maturation age of a species. However, in many cases, the exact form that the delay takes is less easy to measure. In this article, we investigate a mathematical model that describes the dynamical behavior of a nephron, the main functioning unit in a kidney. Particular emphasis is placed on assessing the impact on the dynamical behavior of the exact way in which the delay is modelled. There are around one million nephrons in each human kidney. Blood arrives at each nephron via the afferent arteriole. At the nephron, the blood is filtered by the glomerulus and then leaves via the efferent arteriole. The filtrate continues into the tubes of the nephron, entering first the proximal tubule and then the loop of Henle before leaving the nephron via the distal tubule to enter the collecting duct. In the collecting duct, final adjustments to the concentration of the filtrate are made before it leaves the kidney as urine and passes to the bladder. Feedback between the fluid at the end of the loop of Henle and the radius of the afferent arteriole governing the amount of fluid entering the nephron occurs via the macula densa cells. This is known as the tubuloglomerular feedback (TGF). This feedback mechanism ensures that the pressure in the nephron is kept relatively constant and within the operating range of the nephron. Delay occurs in the TGF mechanism for two main reasons: first the time fluid takes to pass around the loop of Henle and, second, the time taken for a signal to be passed via macula densa cells from the end of the loop of Henle to the afferent arteriole [9]. Experimental observations in rats indicate that, in normotensive rats, the pressure in the proximal tubules either …


1.
Introduction. Delay arises naturally in many mathematical models of biological systems, including population dynamics [7], epidemiology [2], and physiology [6]. In some cases the precise form of the delay is well defined, for example, if it represents the maturation age of a species. However, in many cases, the exact form that the delay takes is less easy to measure. In this article, we investigate a mathematical model that describes the dynamical behavior of a nephron, the main functioning unit in a kidney. Particular emphasis is placed on assessing the impact on the dynamical behavior of the exact way in which the delay is modelled.
There are around one million nephrons in each human kidney. Blood arrives at each nephron via the afferent arteriole. At the nephron, the blood is filtered by the glomerulus and then leaves via the efferent arteriole. The filtrate continues into the tubes of the nephron, entering first the proximal tubule and then the loop of Henle before leaving the nephron via the distal tubule to enter the collecting duct. In the collecting duct, final adjustments to the concentration of the filtrate are made before it leaves the kidney as urine and passes to the bladder. Feedback between the fluid at the end of the loop of Henle and the radius of the afferent arteriole governing the amount of fluid entering the nephron occurs via the macula densa cells. This is known as the tubuloglomerular feedback (TGF). This feedback mechanism ensures that the pressure in the nephron is kept relatively constant and within the operating range of the nephron. Delay occurs in the TGF mechanism for two main reasons: first the time fluid takes to pass around the loop of Henle and, second, the time taken for a signal to be passed via macula densa cells from the end of the loop of Henle to the afferent arteriole [9].
Experimental observations in rats indicate that, in normotensive rats, the pressure in the proximal tubules either can be steady or can undergo regular oscillations, but for hypertensive rats the pressure can be chaotic [12,8]. Detailed mathematical models containing the key physiological elements of the TGF mechanism have sought to demonstrate that for realistic parameter values, steady behavior, regular oscillations, and chaotic phenomena are possible [11,9,1]. Coupling between the TGF mechanism and myogenic mechanisms have also been investigated [13], and more recently, the coupling between neighboring nephrons have been considered both experimentally [14] and through a mathematical model [10]. In all the mathematical models, a third-order distributed delay is used to model the delay as simpler, although less destabilizing, than a discrete delay [9].
In this paper, we examine how the dynamics of the mathematical model presented in [1] vary depending on the order of the delay, taking as an extreme case a discrete delay. We investigate the extent to which quantitative agreement with experimental results can be expected without more detailed modelling of the delay. In section 2, we give a brief description of the mathematical model for the nephron and the types of delay terms we have investigated. In section 3, the results are presented followed, in section 4 by a discussion of the key points.
2. Mathematical model. The physiological basis for the mathematical model of the nephron that we use is discussed in detail in [1], and only a brief overview is given here. The main component of the model is an equation for the time evolution of the pressure in the proximal tubule that results from considering the filtrate that is flowing in from the glomerulus, out to the loop of Henle, and out through the walls of the proximal tubule. This equation for the pressure in the proximal tubule is coupled to an equation for the radius of the afferent arteriole in two ways: first, the radius of the afferent arteriole has a direct bearing on the amount of fluid that flows into the proximal tubule, and second, because the fluid flowing out of the proximal tubule is directly related to the concentration of chloride ions at the end of the loop of Henle. It is the concentration of chloride ions detected by the macula densa cells that determines the radius of the afferent arteriole.
Specifically, the pressure in the proximal tubule P t is modelled by Here, C tub is the elasticity of the proximal tubule, F filt is the filtrate that enters the proximal tubule from the afferent arteriole, F reab is the amount of fluid that leaves through the walls of the proximal tubule and enters the interstitial space between the nephrons, and F hen is the fluid that flows out of the proximal tubule and into the loop of Henle. Typical values for C tub and F filt are given in Table 1. F hen depends on the pressure in the proximal tubule P t along with various measured parameters. The precise dependence is given in Equation (6), and measured values for the parameters are given in Table 1. The amount of incoming fluid F filt depends on both P t and the radius of the flexible part of the afferent arteriole, r, and is given in Equation (7). For further details of the physiological modelling, the reader is referred to [1]. Experiments have indicated that the flexible part of the afferent arteriole behaves like a second-order differential equation where d is a measured damping coefficient, ω is an effective mass to elasticity ratio for the arteriole wall, P av is the average pressure in the arteriole, and P eq is the equilibrium pressure of the arteriole wall. Typical experimental values for d and ω are given in Table 1. The average pressure P av depends on the pressure in the proximal tubule P t in a non trivial way and is given in Equation (11). The equilibrium pressure P eq of the arteriole wall depends not only on the radius of the afferent arteriole r but also on the pressures applied by the macula densa cells. In turn, this depends on the concentration of chloride in the fluid leaving the loop of Henle. An expression for P eq is given in Equation (12). Letting V r = dr/dt, the full equations for the nephron model are three first-order ordinary differential equations, along with ten further relationships that contain the precise details of the functional relationships for the different terms in the equations. The three first-order equations are with the nine auxiliary equations P eq = 1.6(r − 1) + 0.006 exp(10(r − 0.8)) where C e is given by solution of the cubic equation In deriving these equations, careful modelling of all the known physiological processes was made with values for the different constants taken from experiments. The delay is included in the equation for Ψ through F delay . The function Ψ itself describes how sensitive the equilibrium pressure P eq in the afferent arteriole is to the concentration of chloride at the macula densa cells. For the nephron, direct evidence for the precise form of the delay F delay is difficult to come by. Typical transit times for the fluid to pass around the loop of Henle are quoted in [9] as 2 to 3s, with the total time for a change to glomerular capillary pressure after changes in the tubular flows being typically 4 to 5s. In work on coupled nephrons [10], typical delay times of 16s are quoted. If it is assumed that it always takes precisely the same time for the afferent arteriole to respond because of changes to fluid leaving the proximal tubule, then the most appropriate way to model the delay is as a discrete delay, where T is the delay time. However, it is likely that the delay time is not always precisely the same, but has some distribution around a mean non-zero value. In this case it is appropriate to take a distributed delay, A common choice for the kernel f (τ ) is to use the form where θ n = n/T and α n = 1/(n − 1)!. The values of α n and θ n come from setting the mean of each distribution to T and normalizing each to one. The value of n is known as the order of the delay. The kernels for n = 1, 2, 3, 5, and 10 are plotted in Figure 1 for T = 1. The case n = 1 is commonly known as a weak delay and takes its maximum value at τ = 0. In contrast, the higher-order delays have a kernel that approaches zero as τ → 0 and a maximum occurring at non zero τ . The higher the order of the delay, the more strongly the distribution of delay times is located around the mean value. The discrete delay (16) can be formally considered as the limit of the distributed delay (17) as n → ∞, while keeping the mean delay T fixed. Mathematical models for the nephron, such as that discussed by [1], have focused on the third-order delay, n = 3. Note that this still gives significant weighting to very short delay times and that one would expect a higher-order delay to give a better approximation. Our purpose here is to investigate to what extent it matters what order delay is taken. Note that other choices for the kernel can be made, but the main advantage of using those of the form (18) is that the "linear chain trick" can be used to transform the three integro-differential equations (3)-(5) of the nephron model to three ordinary differential equations supplemented by n first-order differential equations for the delay. For example, for the first-order delay n = 1, we have and then differentiating with respect to t leads to where F hen is a function of P t and is given by Equation (6). Hence, for the firstorder delay, the model for the nephron is represented by Equations (3)-(5) coupled to Equation (19) along with the auxiliary Equations (7)-(15). In a similar way, the nth order delay results in supplementing Equations (3)-(5) and (7)-(15) with n first order ordinary differential equations In the following section, we present the results of a numerical investigation for six different cases, five resulting from taking different-order kernels of the form (18) and one from considering the discrete delay (16).

Results.
In each different delay case, the dynamics of the nephron model were explored as a function of the two parameters α and the mean delay T . The value of α, the sensitivity of the macula densa cells to changes in chloride ions, is expected to increase as rats change from normo-to hypertensive. The same values as given in [1] were taken for the other physiological parameters, and these are given in Table  1. We note however that we have taken the numerical values in the expression for P eq from [3].
To compare the different cases, bifurcation sets showing the locus of different types of bifurcations that bound regions of different dynamical behavior as a function of T and α were computed. The complexity of the model equations makes analytic progress difficult, and hence numerical methods were used. In the case of the ordinary differential equation models that result from applying the linear chain trick to the nephron model with the distributed delays of the form (18), results were produced using AUTO [4]. For the discrete delay case, DDE-BIFTOOL [5] was used. AUTO is a well-developed bifurcation and path-following package that has been extensively used to investigate the dynamical behavior of a range of systems. The user supplies the differential equations in a specified format and, by setting the values of various flags, can follow paths of steady-state or periodic solutions in the parameter space. In addition, AUTO finds steady-state, Hopf, saddle-node and period-doubling bifurcation points as a function of user-specified parameters and locates features such as homoclinic orbits. Information for each steady-state solution and any bifurcation point can be saved to provide initial data for subsequent runs. The locus of bifurcation points can be followed in two parameters by starting at a previously computed bifurcation point and setting the flags appropriately. DDE-BIFTOOL is a bifurcation and path-following package for delay differential equations that runs within MATLAB, but is less developed and tested than AUTO. The user supplies the equations in the required MATLAB format and then interactively runs a series of MATLAB routines to follow steadystate solutions, to compute their stability, to find bifurcation points, and to follow the bifurcation points as a function of two parameters.
For all the models of the delay and for sufficiently small values of α, a steady state exists where r and P t do not vary with time. As α is increased, a Hopf bifurcation occurs and stable periodic oscillations in r and P t result. The precise values for the Hopf bifurcation point depend on the particular model of the delay and the mean delay time T ; although as T → 0, all models necessarily behave the same as they become equivalent. When α is increased further, complicated sequences of period-doubling bifurcations leading to chaos can result, depending on the value of T .
The bifurcation sets for the three lowest order kernels are shown in Figures 2, 3  and 6. In the case of Figure 2, the dynamics for the range of parameters investigated is straightforward. For low values of α the steady-state solution with P t and r constant is stable. If α is increased, the steady-state solution becomes unstable and is replaced by a stable periodic solution in which both P t and r oscillate. For the range of α shown and for sufficiently large T , a further increase in α results in a destabilization of this periodic solution, and it is replaced by a period-two oscillation in P t and r. The lowest value of α for which further period-doubling bifurcations occur is α ≈ 17.9 and is therefore off the scale of this figure. Figure 3 shows that the case for the second-order delay is more complicated at low values of α than the first-order, with lines of both period-doubling and saddlenode bifurcations. The period-doubling bifurcations that are plotted represent the first three bifurcations (period one to period two, period two to period four, period four to period eight) in a cascade of such bifurcations. Successive period-doubling bifurcations occur at smaller and smaller intervals in α and, in general, numerical integrations within the regions bounded by the four-eight period-doubling curves indicate complex dynamical behavior with multiple stable solutions and regions of chaos. Different period-doubling sequences occur from the period two orbit depending on the value of T . For example, on increasing α from T = 2.5, the periodic orbit first bifurcates along the line pd 12 , and then the period two orbit bifurcates along the line pd 24a . However, for T = 4, the periodic orbit bifurcates along the line pd 12 , but now the period two orbit bifurcates along the line pd 24b . These different period-doubling bifurcations are part of sequences of such bifurcations that are folded on top of each other, and the dashed curve "sn" on the bifurcation set indicates the path of the fold (saddle-node) lines. To illustrate this folded structure, a bifurcation diagram showing the maximum value of r for α = 12 and for increasing T is given in Figure 4. The maximum value of r rather than the L2norm of the solution or the maximum value of P t is used as a convenient measure that does not produce spurious intersections in the bifurcation diagram. Stable  solutions are indicated by a solid line and unstable ones, by a dashed line. As T is decreased the steady-state solution loses stability at the Hopf bifurcation at H. This Hopf bifurcation is supercritical, and a branch of stable periodic solutions emerges. This is destabilized at pd 12 , corresponding to one side of the "tongue" produced by the one-two period-doubling bifurcation shown in the bifurcation set in Figure 3. This solution is restabilized at the second point marked pd 12 , corresponding to the other side of the tongue. Between these two points, the folding and the perioddoubling cascades occur: note that we have included only bifurcations up to the two-four period-doubling in this diagram. One consequence of the folded structure is that there are regions where different types of period two orbits coexist stably; for example, at α = 12 in the region between T ≈ 4.62 and T ≈ 4.89. Numerical integrations showing the form of the two different period-doubled solutions for the particular case α = 12 and T = 4.8 are shown in Figure 5. Figure 5a corresponds to the upper branch and Figure 5b corresponds to the lower branch.
The folded structure of period-doubling bifurcations also occurs for the thirdorder delay, as shown in Figure 6 and discussed in [1]. In the bifurcation set for the third-order delay, note again only the lines for the first three bifurcations in the various period-doubling cascades are plotted. Although the same general features occur in the bifurcation sets for both the second-and third-order delays, in the latter case the bifurcations are more squashed and the folding more extreme: we have plotted two lines of saddle-node bifurcations, but more exist.
We have explored the bifurcation sets for the higher-order delays of n = 5 and n = 10, and the same general features are seen. In particular, with the exception of n = 1, there is a general trend that the higher the order of the delay the lower the value of α for which oscillations first occur. This trend can be seen in Figure 7, where the Hopf bifurcation lines for n = 2, 3, 5, and 10 are plotted along with the one for the discrete delay. A similar trend indicating that the higher the order of the delay the earlier the onset of more complex dynamics is seen by comparing the position of the line of the first period-doubling bifurcation for different delays (see Fig. 8).   4. Discussion. The precise form of the delay that occurs in this biological system is hard to measure. It consists of two components: (1) the finite time that it takes for fluid to pass around the loop of Henle and (2) the time for the communication between the macula densa cells at the end of the loop of Henle and the afferent arteriole. Probably the form of the delay is distributed to some extent, with the distribution significantly localized around a finite non zero mean value for the delay. In previous work, a third-order delay has been used and the linear chain trick employed to convert the resulting integro-differential equations to ordinary differential equations. This is simpler than considering a discrete delay. Choosing a third-order delay is a compromise between not adding too many ordinary differential equations to the already complex system of three first-order differential equations describing the behavior of the nephron but at the same time allowing for some degree of localization around a mean.
In this paper, we have shown that taking a distributed delay of the form given in equations (17) and (18) has relatively little qualitative impact on the bifurcation set. However, if quantitative comparison with physiological data is the aim, then changing the order of the delay can have a significant impact on the onset of oscillations and, to a lesser extent, on the onset of more complicated dynamical behavior such as chaos. For the parameter ranges we have explored, the larger the mean delay time T , the more significant these differences are.
Finally, it is frequently stated that the stronger the delay, the greater the propensity for oscillations to occur. We note that this is not true for the model of the nephron, as shown by the comparison of the Hopf bifurcation curves for the cases n = 1 and n = 2 in Figure 9.