## Abstract

In this paper, concepts from network automata are adapted and extended to model complex biological systems. Specifically, systems of *nephrons*, the operational units of the kidney, are modelled and the *dynamics* of such systems are explored. Nephron behaviour can fluctuate widely and, under certain conditions, become chaotic. However, the behaviour of the whole kidney remains remarkably stable and blood solute levels are maintained under a wide range of conditions even when many nephrons are damaged or lost. A network model is used to investigate the stability of systems of nephrons and interactions between nephrons. More sophisticated dynamics are explored including the observed oscillations in single nephron filtration rates and the development of stable ionic and osmotic gradients in the inner medulla which contribute to the countercurrent exchange mechanism. We have used the model to explore the effects of changes in input parameters including hydrostatic and osmotic pressures and concentrations of ions, such as sodium and chloride. The intrinsic nephron control, tubuloglomerular feedback, is included and the effects of coupling between nephrons are explored in two-, eight- and 72-nephron models.

## 1. Introduction

The kidney is responsible for maintaining the stability of the body's extracellular environment. It does this by regulating the solute concentrations in the blood, the volume of extracellular fluid, the acid–base balance and blood pressure in the body, as well as excreting urea and the other metabolic waste products. The kidney also plays a part in the body's endocrine system, secreting hormones such as *renin* as part of the *renin–angiotensin* system, *erythropoietin* to regulate red blood cell production and *thrombopoietin* to regulate platelet production.

The human kidney is made up of approximately 1×10^{6} *nephrons* (Guyton & Hall 2006), convoluted *tubules* that adjust the solute levels of the body's blood plasma. Nephrons are supplied with blood by a network of blood vessels, and filter the blood, removing metabolic wastes and reabsorbing the solutes and hormones needed by the body.

The behaviour of individual nephrons can fluctuate widely and can even become chaotic (Marsh *et al*. 2005), but the behaviour of the kidney remains remarkably stable. Kidneys maintain their homeostatic functions over an astonishing range of bodily conditions, such as wide variations in blood pressure, extracellular fluid volume and solute concentrations. Kidneys are also resilient to the effects of renal diseases and the loss of renal tissue, for example through surgical removal, and a human can survive with only one-half of a single kidney.

To understand just how the kidney achieves this, a number of mathematical models have been advanced. The majority of these, however, have focused on specific behaviours of nephron segments, for example, the study of oscillations in glomerular filtration rate (GFR), transport in the distal and proximal tubules, diuresis in the proximal tubule or the urine-concentrating effects of the loop of Henle. Some uncertainty still remains about the details of these effects and the current models do not capture the essence of all observed phenomena (Weinstein 2003). Fewer models studying multi-nephron systems exist.

This paper introduces a new computational kidney model, one that takes a more holistic approach by modelling nephrons as networks of tubules and larger structures in the kidney as networks of nephrons. The focus is on investigating how the behaviours of individual nephrons combine to achieve autoregulation and other important effects of the kidney.

In particular, we view systems of interacting nephrons as an example of *complex systems* (Shalizi 2006). Kidneys are composed of many interconnected parts whose behaviours are both variable and dependent on the behaviour of the other parts. The dependencies often lead to emergent behaviour, which often has properties that can only be examined at higher levels of hierarchy than the individual parts. The difficulty in understanding such systems is that the behaviour of the whole is highly dependent on the individual parts, and each part is highly dependent on other parts. Network models can help to understand emergent *dynamics* by providing an explicit model of the interactions that occur in a system. They are also well suited to simulation, since network automata are well-known models of computation.

Our network models are capable of reproducing experimentally observed behaviour, and their performance characteristics render the simulation of whole-kidney function feasible. In turn, such a computational model of the whole kidney facilitates some novel research possibilities. An example is to study the effects of kidney disease on system dynamics, and to investigate the stability under perturbations due to disease.

