Part II: Mathematical Modelling of
Progressive Renal Disease

So far I have discussed a model motivated by renal and patient survival data. The discrete event simulation uses as its motor three survivor functions: the patient waiting list survivor function, the transplanted patient survivor function and the graft survivor function. These functions, along with a random number generator and copy of the original 434 patient records allowed us to test a variety of what if? scenarios. This type of model is of the Monte Carlo type which are widely employed in the mathematical modelling literature [see, for a discussion, the book by Press et al (1993)]. In essence, what we are doing is realising numerous real world scenarios. Our only assumption is that our database of 434 records captures the essential features of the true (yet unobtainable) distribution of renal graft patients in the West of Scotland.

I'd like to present you with another type of model which is one stage more abstract than the one discussed above. We will use differential equations to describe a mathematical model which is not driven by data, as in our discrete event simulation, but which founds itself on a set of fundamental underlying laws.

Let's draw an analogy with fluid flow. Continuum physical theories suggest that force balances and conservation laws must be obeyed by any fluid be it water, blood or air. Consider the entire atmosphere for example. The atmosphere could be completely described if we knew the wind speed and direction at every position above sea-level. But how can we obtain this? We find that the wind velocity is embedded in a force balance between the inertial, pressure gradient, buoyancy and viscous forces. This isn't yet enough to fully extract the wind speed at every point, but if we ensure that the air mass is conserved we are home and dry.

The theory of fluid dynamics has been much developed over the last one hundred or so years, and the hall mark of the subject must surely be the famous Navier-Stokes equations from which we can deduce the flow speed and direction. The equations are not an end in themselves as many situations arise where they are simply too difficult to solve. Mathematical Biology (and more recently, Mathematical Medicine) is a very new subject. Very often, mathematical biologists are, or have had a training in, subjects such as fluid dynamics. However, there is often no coherent mathematical theory for biological processes analogous to that found for fluids. This is certainly the case in renal physiology although some models have been produced. Generally these models are for single nephrons and consist of systems of partial differential equations which do not look so different from the equations for fluid flow! We are interested in renal disease here, but for a description of some of these physiological models (whose future inclusion in disease models is not precluded), see Chapter 20 of the book by Keener & Sneyd (1998), Layton et al (1995) or the special edition of the Bulletin of Mathematical Biology (1994) which was dedicated to renal physiology.

Unfortunately, very little mathematical modelling has focussed on the nature of progressive renal disease. There are many good reviews in the medical literature, e.g.that by Remuzzi & Bertani (1998) or from the book edited by Mitch et al (1986) where there is agreement that the progession to chronic renal failure is independent of the nature of the initial insult. From a simplistic point of view, serving our purposes, any insult leads to nephron loss which induces the release of growth promoting factors. These factors induce glomerular hypertrophy and hence glomerularsclerosis leading to more nephron loss and more hypertrophy. There is hence a runaway snowball effect until chronic renal failure occurs. A strict low protein diet will minimise increased capillary pressure and flow and slow the disease process. But how long can you hold off the devil? In some instances, intervention can be very successful at slowing down the progression of the disease if it is diagnosed sufficiently early.

One mathematical model of chronic renal disease is that by Chaturvedi & Insana (1997) [CI97]. Their work models the progression of renal disease in an experimental subtotal nephrectomy rat model [Yoshida et al (1989)]. CI97 apply a simple dynamical systems approach and construct a mathematical model of the disease process in terms of the sclerosis index s [Raij et al (1984)], the glomerular diameter g and in terms of the number of remaining non-diseased nephrons n. They consider these quantities only as they vary with time t. By dynamical system they refer to the coupled system of nonlinear ordinary differential equations:


$${\mbox{d} {s}\over \mbox{d}t} = \alpha g,$$ (1)

$${\mbox{d} {g}\over \mbox{d}t} = \Sigma - \beta sg,$$ (2)

$${\mbox{d} {n}\over \mbox{d}t} = -\lambda n.$$ (3)



Here $\Sigma$ is the constant rate of glomerular hypertrophy, $\beta$ is the constant rate at which sclerosis reduces glomerular size and $\lambda$ is the constant disease rate of the nephrons.