The paper is organized as follows. Section 2 briefly introduces the anatomy of the kidney and surveys existing computational nephron models. In §3, we define our single nephron model (§3*b*), and our multi-nephron models (§3*c*). Specifically, in §3*c*, the single nephron model is used to construct systems of interacting nephrons. In §4, we give a brief analysis of the behaviour of our models. In particular, we explore the response of the single nephron model's filtration rates to changes in pressure and sodium concentration and show that it accords with the observed data. We also explore the stability of the sodium gradient in two-nephron and multi-nephron systems and the oscillations in GFRs under perturbations to the model. In §5, we explore future directions and conclude.

## 2. Renal modelling

### (a) The kidney and nephrons

To begin with, the gross structure of the kidney is made up of an outer portion, the *cortex*, and an inner portion, the *medulla*, which consists of an *outer layer* and an *inner layer*. Blood is supplied to the kidney via the *renal artery* and drains into the *renal vein*. Kidneys have a hierarchical organization. Nephrons can be grouped into interacting neighbourhoods; for example, all of the nephrons joined to a single cortical radial artery may be grouped together. Large numbers of neighbourhoods form *renal pyramids* and the kidney consists of a number of such pyramids.

Nephrons are the basic functional units of the kidney and are long segmented tubes consisting of a glomerulus and a series of tubules (figure 1). Nephrons are divided into a number of tubules based on tubule function and the reabsorption and secretion characteristics of the tubule. With reference to figure 1, blood flows through the afferent arteriole (1) and into the glomerulus (2), where a fraction of the blood plasma is filtered into the nephron tubule, and the remaining blood drains into the efferent arteriole and then into the peritubular capillaries (not shown). The filtrate then flows through the proximal tubule (3), the descending and ascending limbs of the loop of Henle (4–6), the distal tubule (7) and finally into the collecting ducts (8). The loop of Henle acts as a countercurrent exchange mechanism that increases the solute concentration in the interstitial fluid as the loop descends. The gradient between the top and bottom of the loop can be as large as 1200 mmol l^{−1} in human kidneys (Stephenson 1992) and over 2000 mmol l^{−1} in other species (Hervy & Thomas 2003).

The kidney maintains a relatively constant renal plasma flow (RPF) and GFR through *autoregulation*. A key mechanism in the nephron contributing to autoregulation is *tubuloglomerular feedback* (TGF). TGF is a negative feedback mechanism that regulates the filtration rate of an individual nephron based on the delivery of solute to the end of the ascending limb of Henle. The top of the ascending limb of the loop of Henle loops back to the nephron's glomerulus where lies the *juxtamedullary apparatus* that contains the *macula densa* (located at (9) in figure 1), a collection of cells that can sense NaCl delivery.

The GFR affects solute delivery—a high delivery indicates that the filtration rate is too fast for solutes to be reabsorbed sufficiently, while a low delivery indicates that the filtration is so slow that too much of the solutes are being reabsorbed. The macula densa senses flow rate through sodium delivery and sends a signal, the TGF signal, to the afferent arteriole to initiate vasoconstriction to regulate the flow rate into the nephron.

The TGF signal is generated by cells in the distal convoluted tubule (DCT) wall, and controls the diameter of the afferent arteriole (located at (1) in figure 1). The signal is not instantaneous, but has a delay of approximately 3–5 s (Sosnovtseva *et al*. 2003; Pitman *et al*. 2004), and there is a further delay of approximately 3–5 s for the filtrate to flow from the glomerulus to the macula densa (Holstein-Rathlou & Marsh 1990; Sosnovtseva *et al*. 2003). This negative feedback mechanism has been the focus of a number of multi-nephron models, and is known to have an inverse sigmoid relationship with sodium delivery at the macula densa (Schnermann *et al*. 1970).

Nephrons are known to *interact* with each other in multi-nephron systems. Two different types of interaction are known to occur (Holstein-Rathlou *et al*. 2001; Sosnovtseva *et al*. 2003; Marsh *et al*. 2007): *vascular signalling*, where the TGF signal from one nephron also reaches the afferent arteriole of a nearby nephron; and *haemodynamic coupling*, where the resistance of an afferent arteriole affects the hydrostatic pressure at a nearby nephron. Experimentally, vascular signalling favours in-phase synchronization, while haemodynamic coupling favours anti-phase synchronization.

### (b) Computational renal modelling

A large number of mathematical nephron models exist and so our review will necessarily be brief, focusing on the salient models in recent years and the most relevant to the current work.

One focus in recent years has been to investigate the role of TGF in autoregulation and to explain the observed oscillations associated with TGF. Oscillations in the pressure, flow rates of plasma in the proximal tube and chloride concentration in the early distal tube have all been experimentally observed in rats (Leyssac & Baumbach 1983; Leyssac & Holstein-Rathlou 1986; Holstein-Rathlou & Marsh 1989). The models by Holstein-Rathlou & Marsh (1990) and Holstein-Rathlou *et al*. (1991) hypothesized that these oscillations were the result of the TGF mechanism and these authors constructed their models to verify this hypothesis. Another interesting set of models are those that attempt to explore the role of vascular signalling and haemodynamic coupling in the production of oscillations (Holstein-Rathlou *et al*. 2001; Sosnovtseva *et al*. 2003; Marsh *et al*. 2007). Vascular signalling occurs when the TGF signal from one nephron reaches nearby nephrons, and haemodynamic coupling occurs when the resistance of an afferent arteriole affects the hydrostatic pressure at nearby nephrons.

Other models have been more concerned with understanding tubular transport. Weinstein's (1986*a*,*b*) proximal tubule model investigates the distribution of solutes in the lumen of a small length of a rat's proximal tube and the reabsorption of solutes and water. Weinstein's (2005*a*,*b*) model of the DCT models transport of solutes and potassium secretion. Chang & Fujita's (2001) model also models water and solute reabsorption, and potassium secretion. These models are able to reproduce a wide range of experimental data, but there still remains considerable uncertainty about transporters along the nephron tubule (Weinstein 2003).

Multi-nephron models have been used to study specific phenomena such as the ability of nephrons to regulate the sodium gradient in the medullary fluid (Stephenson 1992; Thomas 1998, 2000; Hervy & Thomas 2003; Layton & Layton 2005*a*,*b*; Thomas *et al*. 2006) and in-phase synchronizations in filtration rate, which can arise due to pressure oscillations in blood flow (Holstein-Rathlou *et al*. 2001; Sosnovtseva *et al*. 2003; Pitman *et al*. 2004; Marsh *et al*. 2007). These models, however, do not model the entire nephron. Multi-nephron models are also capable of reproducing a range of experimentally observed phenomena. However, interactions that may arise via other channels, such as changes in the interstitial fluid due to solute secretion or reabsorption by another segment of the nephron tubule, have not yet been explored.

Two observations can be made about these models: (i) they are highly specialized models of nephron tubules, and as such focus on specific nephron behaviours and (ii) the interactions between different nephron tubules, and indeed between nephrons and their surrounding environment, are tacitly embedded within the models, usually as coupling between equations.

## 3. Models

In this section, we begin by briefly defining what we mean by a network automaton. With this done, we can then introduce our nephron model, its underlying graph structure, state space and update rules. Once we have an understanding of the model of a nephron, we can use it to create systems of nephrons using the same principles of model construction as for the single nephron.

### (a) Network automata

Intuitively, a network automaton is a graph in which each node in the graph contains a number of state variables whose values determine the current state for the node and an update rule that specifies how the values of the state variables change over time.

More formally, let *G*=(*V*, *E*) be a graph with vertices *V* and edges *E*. A network automaton on a graph *G* is given by a quintuple *A*={*G*, *S*, *C*, *N*, {*f*_{v}}_{v∈V}} and consists of a finite or countable state space *S*, a configuration of the vertices *C*∈*S*^{V}, a neighbourhood function *N*:*V*→2^{S} and an update rule *f*_{v}:*S*_{v}→*S* for each *v*∈*V*, and where *S*_{v}=*N*(*v*) is the domain of *f*_{v}. The neighbourhood of a node *v* is the set of all vertices that are connected to *v* by an edge.

The network evolves in discrete time. A configuration of the network *C*_{t} at time *t* is given by the value of the state variables of each node at time *t*. The configuration of the automaton at time *t*+Δ*t* is obtained from the configuration at time *t* by applying every update rule for every node in *C*_{t}. This is the *computation rule* for the automaton.