An important point to note is that the equation for $n$ in (3) is actually decoupled from the other two. Given a set of initial conditions:


$$s(0)=s_0,\ g(0)=g_0,\ n(0)=n_0,\quad (\mbox{for constants }s_0,g_0,n_0>0);$$ (4)



we see easily that $n(t)=n_0\exp(-\lambda t)$. It thus appears that the number of non-diseased nephrons has nothing at all to do with sclerosis nor hypertrophy! This is not so. CI97 consider an initial population of $n_0=50$ nephrons. Each nephron is tracked through its individual disease process by equations (1) and (2) when and only when it becomes diseased. This process is governed by (3). For example, when $n$ drops by one from 23 to 22, say, nephron number 23 enters the disease process. It is given a random set of initial conditions $s_0$ and $g_0$ selected au hasard from exponential and gaussian distributions (respectively) and tracked until the sclerosis index reaches its maximum sclerosed state. CI97 can then plot $s$ versus $g$ on a graph at equal time intervals after the beginning of the experiment.

CI97 have enjoyed some success with their model. They manage to capture the instantaneous parabolic relationship between hypertrophy and sclerosis projected by Yoshida et al (1989). However, the model is quite unsatisfactory in many ways. There are two independent disease mechanisms at work: the first is through equation (3) and the second through the interplay between $s$ and $g$. Surely, $s$, $g$ and $n$ should all interact simultaneously? Why not simply do without equation (3) if it is indeed superfluous? CI97 can't! If CI97 simply placed their whole population of nephrons on the plot $(s,g)$ at time $t=0$ all nephrons would become immediately diseased and proceed to total sclerosis regardless of the initial conditions! Even for a healthy kidney where all the nephrons have initial conditions of $s(0)=0$ the model must necessarily lead to chronic renal failure! Things become even more bizarre as we reach the maximum sclerosed state of $s=4$. There, sclerosis continues to increase! In fact, $s\to\infty$ as $t\to\infty$.

What we are in fact alluding to is that there is no equilibrium configuration for nephrons in the model of CI97. Equilibrium or steady state configurations occur when $s$, $g$ and $n$ cease to vary with time $t$. There, the derivatives must necessarily be zero. This leads us to a set of equilibrium equations for (1), (2) and (3):


$$0 = \alpha g,$$ (5)

$$0 = \Sigma - \beta sg,$$ (6)

$$0 = -\lambda n$$ (7)



and there can only be a solution to these equations if $\Sigma = 0$. I am sure that CI97 did not intend for this to be the case! What we require for each nephron is a healthy state and a diseased state. Healthy nephrons should stay healthy unless they are perturbed, in which case we know that they should proceed along a disease path leading to complete sclerosis. We need a new dynamical system that can accommodate two such states. I will now discuss the model which I have developed.

Let's assume that the sclerosis index $s$ is the direct measure of the progressive disease process in a single nephron. Suppose that any nephron can be in one of two equilibrium states: completely healthy ($s=0$) and completely sclerosed ($s=s_{max}$). Then, if, for a particular nephron, $0<s<s_{max}$ is true, then the nephron must be in a transitory state moving between equilibrium states. We know that once sclerosis has begun we may slow it down with a low-protein diet, for example, so let us assume that once perturbed from a healthy equilibrium state, the only path is toward a completely sclerosed state with $s=s_{max}$. A good mathematical model of this scenario is the Verlhust-Pearl equation which has been used extensively as a nonlinear population model. However, here its usefulness is in its two equilibrium states. Let's model the sclerosis index as:


$${\mbox{d} {s}\over \mbox{d}t} = s\left(1-{s \over s_{max}}\right).$$ (8)



Here, we see that the equilibrium states are $s=0$ and $s=s_{max}$ by setting the left hand side of (8) to zero. Linearisation around the two equilibrium states shows that $s=0$ is an unstable state and $s=s_{max}$ a stable state as we require.

We must now build a model for the glomerular diameter's response to the disease process modelled in $s$. For simplicity, let's redefine the variable $g$ in the formula:


$$G(t) = g_0+g(t).$$ (9)