### (b) The single nephron model

The single nephron model is a network automaton. The model structure for the single nephron model is shown in figure 2. The nodes of the network model tubules and the edges model interactions between tubules. Nephrons are immersed in interstitial fluid and peritubular capillaries with which they exchange solutes. The single nephron model must therefore capture the exchange of solutes and water with the interstitial fluid and the peritubular capillaries. To facilitate the modelling of transport of solutes between the nephron and the surrounding interstitial fluid, we add an additional node to represent each layer of the interstitial fluid (not shown in figure 2).

Edges in the model represent significant interactions between tubule segments, which are (i) the flow of a volume of fluid from one tubule segment to the next, (ii) the signal between the macula densa and the glomerulus, and (iii) the transport of solutes and water between nephron segments and the interstitial fluid nodes (also not shown in figure 2).

The dynamics of the system is governed by the update rules. The update rules for the tubules are equations that determine the flow of blood from one tubule to another and the transport of water and solutes between the tubule and its environment.

#### (i) The state and update rule for the single nephron model

The *state* of every node contains variables that represent the quantities of solutes in the blood flowing into that node. Thus, the state space of a node is at least *R*^{N} if *N* solutes are being modelled. In our model, water and sodium are always present as solutes. Different nodes may have state variables in addition to those of the blood. In the model below, there are two groups of such additional state variables: those for the glomerulus and those for the ascending thick loop of Henle. The additional state variables for the glomerulus node are specified in table 1. The ascending thick limb of the loop of Henle has only one additional state variable Na_{M}(*t*), which is the history of sodium delivery to the macula densa and is calculated from the composition of the filtrate leaving the thick ascending limb.

#### (ii) Update rule and filtration rate for the glomerulus node

The update rules for the glomerulus are given by the following set of equations:(3.1)(3.2)(3.3)(3.4)and(3.5)where *Γ* is the viscous resistance coefficient; *H*_{A} is the average afferent haematocrit; *P*_{A} is the arterial pressure; *P*_{V} is the venous pressure; *ζ* is the damping coefficient in the TGF signal; *ω* is the natural frequency of the TGF mechanism; Na_{1/2} is the sodium delivery for which the TGF response is half maximal; *D*_{TGF} is the delay in sending the TGF signal; and *κ* is the range constant for the TGF signal. The glomerulus state variables, and their associated update rules, are discrete forms of the equations for the same variables in Holstein-Rathlou & Marsh (1990), while the formula for calculating oncotic pressure *Π*(*C*) (and the associated constants *a* and *b*) come from Deen *et al*. (1972). The variables *C*_{A}(*t*) and *C*_{E}(*t*) are calculated from the composition of the plasma at the afferent and efferent arterioles, respectively.

Finally, the single nephron glomerular filtration rate (SNGFR) at the glomerulus comes from the standard filtration equation(3.6)where *P*_{BC} is the pressure at the start of the tubule and *K*_{f} is the permeability constant for the glomerulus (Guyton & Hall 2006).

#### (iii) Update rule for the proximal tube node

The proximal tube reabsorption is essentially uniform (Thomas *et al*. 2006), meaning that sodium and water are reabsorbed in approximately equal proportion. Furthermore, the proximal tubule normally reabsorbs approximately two-thirds of the filtered load of sodium and water (Guyton & Hall 2006). The update rules only require that we capture reabsorption of water and sodium,(3.7)and(3.8)This constant behaviour is clearly far simpler than the existing proximal tubule models, but it is able to capture some of the essential dynamics of the tubule segment.

#### (iv) Update rules for the loop of Henle

The loop of Henle is more difficult than the other nodes because it is here that some critical interactions between nephrons take place. The descending limb only permits the diffusion of water and the ascending limb only the active transport of sodium. A large number of other transport mechanisms are known to occur in the limb (Thomas 1998; Hervy & Thomas 2003), but are currently not modelled.

We cannot simply model the ascending and descending limbs using a single node for each. This is because such a topology would not permit us to model the countercurrent exchange mechanism, and in particular the salt gradient along the length of the loop of Henle. Instead, we decompose the nodes for the ascending limb, descending limb and interstitial fluid into a number of layers. Consider a horizontal slice of the loop of Henle at time *t*, as shown in figure 3. At each time step, the model simulates the diffusion of water from the descending limb and the active transport of sodium from the ascending limb. The model assumes that both diffusion and active transport reach equilibrium during the time step. This means that: (i) enough water is diffused to balance the ratio of sodium to water in the interstitial fluid and in the descending limb and (ii) enough sodium is actively transported to maintain the gradient Na_{G} between the interstitial fluid and the ascending limb.

In the equations below, the subscript *M* refers to a quantity in the interstitial fluid, the subscript *A* to a quantity in the ascending limb, the subscript *D* to the descending limb, the subscript *D+M* to a combination of the quantity in the descending limb and interstitial fluid, and the subscript *M+A* to a combination of the quantity in the interstitial fluid and ascending limb.

Firstly, the model needs to conserve the total water and sodium in the system, so we have(3.9)The descending limb reaches an equilibrium with the medullary fluid, maintaining a sodium concentration that is identical to the sodium concentration in the medullary fluid,(3.10)(3.11)and(3.12)The ascending limb reaches an equilibrium with the medullary fluid, maintaining a sodium concentration that is Na_{G} less than the sodium concentration in the medullary fluid,(3.13)(3.14)and(3.15)

There is, however, a problem with this model when it comes to simulations with two (or more) nephron models. Simulations result in large oscillations in the medullary sodium concentration as both the loops of Henle compensate for any imbalance in the equilibrium. The assumption in the equations is that each loop of Henle reaches equilibrium in isolation. The solution is to ensure that each layer of interstitial fluid reaches equilibrium with all of the descending and ascending limbs of Henle interacting.

For a system of *n* nephrons in a *single column* of interstitial fluid, we assume that the *n* nephrons cooperatively maintain the sodium gradient. Each nephron reabsorbs some volume *w*_{i} of water from the descending limb D_{i}, and some amount *s*_{i} of salt from the ascending limb A_{i} for each *i*∈{1, …, *n*}. At equilibrium, we obtain three equations for the water transport in the descending limbs and three equations for the sodium transport in the ascending limb to give(3.16)(3.17)and(3.18)For water and sodium, we have(3.19)(3.20)and(3.21)These equations need to be solved simultaneously for the time step. To do this, the equations for water transport and those for sodium transport can be rewritten in matrix form as *Ax=b* to solve for *x*, where *x*^{T}=[*w*_{1}, …, *w*_{n}] or *x*^{T}=[*s*_{1}, …, *s*_{n}] depending on whether we are solving for water transport or sodium transport. Currently, this is done with Gaussian elimination.

#### (v) Update rules for the distal tube and the collecting ducts

The *early distal tubule* is similar in function to the thick ascending limb and is dealt with in the upper loop of Henle. The distal tubule node models the dynamics of the *late distal tubule*. Here, sodium is reabsorbed at a rate that is regulated by hormones. In particular, the reabsorption rate is proportional to the aldosterone concentration in the filtrate. The permeability of the tubule to water is controlled by the concentration of antidiuretic hormone (ADH). When the ADH concentration is high, the tubule is permeable to water, but in the absence of ADH the tubule is essentially impermeable to water. Consequently, the update rules for the later distal tubule depend on the concentration of ADH and aldosterone in the filtrate,(3.22)(3.23)(3.24)and(3.25)where Na(*t*)_{R} is the volume of filtrate sodium Na(*t*)_{I} multiplied by the tubule response *R*_{ALD}(*t*) and a parameter Na_{max}. This parameter represents the maximal amount of sodium reabsorption that can be achieved by the tubule, and is currently fixed at 0.9.

The cortical collecting ducts have functional characteristics similar to those of the late distal tubule. In particular, they exhibit the same reabsorption characteristics for sodium and water. The medullary collecting ducts, while differing in function from the cortical ducts in several ways, also exhibit the same reabsorption characteristics for sodium and water. Accordingly, the update rules for the collecting ducts are the same as those for the distal tubule.

### (c) Multi-nephron systems