Here, $G$ is the actual glomerular diameter and $g$ measures the deviation in size from the standard healthy diameter of $g_0$. Then a value of $g=0$ represents a healthy nephron. Yoshida et al (1989) have shown that during the disease process, nephrons initially increase their diameter up to some critical width before shrinking back down. For $g$ then, we wish an unstable healthy state of $g=0$ and $s=0$ leading, after a small perturbation, to an increase and then a decrease of $g$ through to a final stable sclerosed state at $g=0$ and $s=s_{max}$. A variation on the Verhulst-Pearl equation will work well:


$${\mbox{d} {g}\over \mbox{d}t} = g\left(1-{g \over g_{max}}\right)-{2sg \over s_{max}}.$$ (10)



What does (10) mean? It behaves in much the same way as the equation for sclerosis (8) until $s$ becomes large. It is then the action of the sclerosis on the increasing hypertrophy that must eventually reduce $g$. If we set the left hand side of (10) to zero we see that either $g=0$ is a solution or $g=g_{max}(1-s/s_{max})$. However, in the latter case, we know that $s=s_{max}$ yields a stable sclerosed state, hence $g=0$ is the value to which this corresponds.

Thus equations (8) and (10) constitute our single nephron disease process model. Any perturbation from the healthy state (an insult) will send the nephron on a disease path to the sclerosed state as we require. There are in fact two other equilibrium solutions to our system at $g=-g_{max}$, $s=s_{max}$ and at $g=g_{max}$ and $s=0$. These are what we call saddle points and we will say some more about them later.

Now that we have our single nephron disease model, how should we make it work for a population of nephrons or for a whole kidney? We need to bring in the idea of the glomerular filtration rate (GFR). The GFR is proportional to the total available surface area presented by all the nephrons. Let's call this area $A(t)$ and further suppose that there are initially a population of $n_0$ nephrons in a kidney. Prior to insult, the total glomerular surface area is $n_0g_0^2$. However, after an insult (e.g.a subtotal nephrectomy) the new available surface area will be:


$$A(t) = \sum_{j=1}^{n(t)}[g_0+g_j(t)]^2-n_0g_0^2.$$ (11)



where a subscript $j$ simply indexes an individual nephron. As time progresses in our model simulation, if we ensure that at time $t$ a number


$$n_0 \left( 1-{A(t) \over n_0g_0^2} \right)$$ (12)



of nephrons are given random but small perturbations $s>0$ and $g>0$ to bring them onto a disease path, then at least initially the hypertrophy will stop the available surface area from reducing and hence maintain a ``healthy'' GFR. Eventually, however, most (or all) of the nephrons will have become sclerosed and there will be nothing left to stop the GFR from tumbling toward zero and ultimate renal failure.

Figure 5: A variety of single nephron disease paths that result from the following initial perturbations in a plot of sclerosis $s$ versus hypertrophy $g$: ``A'' corresponds to a perturbation $s=.001$, $g=.001$; ``B'' to $s=.001$, $g=.01$; and ``C'' to $s=.01$, $g=.001$.



In terms of what perturbations to give nephrons, I have illustrated a variety of single nephron disease paths in a plot of $s$ versus $g$ in Figure 5 where I have taken $s_{max}$ and $g_{max}$ both equal to 1. Note that how, depending on the random orientation of perturbation, nephrons can become quite hypertrophied before succumbing to sclerosis or only lightly hypertrophied. Here I suggest that nephrons be selected in a random order onto the disease path and that the degree of the initial perturbation be of the magnitude of $10^{-2}$ but in an random orientation with the $s-$ and $g-$axes.

One unsatisfactory point of this model is that the orientation of the perturbations should not have $s=0$. Then, the model will simply increase the nephron's diameter to $g=g_{max}$ without increasing $s$. This is due to a saddle-type equilibrium point which exists for this at the point $s=0$ and $g=g_{max}$. All we need do to avoid this problem is to ensure that positive perturbations (no matter how small) are given to both $s$ and $g$. Finally, one other saddle point exists for my system at the point $s=s_{max}$ and $g=-g_{max}$. However, for positive perturbations, negative values of $g$ or $s$ can never be reached.