The extension of the single nephron model to a multi-nephron model is achieved by connecting up nephrons into a network. This is done by employing the same principles as in the construction of the single nephron model, where the nodes model the components of the system and the edges model the interactions between components, in this case including interactions between nephrons. For the purposes of this paper, nephrons are assumed to interact with nephrons in their own renal pyramid and not with nephrons in other renal pyramids.

Firstly, a two-nephron model is obtained by connecting two nephrons to the same interstitial fluid nodes. We refer to the interstitial nodes collectively as a column of fluid. This connecting up process can be generalized to create a single-column multi-nephron model. For small numbers of nephrons (typically up to approx. 72), it is sufficient to model the nephrons as parallel resistors (resisting the flow of blood) placed across a potential equal to the difference between the arterial and vascular pressures as shown in figure 4*b* (Marsh *et al*. 2007). In our model, this resistance is the resistance of the afferent arteriole. Thus, by changing the resistance of its afferent arteriole, a nephron can influence the GFRs of the nephrons in the same column of fluid. Changes in resistance model the effects of haemodynamic coupling.

Multi-column, multi-nephron networks are also possible. Again, this is achieved by connecting up components, in this case single-column multi-nephron models, into a network. The main interaction between columns is through blood flow, and edges model the arterial trees that supply blood to the single-column multi-nephron models as shown in figure 4*c*. Doing this also separates nephrons into local neighbourhoods, and allows the different columns of interstitial fluid to evolve inhomogeneously.

To simulate a single time step for this model, we simulate the time step for each column of renal fluid, and then permit diffusion between neighbouring columns. We make two assumptions concerning the nature of this diffusion. Firstly, we only permit horizontal diffusion, i.e. there is no diffusion between different layers of the interstitial fluid. Secondly, we only permit the diffusion to act on one-quarter of the fluid in each column, as this prevents any portion of the interstitial fluid from undergoing intercolumn diffusion with more than one column—each quarter of the fluid diffuses with the neighbour on the nearest side (if any).

## 4. Model dynamics

The single nephron system reproduced experimentally observed results, such as the production of regular oscillations in filtration rate in response to the TGF mechanism (figure 5*a*), and the construction of a stable solute gradient in the medullary fluid (figure 5*b*). A two-nephron system also reproduced the experimentally observed behaviour: (i) synchronizations arose between the filtration rates, (ii) the effects of haemodynamic coupling and vascular signalling were consistent with those observed in existing two-nephron models, and (iii) the system maintained a larger, more stable solute gradient in the medullary fluid than that maintained by the single nephron system (figure 5*c*). This is particularly interesting in view of the role of nephrons in autoregulation.

Larger systems of eight and 72 nephrons have also been explored in which the dynamics was affected by haemodynamic coupling and vascular signalling. Power spectral densities were used to determine the dominant frequencies in the whole-system filtration rate (WSFR), and demonstrated that the degree of vascular signalling not only affected the average WSFR but also altered the dominant frequencies. In the 72-nephron system, the presence of vascular signalling resulted in peaks at higher frequencies than when no vascular signalling was present.

More interesting is that the multi-nephron model has been used to investigate the effects of perturbations and impaired nephron function. Consider a hypothetical impairment of kidney function based on *proteinuria*, which is the presence of abnormally large amounts of protein in the urine and a sign of potential renal damage. One effect of this symptom is an increased oncotic pressure in the nephron tubule and this impedes the diffusion of water into the interstitial fluid. This effect is modelled as an impairment that reduces the reabsorption of water in the descending limb of Henle by some fixed amount *σ*. We outline the results of four experiments: two on an eight-nephron system and two on a 72-nephron system, where the average system filtration rate and the range in system filtration rate were used as measures of stability.

Consider a system with eight nephrons that are *impaired* with reduced water reabsorption levels and containing four pairs of nephrons, with each pair having longer loops of Henle than the previous pair. This system was most strongly affected when the nephrons with the longest loops of Henle were impaired. However, even when the individual nephron function was impaired to 50 per cent of normal levels (*σ*=0.5), the effect on the whole-system dynamics was measured at approximately 0.7 per cent of the unimpaired system filtration rate.

More experiments were conducted on a 72-nephron system that consisted of nine columns—each containing eight nephrons having the same structure as the previous eight-nephron system. The first experiment investigated the effects of two diseased nephrons in each column simultaneously, and resulted in a maximum effect of 1.3 per cent on the system. In the second experiment, all of the nephrons in a single column were diseased simultaneously. This experiment demonstrated the strong averaging effect of the inter-column diffusion mechanism, as the effect on the system was only 0.3 per cent.

In conclusion, this preliminary investigation into the effects of impaired function on the whole-system dynamics demonstrated that the models are highly resilient to such impairments. Despite the severe simplicity of the impairment model that was used, we contend that this resilience is a feature of the model, and of the kidney proper.

## 5. Future research directions and opportunities

The models and simulation experiments discussed in §§3 and 4 clearly provide supporting evidence that a hierarchical dynamical network-based approach renders a whole-kidney simulation of 10^{6} nephrons possible. We now outline a number of exciting research opportunities that can be addressed using our computational multi-nephron model as a starting point.

The first problem is to complete the model to more accurately predict whole-kidney function. To do this, the model must be extended to capture the reabsorption and secretion of many more solutes, including potassium, chloride, calcium, urea and glucose. In addition, the synthesis of hormones by the nephron must also be incorporated into the model. Update rules to capture these dynamics can be based on the existing, highly detailed models of these processes.

The next problem for any large network model, especially one of 10^{6} nephrons, is to create the model in the first place. Models of 72 nephrons can be created by hand, but the construction of larger models takes significant time. One method for automating the construction of large models is suggested in Moss (2009), where nephron systems are *grown*. The method uses a parametrized *graph rewriting system*, a Lindenmayer (1968) system, to generate the model structure for the arterial tree, and uses *stochastic constraints* to select and connect individual nephron models to this tree. Studies of the fractal dimensions of the branching properties of renal arterial trees (Zamir 2001) have been shown to exhibit variations in branching parameters, resulting in multiple fractal dimensions for a single arterial tree. Furthermore, structural morphology such as radial, length and connectivity data for the arterial and venous trees of the rat kidney is available (Nordsletten *et al*. 2006). The fractal dimensions and morphological data can be used as parameters to the Lindenmayer system to grow *realistic* kidney models.

Network models can also be used to study multi-organ systems, such as two interacting kidneys, or the endocrine system. The modelling methodology presented in §3 can be easily extended to such systems, by adding additional layers of hierarchy to the model, and capturing the structure and interactions between the multiple organs at the highest level of the hierarchy. Hierarchical dynamical networks lend themselves to running on distributed systems. The multi-nephron models can easily be made to run concurrently by distributing blocks of renal fluid and nephrons to other processors.

Section 4 reports some of the effects of impairments on eight- and 72-nephron models and models just some of the effects of a renal disease such as proteinuria. Of course, identical experiments can be conducted on whole-kidney models when these become available. Exploring the onset and effects of renal diseases opens up a particularly interesting venue of research for our network models.

The causes and mechanisms of renal disease have been investigated at the cellular level and have been the subject of many clinical experiments and case studies, but there appear to be few, if any, attempts to model the onset, spread or impact of renal disease in multi-nephron systems. This lack of disease in multi-nephron systems precludes the study of emergent kidney dynamics in the presence of disease. Our hierarchical dynamical networks model offers perhaps the first such opportunity.

## 6. Conclusion

We have presented a novel approach to modelling multi-nephron systems, which combines graph automata and complex network approaches into a single model, allowing the connectivity and interactions in the model to be altered interactively. This model has reproduced the experimentally observed behaviour, and has been shown to scale to whole-kidney models. This ability to model large multi-nephron systems permits exciting new research possibilities, based on studying how the system dynamics emerge from and are affected by low-level couplings and interactions. A preliminary investigation into the effects of localized impairment on the whole-system dynamics and stability was given as an example of such a study.

## Acknowledgments

The authors would like to thank Donald Marsh and Randy Thomas for their feedback and commentary.

## Footnotes

One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications II’.

- © 2009 The Royal Society