- Email: [email protected]
A complete steady state model of solute and water transport in the kidney
COMPUTER MODELLING PERGAMON
Mathematical and Computer Modelling 29 (1999) 63-82
A Complete Steady State Model of Solute and Water Transport in the Kidney N. E. KOTTLER Department of Mathematics, Box 8205 North Carolina State University Raleigh, NC 27695-8205, U.S.A. Ii. T. TRAN* Center for Research in Scientific Computation, North Carolina State University Raleigh, NC 27695-8205 U.S.A. trana)control.math.ncsu.edu
Box 8205
D. E. WESSELL Department of Mathematics, Box 8205 North Carolina State University Raleigh, NC 27695-8205, U.S.A. (Received June 1998; accepted July 1998) Abstract-The purpose of this paper is to incorporate a detailed model, along with an optimized set of parameters for the proximal tubule, into J. L. Stephenson’s current central core model of the nephron. In this model a set of equations for the proximal tubule are combined with Stephenson’s equations for the remaining four tubules and interstitium, to form a complete nonlinear system of 34 ordinary differential and algebraic equations governing fluid and solute flow in the kidney. These equations are then discretised by the Crank-Nicholson scheme to form an algebraic system of nonlinear equations for the unknown concentrations, flows, hydrostatic pressure, and potentials. The resulting system is solved via factored secant update with a finite-difference approximation to the Jacobian. Finally, numerical simulations performed on the model showed that the modeled behavior approximates, in a general way, the physiological mechanisms of solvent and solute flow in the kidney. @ 1999 Elsevier Science Ltd. All rights reserved. KeywordsTransport model, Biomathemat~cs,Differential algebraic system, Inserve problem, Numerical methods.
1. INTRODUCTION In the late 1960s and early 1970s many researchers began to develop mathematical models of the kidney [l-4]. This research has progressed from the mid-1970s to this day with many variations, additions, and improvements [5-lo]. Some of these models have simply tried to fit the empirical data without regard to modeling the actual physiological processes. Researchers who were trying to develop physiologically accurate models encountered a substantial problem in modeling the mechanism that returns recycled water and salt from the kidney to the generaL circulation; this function is performed by a network of blood vessels in the kidney called the vasa recta. The central core model of the urine concentrating mechanism, which was initially proposed by Stephenson [4] *This research was supported in part by the National Aeronautics and Space Administration Grant NAGW-2865. 0895-7177/99/$ - see front matter @ 1999 Elsevier Science Ltd. All rights reserved. PII:SO895-7177(99)00055-2
Typeset by AM-T@
N. E. KOTTLERet
64
d.
in 1973, greatly simplified the modeling of the vasa recta while retaining all of the important physiological aspects of kidney function. This model quickly became the standard. In the past In 1987, Stephenson
20 years there have been many enhancements to the original model [ll-171.
et al. [16] considered for the first time electrolyte transport and introduced electric potential as a driving force for passive solute movement.
Other enhancements have been primarily due to
a greater understanding of renal physiology and to advances in computer technology that have allowed us to solve the model equations numerically [18]. As advanced as the model has become, however, there has been a disregard of the proximal tubule. In this paper, we extend Stephenson’s central core model of the nephron to include the proximal tubule. The physiologically based model contains a set of equations and parameters for all five nephron tubules, the medullary interstitium, and constants for the cortical interstitium. Because the model is an extension of that presented by Stephenson et al., parameter data is taken directly from Stephenson’s papers [16,17]. In addition, since the parameters for the proximal tubule were unavailable, they are estimated based on a nonlinear least-squares method with the cost function derived from the known behavior of the proximal tubule. The organization of the paper is as follows. We begin in Section 2 with a brief description of renal physiology. This summary of the renal physiology is important for the understanding of the mathematical model developed in subsequent sections. In Section 3, we first develop a set of differential and algebraic equations governing the flow of solvent and solutes in an arbitrary section of the nephron, then we customize these equations to each nephron tubule and to the interstitium. Section 4 presents our efforts on improving the model of the proximal tubule. The parameters used in the proximal tubule model are estimated based on a nonlinear least-squares method. The cost function used in the least-squares method is derived from the physiologically known behavior of the proximal tubule. Section 5 summarizes the steady-state nephron model of the kidney and discusses the numerical methods used to solve the nephron model and presents its numerical solution. Finally, Section 6 contains our concluding remarks including a comp~~on of our solution with the known behavior of the kidney.
2. THE KIDNEY The functional unit of the kidney, the nephron, is composed of a glomerulus and five connected tubules, the proximal tubule (PT), the descending limb of the loop of Henle (DHL), the ascending limb of the loop of Henle (AHL), the distal tubule (DT), and the collecting duct (CD). This general structure of the nephron is depicted in Figure 1. The tubules form a continuous system which modifies the filtrate and transports it to a cavity in the kidney called the renal pelvis. This filtrate which passes through the glomerulus is referred to as glomerular ultrafiltrate (GU). The combined effects of each tubule enable the kidney to form a concentrated waste product while maintaining homeostasis in the body-the primary function of the kidney. Each tubule has its own role integral to the success of this process. 2.1. Interstitium There are two interstitia within which the nephron lies, the cortical interstitium (CI), and the medullary interstitium (MI). While the MI has a set of equations modeling its compIex behavior, the CI is simply assumed to be a well-mixed bath of constant volume. The PT and DT lie in the CI; the DHL, AHL, and CD lie in the MI. 2.2. Proximal
Tubule
The PT acts as a large sieve, returning the majority of what it is delivered from the glomerulus back to the general circulation. This is accomplished by a sodium potassium pump which actively transports sodium and potassium ions into the CI. The active transport favors a passive movement of chloride ions (via the potential gradient) and water (via the osmotic gradient) into the CI.
SoIute and Water T’rsnsport
Figure 1. General structure
Because
the net flux of PIat- is matched
by water,
65
af the nephron.
sodium
absorption
is isotonic
along the length
of the PT, Movement of urea is primarily due to the effects of its concentration gradient and solvent drag. At the end of the proximal tubule, approximately 70% of the Hz0 and Na+, 85% of the K+, 42% of the Cl-, and 50% of the urea that entered the proximal tubule have been returned to the general circulation via the cortical interstitium. The goal of forming a concentrated waste product (i.e., urine) has already begun in the proximal tubule. The function af the remaining tubules of the nephron is to enhance this concentrating process. 2.3.
The Descending Limb of Henle’s Loop
Active transport of Na+ in the AHL causes the formation of a sodium in the MI. As the GU traverses down the DHL, the exterior environment more concentrated. In fact, by the bottom of the DHL, in the deepest concentration levels in the MI reach three to four times what they are The walls of the DHL are relatively impermeable to the passive diffusion very smalf amounts of Na+ and Cl- enter the DHL to equilibrate this
~oneentrati~~ gradient becomes increasingly region of the meduila, in the upper medulla. of NaCl. Hence, only concentration gradient.
The DHL, however, is very permeable to water. Thus water is drawn out of the DHL by osmosis and enters the capillary network to return to the vascular system. There is no active transport in the DHL; and solvent drag, concentration gradients, and electrical gradients lead only to minor movements of Nai, Cl-, and IS+. Simil~ly~ movement of urea is affected by solvent drag and the urea concentration gradient. In essence, the descending limb of Henle’s Ioop accomplishes two main tasks: the extrusion of an additional 50% of the water, and hence, the equilibration of the concentration gradient between the GU and the meduliary interstitium. Thus, by the end of the descending limb the tubular fluid is isosmotic with the exterior enviranment, and a highly concentrated filtrate is delivered to the ascending limb of Henle’s loop. 2.4. The Ascending Limb of Henle’s Loop The system most commonly regarded nephron, counter-current multiplication,
as the backbone of the urine ~on&entrat~ng ability of the resides in the loop of Henle. Saft is actively pumped
66
N. E. KOTTLER
et al.
out of the AHL creating a concentration gradient in the MI. As the filtrate in the DHL passes through the concentrated MI, water flows out of the filtrate by osmosis. Hence, the DHL delivers a solute concentrated filtrate to the AHL, allowing for the extrusion of even more salt, and so on. Protecting the concentration gradient in the MI, the walls of the AHL are impermeable to water. Therefore flow remains constant in the AHL. In addition, because there is no fiux of water out of the AHL, there is no solvent drag. Thus, urea is moved only by the urea concentration
gradient.
In essence, the counter-current mechanism traps large amounts of sodium in the medullary interstitium, enabling large amounts of water to be reabsorbed in the descending limb. It also lowers the osmolality
of the filtrate
as it moves up the ascending
limb toward
the distal convoluted
tubule. 2.5. The Distal Convoluted As a result
of the counter-current
Tubule multiplier
system,
the osmolality
of the GU is slightly
lower
than that of the CI, causing water to leave the DT by osmosis. The two forces moving urea are the urea concentration gradient and solvent drag. The reabsorption of Na+ and K+ is controlled by the hormone aldosterone. Aldosterone is secreted by the adrenal cortex in response to a high K+ concentration in the blood. An increase in aldosterone secretion yields an increase in Naf reabsorption (out of the DT), and an increase in K+ secretion (into the distal tubule). In fact, without any aldosterone secretion, 100% of the K+ is reabsorbed into the vascular system; whereas with aldosterone secretion at its highest, 100% of the Naf is reabsorbed. Thus, aldosterone can play a powerful role in maintaining homeostasis. For the purposes of this paper, an average aldosterone secretion, which would occur under hom~static conditions, is assumed. For a more complete and accurate model, a control system representing aldosterone secretion must be incorporated. 2.6.
The Collecting Duct
As a result of the counter-current multiplier system, the exterior environment becomes increasingly hypertonic as the filtrate traverses down the CD. In addition, the counter-current multiplier system has caused the fluid entering the CD to be hypotonic (despite the efforts of the DT). Therefore, counter-current multiplication creates a large concentration gradient across the wall of the CD. Salt does not enter the CD to relieve the concentration gradient because the walls of the CD are impermeable to salt. Furthermore, there is no active transport of any solute in the CD. Hence, strong osmotic effects cause water to be drawn into the MI, where it enters the capillary network and is returned to the general circulation. The rate at which water is reabsorbed is controlled by ADH, antidiuretic hormone. ADH, produced by the neurons in the hypothalamus and secreted by the posterior pituitary gland, increases the permeability of the collecting duct to water. For example, in a state of severe dehydration ADH is secreted, and up to 99.8% of the original GU is reabsorbed (the remaining 0.2% is the obligatory water loss). In a well hydrated person, ADH would not be secreted, and 99.2% of the original GU is reabsorbed. This control translates to urine volumes ranging from as small as 0.4L (dehydrated person), to 1.5 L (well hydrated person). In this paper, an average ADH concentration is assumed. Therefore, while this model may be accurate under homeostatic conditions, it will not accurately represent any deviations. As an expansion of this model, an ADH control system should be modeled. At this point, the processing by the kidney is complete. The urine flows from the collecting duct via the ureters to the bladder. The entire process is summarized in Tables 1 and 2. 3.
MATHEMATICAL
DESCFUPTION
The solutes we will model are Na’, Cl-, K+, and urea. Figure 2 shows a section of a genera1 tubule of the nephron. The mass balance equation for solute k in the differential volume element
67
Solute and Water Transport Table 1. Relative flow rates of water and relative amounts of sodium, urea, and potassium in different segments of the nephron of a rat kidney.
I
Nephron Segments Gtomerulus
Percentage Remaining Potassium Sodium Urea
Water 190
1ilQ
100
IQ6
End of Proximal Tubule End of Descending Limb
23.0 17.0
23.0 47.1
44.7 605.3
End of Ascending Limb
17.0
9.5
133.0
3.7
2.9 0.1
77.1 26.3
26.7 11.9
End of Distal Tubule End of Collecting Duct
6.0 0.24
16.4 86.9
Table 2. Concentrations of sodium, urea and potassium in different segments of the nephron of a rat kidney.
I
Concentration in mmol,/L
Nephron Segments
Sodium
Urea
Potassium
Gfomerulus End of Proximal Tubule
140.0 140.0
3.0 5.0
6.3 4.5
End of Descending Limb
294.4
31.8
10.3
End of Ascending Limb End of Distal Tubule
101.0 89.7
35.4 87.8
3.4 8.3
End of Collecting Duct
176.2
243.3
71.8
J J
fP f+
i
t
1
)
!
t,
*
!
x,=2
i f !
a
X,=L,
Figure 2. A general flow tubule. of tubule i shown in Figure
2 is given by
where AZ; is the length of the dierentiaf vdume element, Ai is the cross-sectional area of the differentid volume element, C:kfzi, t> is the ~~n~~nt~at~~nof solute k in tubule i, Fikfz+, t) is
N.
68
the axial solute flow of solute tubule
destroyed in tubule The corresponding
i, Jik is the net outward
k in tubule
i, and Sik is the average
E. KOTTLER et al.
net rate per unit volume
transmural
at which solute
flux of solute
k in
k is being created
or
i. mass balance d [AziAi]
for volume
(i.e., water)
flow in tubule
i is given by
= Fiv (xi, t) - Fiv (xi + A~ir t) - Axi Jiv,
at where Fiv is the axial volume
equation
i, and Jiv is the transmural
flow in tubule
t) with its Taylor series expansion,
In both cases, replacing the term F~,(xi+Az:i, and letting Axi ---f 0, yields, respectively,
a[Ad& (Xi,t)] = _ aFik(xi7t, _ Jik + at and
volume
&ik
aFiv (xi,t) ax
Equation (3) gives a partial differential equation arbitrary flow tubule. Equation (4) gives a partial in an arbitrary flow tubule.
dividing
i.
by Axi,
(3)
dX
aAi -=at
flux in tubule
_ J,
ZV.
for mass balance of an arbitrary solute in an differential equation for mass balance of water
For the purposes of this paper, we are interested only in the steady state problem, and thus, the left-hand sides of both equations (3) and (4) become zero. In addition, we make the simplifying assumption that no solute is being produced or destroyed in any tubule (i.e., sik = 0). Hence, the mass balance equations simplify to
(Xi,t> = _ J,
a&k
tk
ax and
8FiV (xilt) = _ ax
The transmural pressure gradients written as
J,
zv.
(6)
volume flux, Jiv, is determined by a combination that exist between the tubule and the interstitium.
Jiv = hiv
2
RT (elk - Cik) gik f Pi - PI
k=l
of the concentration and In general Jiv can be
1,
(7)
where hiv is the hydraulic permeability coefficient per unit length for tubule i, R is the universal gas constant, T is the absolute temperature, CIk is the concentration of solute k in the interstitium surrounding tubule i, C& is the concentration of solute k in tubule i, oik is the Staverman reflection coefficient of the wall of the ith tubule to solute k, Pi is the pressure in tubule i, and PI is the pressure in the interstitium surrounding tubule i. The transmural solute flux, Jik, is determined by a combination of solvent drag, Fick’s Law, electronic attraction, and active transport. Hence, Jik can be written as
Jik = J~v (1 -
aik>
CIk
i- Cik n L
+
hik (CIk - Cik) +
‘Ik; ‘jkg
(‘Pi -
e,,] +
(8)
J$,
where hik is the passive permeability of tubule i to solute k, zk is the valence of solute k, F is Faraday’s constant, 8i is the potential in tubule i, *I is the potential in the interstitium surrounding tubule i, and J,“kis the active transport of solute k from tubule i.
69
Solute and Water Transport
The active transport, J$, is a metabolically Menten kinetics. that is.
driven transport that is consistent with Michaelis-
where o;k is the m~mum rate of active transport, and &k is the ~on~entratian at which 2% reaches onehalf its m~mum value (Le., i%&is the Michaelis constant). The Nernst-Plank equation relates solute flow, Fit,, to v&me flow, &v, concentration, C&-l and potential !3Zrias follows:
where L&k is the diffusion coefficient of solute k in tubule a, and T&kis the mobility of solute k in tubule i. Furthermore, diffusivity and mobility are related by
Assuming Poiseuille flow, pressure drop along tubule i is given by the equation
where Pi is the hydrostatic pressure in tubule i; Rg = 8qfrf is the hydrodyn~~~ flow resiistance of tubule i; where 77is the viscosity of the solvent (water) and T is the radius of the cylindrical tubule. Since Na+, K+, and Cl- are ions, electroneutrality must be maintained. This condition yields the final equations, 3 c k=l
and
=
0
(13)
3 c k=l
Using
zkC$k
equation
(f3), the concentration
Zks&k =
0.
(14)
of Cf-, a negatively charged ion, can be written in
terms of the ~on~ent~at~on~ of Na+ and K+, the positively charged ions. Hence,
where Cis is the concentration of chloride in tubule i; C+r is the concentration of sodium in tubule i, and Ci3 is the concentration of potassium in tubule i. Finally, substituting equation (8) into equation (I$) and solving for 9i - 91 yieids
i&.-g,
=
5 Zk f&V (1 - %k) ((elk + cik) 12) f hik (elk 5 +ik ((Clk +Cik) 12)(F/W
Cik) + J$]
_k=l
k=l
Our general system can be summarized as follows: in each tubule, i, we have eleven dependent variables: volume flow, Fgv; sodium flow, 2%; sodium concentrations Gil; chloride flow, Fi2; chloride concentration, i&; potassium Bow, S’i3; potassium concentration, CS; urea flow, F&; urea concentrations C&; pressure, Pi; and potential, 96. These variables are matched by nine differential equations: (5), (6), (lo), and (12); and twa algebraic equations: (15) and (IS).
70
N. E. KOTTLER
&al.
In all tubules except the medullary interstitium, we assume that the diffusivity and mobility terms are negligible as compared to the convective flow term. Thus, in these tubules, equation (10) simplifies to the algebraic relationship,
Fik = FivCik. Taking the derivative on both sides and solving for 2
(17) yields
1
-=-
(18)
.
Substituting equations (6) and (5) into (18) gives an equation for the change in concentration of each solute, k. Namely,
JiVCi/c- Jik Fiv ’
acik -= ax
(19)
Once the four solute concentrations are found, equation (17) gives values for the four solute flows, Fi1, Fis, Fi3, Fi4. Since the solute flows can be computed algebraically, the four equations represented by differential equation (5) do not need to be evaluated. Hence, the total number of differential equations in each tubule, except for the medullary interstitium, is reduced from nine to five. These are: -
the volume flow equation a&v -
ax
-
ax
(21)
JivCi3- Ji3 Fiv ’
(22)
the urea concentration equation aCi4 -=
ax
-
JivCil - Jil Fiv ’
the potassium concentration equation aCi3 -= ax
-
(20)
the sodium concentration equation acil -=
-
= -Jiv,
JivCi4
-
Fiv
Ji4
(23)
’
the pressure equation
(24) The remaining values to be computed are chloride concentration, C~Z and potential, +. The chloride concentration, CQ, can be found from equation (15); and Qi can be computed using the algebraic relationship expressed in equation (16). Now we can adapt our general system of equations to each tubule by determining the appropriate forms for Jiv and Jik. 3.1. The Proximal Tubule Model Using equation (7), we can write the equation for volume flux in the proximal tubule,
Jov = hov f: RT (Cck - COk) ffOk + PO- PC , k=l
I
(25)
71
Solute and Water Transport
where CCk is the concentration of solute k in the cortical interstitium, and PC is the pressure in the cortical interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium. cc1
-t Co1
Joi = Jov (1 - COk)
(Cc1
2
-co1>+
“I
+ 2 ‘01
g
(%I
-
Ic)]
+
J&;
(26)
cc3
+ ‘03 2
g
(go
-
qc)]
+
Jg3;
(27)
for potassium,
Jo3
=
Jov
(1 -
cc3
ao3)
+ co3 2
and for urea, Jo4
=
Jov
(1. -
(Cc4
+ ho4
cc4
co4)
-
+ co4 2
Co4)
f
G4
+ G4 2
g
(ipo
-
Bc)]
_
[
Note that the flux equation for urea does not contain an active transport actively transported in the proximal tubule.
term since urea is not
3.2. The Descending Limb of Henle’s Loop Model Using equation (7), we can write the equation for volume flux in the descending limb,
Jlv = hv
&WIX
-
Clk)fllk
+ pl
1 ,
-PI
k=l
(2%
where Cfk is the concentration of solute k in the medullary interstitium; and PI is the pressure in the medullary interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium, Jll
=
&V
(1 -
Cl1
elk)
+ Cl1 ,, L
(30)
for potassium, A3
=
JlV
(1 -
CI3
613)
+ Cl3 2
(31) (C13
+ h13
-
G3)
+
cI3
+2 c13
g
(i&l
_
~~)I
;
+2 (24
g
(\El
_
*,,I
.
[
and for urea, J14 = JIV (1 -
+ hi4
CI4
~14)
((714
-+- Cl4
2 -
64)
(32) +
a4
Note that none of the flux equations contains an active transport transport in the descending limb.
term since there is no active
N. E. KOTTLER et al.
72 3.3.
The Ascending
Limb
of Henle’s
Loop
Model
Since the ascending limb is not permeable to water, there is no volume flux into or out of the ascending limb. Hence, the equation for volume flux is simply, (33)
Jzv = 0. The equations for solute flux in the ascending limb are, for sodium,
+.
J21 = h21 (Cr1 - C21) +
J& ;
(34)
for potassium, J23 = h23 (C13 - C23) + ‘I3 + c23 g 2
(‘P2 - QI)] + J2a33;
and for urea, - C24) + c14 +2 c24 g
J24 = h24 (64
(q2 - I,)]
.
(36)
Note that the flux equations have no solvent drag term, because J2v = 0. Also, the flux equation for urea does not contain an active transport term since urea is not actively transported in the ascending limb. 3.4. The Distal
Tubule
Model
The equation for volume flux in the distal tubule can be expressed as
J3v = h3v
f:
RT (Ccrc -
C3k)
c3k
+
p3 - PC
k=l
1,
(37)
where 63k = 1, V k. Hence, water is a strong force equilibrating the concentration gradient which exists between the distal tubule and the exterior environment, the cortical interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium, J31 = h31 (Cc1 - C31) + cc1 + c31 g 2
(‘P3 - SC)]
+ J&l;
for potassium, J33 = h33 (Cc3 - C33) + co
+ c33 g 2
(\k3- Qc)]
;
(39)
and for urea, J34 = 0.
(40)
Notice that none of the flux equations contains a solvent drag term. The reason for their absence is that the reflection coefficients for all of the solutes in the distal tubule are at their maximum value (i.e., U3k = 1, Vk). Also, there is no active transport of potassium or urea in the distal tubule; thus the absence of the Michaelis-Menton terms in those equations. Finally, the distal tubule is impermeable to urea. Combining this fact with the lack of urea transport and with the reflection of urea at its maximum, there is no net flux of urea into or out of the distal tubule.
73
Solute and Water Transport
3.5.
The Collecting
Duct
Model
The equation for volume flux in the collecting duct can be expressed as
J4v
= h4v
2
RT
(CI~
-
C4k) g4k + p4 -
k=l
PI
1,
(41)
where U& = 1, V k. Similarly, we can write the equations for solute flux. These are, for sodium,
J41 = h41
(Crl
J43 = h43
(C13 -
C41) + cIl
-
+2 c41 g
(Q4 -
*I)]
;
+2 c43 g
(Q4 -
QI)]
;
for potassium, + c13
C43)
(43)
and for urea, J44 = 0.
(44)
Similar to the distal tubule, the reflection coefficients in the collecting duct are at their maximum. Hence, solvent drag is not a factor as all of the solutes are reflected back into the tubule at the walls of the collecting duct. Also, there is no active transport of any solute in the collecting duct, including sodium. Finally, the collecting duct is relatively impermeable to urea. Thus, ail of the urea delivered to the collecting duct remains inside and is transported as waste out of the kidney. 3.6. The Medullary
Interstitium
Model
Diffusivity and mobility terms can no longer be considered negligible as compared to convective flow in the medullary interstitium. Therefore, we no longer have an algebraic relationship between solute flow and solute concentration. However, we know that in the medullary interstitium, the sum of the solute flows at any level is equal to the final excretion of that solute. In other words, we can use the following equation to compute solute flows in the medullary interstitium:
c
Fik(z)
= F4k(N),
i
where i is summed over all tubules which lie in the medullary interstitium (i.e., i = I, 1,2,4); Fik is the flow of solute k in tubule i; N is the final discretized position in any tubule; and Fdk(N) is the final flow of solute k in the collecting duct. We can use equation (17) to compute the solute flows in each tubule except for the interstitium. Incorporating these into equation (45) yields solute flow equations for sodium, FII(~)
+ F~v(~)Cll(z)
+ F2v(z)C21(z)
+ F4v(z)C41(2)
-
F4v(N)C41(N)
= o;
(46)
+ Flv(z)Cl3(z)
+ FZV(z)C23(z)
+ F4v(x)C43(5)
-
F4v(N)C43(N)
= 0;
(47)
+ F1v(x)C14(z)
+ FZV(z)C24(5)
+ F4v(z)C44(2)
-
F4v(N)C44(N)
= 0.
(48)
for potassium, FIN
and for urea, F14(2)
We have a similar condition for solvent flow in the medullary interstitium. Namely, the sum of the volume flows at any level of the medullary interstitium is equal to the final flow in the
N. E. KOTTLER
74
et al.
collecting duct. This condition translates to the following algebraic equation for solvent flow in the meduliary interstitium: E&V(X)
=
F4V(N).
(49)
i
Equivalently, FlV(Z) + Applying the electroneutrality
FlV(X)
+ J72v(s)
+ F4v@)
-
iJ4v(N)
=
0.
(56)
condition (14) to the medullary interstitium, we have our final
equation, QJrk = 6,
e
(511
k=l
where Jlk: = - E
Jirc
(52)
and Jik is the net transmural flux out of each tubuk which lies in the medullary interstitium (i.e., i = 1,2,4). Hence, ZI (-JII
- J21 - J41) +
22
c-J12
- J2z - J42) + 23 t-513
- 523 - J43f = 0.
(53)
However, the flux of chloride out of tubule i was never explicitly calculated because chloride concentration and hence, flow, could be algebraically computed in all other tubules by the electroneutrality condition, (13). Therefore, we need an equation governing the net flux of chloride out of tubule i. This equation comes from equation (5), where Fi2 = FivCi2 in all tubules except for the interstitium. Thus we have WvG2 3:
=3:_ J, 22.
(54)
Differentiating and using the product rule, we obtain
Equation (55) holds in all tubules except for the interstitium. Hence, we can compute the flux of chloride out of the descending limb, Ji2, the ascending limb, J22, and the collecting duct, 542. Using these net chloride fluxes, we can complete equation (53). The remaining four equations for pressure and solute concentrations are as previously stated. However, it should be once again noted that the entire Nernst-Plank equation must be used to compute solute concentrations in the medullary interstitium. 3.7. The Cortical
Interstitium
Constants
Keeping in line with Stephenson’s model, we will assume that the cortical interstitium is a well-mixed bath of constant volume. Thus, concentration of each solute, (Ccl, CCZ, Cca, and Cc*); pressure, (PC); and potential, (tkc), are all constant in the cortical interstitium. Actual values for these constants can be found in Table 7. 4.
PARAMETERS
AND
BOUNDARY
CONDITIONS
In order to complete our model, initial conditions for the five independent variables (volume flow; sodium, potassium, and urea concentrations; and pressure) as well as parameter values are required. When available, these initial conditions and parameters were taken from Stephenson’s papers [16,17]. However, because Stephenson did not model the PT, several modifications and parameter identifications were necessary. Specifically, Stephenson’s initial conditions are those for the DHL. In order to get our initial conditions to match Stephenson’s at the PT-DHL boundary,
Solute and Water Transport
75
we extrapolated his values based on the physiology of the PT (see Tables 1 and 2). For example, Stephenson’s input volume flow into the descending limb is lO.OnL/min.
Table 1 indicates that
this value is 23% of the initial flow in the proximal tubule (in the literature, values range from 20% to 35%). Therefore, a flow of approximately 33 nL/min is used as the input volume flow for this problem. We also note that the initial pressure is a parameter (its calculation will be discussed subsequently).
Table 6 presents the initial conditions for each differential equation. The system
parameters and their values are listed in Table 8. In order to remedy problems associated with varying tubule length and poor scaling (volume flows are on the order of 10-i’ (L/s) while pressure values are on the order of 10 (mmHg)), the system parameters are normalized with respect to tubule length and entering volume flow as follows (note: Z’ denotes the normalized parameter x). For solute permeabilities,
hi/c&
hk= piv(o)
(56)
’
where hik is the nonnormalized permeability of tubule i to solute k; Li is the length of tubule i; and Fiv(O) M 5.5 x lop7 mL/sec is the entering flow rate in the proximal tubule. For hydraulic permeability, hi,
=
$-$
(57)
EV
where hiv is the nonnormalized hydraulic permeability of tubule i. For maximum active transport rates, I aik
=
aikLi ~~~(0)
(58)
’ .
where aik is the nonnormalized maximum active transport rate of solute k in tubule i. For flow resistance, Rip = 7.5062 X 10m4qLiFiv(0),
(59)
where RiF is the nonnormalized flow resistance of tubule i; and v is the viscosity of water in poise (note: 7.5062 x 10m4 is a conversion factor for converting poise to mmHg s). It is necessary to normalize neither the reflection coefficients, gik, which are dimensionless, nor the Michaelis constants, kik, which are independent of flow and tubule length. Nonnormalized parameter values for the DHL, AHL, DT, CD, and MI were taken directly from Stephenson [16,17]. However, no micropuncture data was available for PT parameters. Hence, an ‘optimal’ set of parameters for the PT was determined independently of the other four tubules and MI. This optimization was done so that the PT model behaved as closely as possible to the known behavior of the PT. The parameter estimation problem was formulated as follows: let Q = (o&i,koi,abs, ko3, hb,, h&, hbar hb4,uo4, PO(O)) be the PT parameter set. In order to limit the number of parameters being estimated, we do not estimate the reflection coefficients for Na+, Cl-, or K+; ph ysiological evidence shows these parameters are close to one [14]. In addition, the initial pressure is listed as a parameter. This inclusion is due to the fact that initial pressure varies greatly in the proximal tubule (range 20-40 mmHg); however, the final pressure should be very near 10.75mmHg for a smooth transition to the DHL. Thus, final pressure is taken into account in our cost function and we allow initial pressure to be varied. The parameter estimation problem is to find the set of parameters 4 that minimizes
J (q) = wi fl: [Co1 (Zj) - 0.14012 +
w2
[Fov
(ZN)
-
0.312
j=l + w3 [CO3 (ZN)
-
o.oo45]2 + w4 [Co4 (ZJV) - 0.00512
+ ws [PO(ZN) - 10.7512,
(60)
76
N. E. KOTTLER et al. Table 3. Weights for cost function J. Value
Weight constant sodium concentration
weight, wi
2.0
final volume flow weight, ws final potassium
1.0
cont. weight, w3
1000.0 1000.0
final urea cont. weight, w4 final pressure weight, ws
0.001
Table 4. Parameter bounds for the proximal tubule model. Parameter active transport
Lower Bound
of sodium, a&
._
Upper Bound
0.0
10.0
Michaelis constant for sodium, kor
0.0
10.0
active transport
0.0
10.0
0.0
10.0
of potassium,
a&s
Michaelis constant for potassium, hydraulic permeability,
h&,
sodium & chloride permeability, potassium permeability, urea permeability,
ko3
hb, , hb,
hbs
hb4
urea reflection coefficient,
of34
initial pressure, PO(O),
0.0
1.0
0.0
100.0
0.0
100.0
0.0
100.0
0.0
1.0
20.0
40.0
Table 5. Initial and final parameter values for the proximal tubule model. Parameter active transport
Initial Guess
of sodium, a&
Final Value
0.10
0.41
Michaelis constant for sodium, kol
0.15
0.15
active transport
0.10
0.27
of potassium,
Michaelis constant
for potassium,
hydraulic permeability, permeability,
urea permeability,
h&s
hb4
urea reflection coefficient, initial pressure, PO(O) value of the cost function numbers of iterations
ko3
hb,
sodium & chloride permeability, potassium
a&
004
h& , h&
0.15
0.15
0.001
0.04
1.5
1.4
1.5
2.0
1.5
10.0
0.50 35.0 0.345525
0.20 25.35 0.2968383-
04
22
where N is the final mesh point and Cor(zj), Fcv(srv), Cc3(5~), COG, and Po(z:N) satisfy the system differential equations (20)-(24) d escribing the proximal tubule. The summation term in equation (60) requires that sodium concentration remain constant at O.l40mmol/mL; final volume flow be 30% of its initial value; final potassium and urea concentrations be, respectively, 0.0045 and 0.005 mmol/L; and final pressure be 10.75 mmHg. The weights, 2uj are chosen so that a I% error in any of the five terms adds an equal amount to the cost function, J. The appropriate weights are listed in Table 3. The minimization of the cost function J is done using a finite difference quasi-Newton method with an active set strategy since the problem is one with simple bounds. This scheme is implemented via the IMSL routine DBCONF 1191.For a detailed description of the numerical methods, see [ZO]. The parameters in the vector 4 are physical quantities that are bounded below by zero (except PO) and above by different values as listed in Table 4. The optimized parameter values (in the column “Final Value”) are given in Table 5 along with the initial value and final value of the cost function J. The initial guesses of these parameters
Solute
are taken
to be average
determined
values
by Stephenson
and Water
of the similar
parameters
[16]. We note that
77
Transport
in other
the optimized
segments
parameter
values
of the cost function by an order of magnitude of 4. These final values for the PT parameters are used as approximations in our larger system.
5. THE
NEPHRON
The complete
steady-state
MODEL
AND
34 ordinary differential equations and algebraic equations (DAE) is summarized below.
(0
For each of the five nephron
tubules,
of the kidney
constraints.
there
as
the value
SOLUTION
consists
The system
are five ordinary
reduced
for those PT parameters
ITS NUMERICAL
model of a single nephron
of the nephron
of a coupled
of differential
differential
set of
algebraic
equations,
equa-
of volume (i.e., water) flow, sodium concentrations (20)-(24), d escribing the distributions tion, potassium concentration, urea concentration, and the pressure. In these equations, the expressions for the transmural volume and solute fluxes, J~v and Jik, are given in Sections 3.1-3.5. the concentrations of sodium, potassium, and urea are solved (ii) In the medullary interstitium, from the Nernst-Plank equation (10). In addition, the volume flow and solute (sodium, potassium, and urea) flows satisfy the algebraic constraints (46)-(49). The pressure still satisfies the pressure equation (24) and the transmural solute flux in the medularry interstitium, JIk, satisfies the electroneutr~ity condition (51) or equivalently (53). The total system of 34 differential algebraic equations (29 differenti~ equations and five algebraic equations) is solved numerically by a quasi-Newton method, implemented in the IMSL routine, DNEQBF [19]. For the analysis of the numerical method implemented here and other possible numerical methods for DAE, we refer the reader to Brenan et al. [21] (and the numerous references therein). First, each tubule is divided into N segments. Because our system is normalized with respect to tubule length, each discretized segment has the same length (l/N cm). Second, the system of differential equations is replaced by the corresponding set of finite difference equations. For example, using the Crank-Nicholson scheme,
becomes
Fov(j + 1) - Fov(d = -Jo&i + 1) - Jov(~)
l/N
(6%
2
Equivalently,
IJo4.i + 1) + JovfAl = o
Fov(j + 1) - Fov(jj +
2N
1
(63)
where Few(j) is the volume flow in the proximal tubule at the point, x = j/N; and Jov(j) is the net volume flux out of the proximal tubule at the point, x = j/N. It is well known that the above discretization method is unconditionally stable and has second-order accuracy [22]. Similarly, discretized equations can be found for the remaining four differential equations in each nephron tubule and the nine differential and algebraic equations in the medullary interstitium. Hence, the general discretized system can be summarized as follows: f(x)
= 0,
(64)
where f is a vector representing the set of 34N discretized equations; and x is a vector representing the 34N discrete values of flows, concentrations, pressures, and potential. A solution to this problem
is found
by solving
the optimization
problem
m,‘n ifT(s)f(x).
(65)
In order to sofve the minimization numerically, the quadratic function, F(z) mated by the first two terms of its Taytor expansion
:= fTf,
is approxi-
where (67) Furthermore, the IMSL routine which was imp~emented~ DNEQBF, uses a trust region for the Taylor expansion of F(s), Hence, the problem becomes
where ~5is updated at each iterate. That is, we are going to trust the quadratic model (66) in a small, variable neighborhood of the point zC. In a qua&Newton
method, the update on x is given by
where Q is a doubfe-dogleg step; the Jacobian, J(x), is approximated by z1 finite difference; and S(r,) is the BFGS appro~mat~on to the Hessian matrix. For a detailed description of the numerical methods, see [ZO]. The initial conditions for each differential equation are given in Tabie 6, the parameters in Table 8, and the CI constants in Table 7, With these initial conditions and parameters, the system is ready to be solved using the quasi-Newton method outlined above. The tubule length is discretized into 20 intervals {i.e., N = 20). Hence, a system of 680 equations and 680 unknowns must be solved. Table 6. Initial conditions pressure.
for flow; sodium, patassium, and urea concentrations;
Initial Value
Independent Variable volume fhv,
1.0 (normalized)
F&(O)
sodium concentration, potassium
concentration,
urea concentration, P-=-=%
0.140 mmof/mL
C&(O)
0.0063 mmoI/mL
Cm(O)
0.0026 mmol/mL
C&4(0)
25.35 mmHg
fofOf Table 7. Cortical
interstitium constants. Value
Cortical Variable cortical sodium concentration,
Cc1
cortical potassium concentration, cortical urea concentration,
Cc4
ccs
0.140 mmol/mL 0.0045 mmol/mL 0.0026 mmol/mL
cortical pressure, PC
0.0 mmHg
cortical potential,
OmV
Qc
1
and
Solute and Water Transport Table 8. Parameter stitium.
values for each tubule of the nephron and the medullary
DHL
, ai1 I
79
0.41
%2 I %3
0.27
G,
0.15
G2
AHL
DT
0.0321
0.0127
0.0642
0.0127
0.15
0.075
0.15
0.075
CD
k/3
0.15
hi”
0.04
0.000814
0.0
0.000107
0.000357
‘&
1.4
0.058
0.284
0.0136
0.0217
hi2
1.4
0.058
0.048
0.0136
0.0217
x3
2.0
0.09
0.855
0.0041
0.0326
10.0
0.054
0.0389
0.0
0.0
h:4 Oil
1.00
0.96
1.0
1.0
1.0
Qi2
1.0
0.96
1.0
1.0
1.0
oi3
1.0
0.96
1.0
1.0
1.0
oi4
0.20
0.95
1.0
1.0
1.0
inter-
MI
%
0.0221
%3
0.0133
0;4
0.0353
uil
0.000827
ui3
0.001321
WI4 pow
0.0 25.35 VOLUME FLOW ALONG LENGTH OF NEPHRON
POSITION
Figure 3.
Due to the large size of the problem, an accurate initial guess is important. With a poor initial guess, the program will converge to a physiologically inaccurate solution. Because much care was taken to have this model be an extension of the model presented by Stephenson [ll], Stephenson’s solutions, averaged over each tubule, were used as a flat initial guess for this system. All computations were performed on an IBM RS/SOOO at the Center for Research in Scientific Computation using software packages written in Fortran and double precision arithmetic with 16 significant digits. A typical run time to solve the large optimization problem (65) is about two hours. Figures 3-7 present the solution in graphical format. All variables are graphed in the direction of flow.
N. E. KOTTLERet al.
CONCENTRATION
PT
1
DHL
OF SODIUM ALONG LENGTH OF NEPHRON
2
AHL 3 DT POSITION
4
CD
5
MI
Figure 4. CONCENTRATION
PT
1
DHL
OF POTASSlUM ALONG LENGTH OF NEPHRON
2
AHplos;nofT
4
CD
5
MI
6
Figure 5.
6. CONCLUSIONS Comparing the results of the model to the known behavior of the nephron, we notice many consistencies, along with a few inconsistencies. Sodium reabsorption in the proximal tubule is isotonic (i.e., sodium concentration is nearly constant along the proximal tubule); flow in the proximal tubule reduces to appro~mately 19% of its original value (slightly smaller than the accepted flow values which range from 20% to 35%); and solute concentrations and pressure reduce accurately. In general, the proximal tubule continues to behave accurately when incorporated into the rest of the model. Volume flow decreases and solute concentrations increase in the DHL, reflecting the passive extrusion of water into the medullary interstitium (counter-current multiplication). Sodium and chloride concentrations fall in AHL as expected, due to the active extrusion of salt, while volume flow remains constant, due to the impermeability of the AHL to water. Because the direction of flow in the AHL is toward the cortical interstitium, opposite that in other tubules, volume flow in the AHL is negative. In the DT and CD, volume flow reduces slightly and there is a net increase in the concentration of each solute. Pressure decreases accurately in
Solute and Water Transport
CONCENTRATION
81
OF UREA ALONG LENGTH OF NEPHRON
POSITION
Figure 6. PRESSURE
1
DHL
2
ALONG LENGTH OF NEPHRON
AHL 3 DT POSITION
4
CD
5
MI
6
Figure 7.
the DT and CD, as well. Flow in the MI is negative, reflecting that the direction of MI flow is toward the CI. Notice also the sodium concentration gradient which forms in the MI. This gradient is the driving force of the counter-current multiplier system, and its existence is essential in any accurate renal model. Furthermore, urea concentration increases almost continuously throughout the nephron. This behavior appropriately reflects one of the ultimate goals of the kidney, to create a waste highly concentrated in urea. Unfortunately, the physiological properties used to measure the accuracy of the results are only guidelines. With access to more exact flow values, the solution could be compared quantitatively. The results were also compared with Stephenson’s [17]. Again, because the two systems were normalized with respect to different flow values, his solutions could only represent directions to follow, rather than exact values that must be achieved. In general, the results compared favorably with Stephenson’s, with the following exceptions. There is a slight oscillatory behavior in the PT volume flow and in the sodium and potassium concentrations (we believe that these mild oscillations are numerical artifacts); Stephenson’s sodium concentration gradientin the MI
N. E. KOTTLER et al.
a2
is slightly more pronounced; we experienced a flat pressure profile and flat urea concentration gradients in the MI. The lack of a urea gradient, however, was also experienced by Stephenson [I?‘]. Therefore, this problem accurate
as a primary
improvement,
must be found.
initial
a more robust
Furthermore,
guess for this problem,
detailed
optimization
routine,
capable
of solving
renal flow data, which would present
must also be obtained.
Other
extensions
a more
to this research
include solving the parameter identification problem, assuming first constant, then nonconstant parameter values along each tubule; and implementing control to mimic the effects of both aldosterone and ADH. Hopefully, research in this area will be continued. what will someday ments
be a complete
this goal will be achieved
and accurate
At this time, we can only catch a glimpse
renal model.
However,
with continued
and will open the door to a more complete
of
improve-
understanding
of
renal physiology.
REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.
16. 17. 18. 19. 20. 21. 22.
J.J. Jacquez, B. Carnahan and P. Abbrecht, A model of the renal cortex and medulla, Mathematical Biosciences, An International Journal 1, 227-261 (1967). E. Koushanpour, R.R. Tarica and W.F. Stevens, Mathematical simulation of normal nephron function in rat and man, Journal of Theoretical Biology 31, 177-214 (1971). J.L. Stephenson, Transient behavior of the single loop cycling countercurrent multiplier, Bulletin of Mathematical Biology 35, 183-195 (1973). J.L. Stephenson, Concentrating engines and the kidney I. Central core model of the renal medulla, Biophysical Journal 13, 512-545 (1973). J.J. Jacquez and E. Daniels, Solute concentration in the kidney-I. A model of the renal medulla and its limit cases, Mathematical Biosciences, An International Journal 32, 307-335 (1976). D. Foster, J.J. Jacquez and E. Daniels, Solute concentration in the kidney-II. Input-output studies on a central core model, Mathematical Biosciences an International JountaZ32, 337-360 (1976). G.L. Barrett and J.S. Packer, Dynamic simulation of the renal medulla, Medical and Biological Engineering and Computing 21, 324-332 (983). A.S. Wexler, R.E. Kalaba and D.J. Marsh, Three-dimensional anatomy and renal concentrating mechanism. I. Modeling results, American Journal of Physiology 260, F368-F383 (1991). A.S. Wexler, R.E. Kalaba and D.J. Marsh, Three-dimensional anatomy and renal concentrating mechanism. II. Sensitivity results, American Jou7714l of Physiology 260, F384-F394 (1991). F.C. Hoppensteaclt and C.S. Peskin, Mathematics in Medicine and Life Sciences, Springer-Verlag, New York, (1992). J.L. Stephenson, Analysis of the transient behavior of kidney models, Bulletin of Mathematical Biology 40, 211-221 (1978). J.L. Stephenson, J.B. Garner and KS. Crump, Transient behavior of the single loop solute cycling model of the renal medulla, Bulletin of Mathematical Biology 40, 274-300 (1978). R. Mejia and J.L. Stephenson, Numerical solution of multinephron kidney equations, Jownal of Computational Physics 32, 235-246 (1979). J.B. Garner, R.B. Kellogg and J.L. Stephenson, Mathematical analysis of a model for the renal concentrating mechanism, Mathematical Biosciences, An International Journal 65, 125-150 (1983). R. Mejia and J.L. Stephenson, Solution of a multinephron, multisolute model of the mammalian kidney by Newton and continuation methods, Mathematical Biosciences, An International Journal 68, 279-298 (1984). J.L. Stephenson, Y. Zhang, A. Eftekhari and R. Tewarson, Electrolyte transport in a central core model of the renal medulla, American Journal of Physiology 253, F982-F997 (1987). J.L. Stephenson, Y. Zhang and R. Tewarson, Electrolyte, urea, and water transport in a two-nephron central core model of the renal medulla, American Journal of Physiology 257, F399-F413 (1989). D.E. Wessell, Parameter estimation for the proximal tubule in a single nephron steady state model of the kidney, Master’s Thesis, North Carolina State University, Raleigh, NC, (1993). IMSL, User’s Manual, IMSL Math/Libmry, FORTRAN Subroutines for Mathematical Applications, IMSL, Inc., Houston, TX, (1989). J.E. Dennis, Jr. and R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice-Hall, Englewood Cliffs, NJ, (1983). K.E. Brenan, S.L. Campbell and L.R. Petsold, Numerical Solution of InitiaG Value Problems in Dinerential-Algebraic Equations, SIAM Classics Series in Applied Mathematics, Philadelphia, PA, (1996). R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, New York, (1984).
Mathematical and Computer Modelling 29 (1999) 63-82
A Complete Steady State Model of Solute and Water Transport in the Kidney N. E. KOTTLER Department of Mathematics, Box 8205 North Carolina State University Raleigh, NC 27695-8205, U.S.A. Ii. T. TRAN* Center for Research in Scientific Computation, North Carolina State University Raleigh, NC 27695-8205 U.S.A. trana)control.math.ncsu.edu
Box 8205
D. E. WESSELL Department of Mathematics, Box 8205 North Carolina State University Raleigh, NC 27695-8205, U.S.A. (Received June 1998; accepted July 1998) Abstract-The purpose of this paper is to incorporate a detailed model, along with an optimized set of parameters for the proximal tubule, into J. L. Stephenson’s current central core model of the nephron. In this model a set of equations for the proximal tubule are combined with Stephenson’s equations for the remaining four tubules and interstitium, to form a complete nonlinear system of 34 ordinary differential and algebraic equations governing fluid and solute flow in the kidney. These equations are then discretised by the Crank-Nicholson scheme to form an algebraic system of nonlinear equations for the unknown concentrations, flows, hydrostatic pressure, and potentials. The resulting system is solved via factored secant update with a finite-difference approximation to the Jacobian. Finally, numerical simulations performed on the model showed that the modeled behavior approximates, in a general way, the physiological mechanisms of solvent and solute flow in the kidney. @ 1999 Elsevier Science Ltd. All rights reserved. KeywordsTransport model, Biomathemat~cs,Differential algebraic system, Inserve problem, Numerical methods.
1. INTRODUCTION In the late 1960s and early 1970s many researchers began to develop mathematical models of the kidney [l-4]. This research has progressed from the mid-1970s to this day with many variations, additions, and improvements [5-lo]. Some of these models have simply tried to fit the empirical data without regard to modeling the actual physiological processes. Researchers who were trying to develop physiologically accurate models encountered a substantial problem in modeling the mechanism that returns recycled water and salt from the kidney to the generaL circulation; this function is performed by a network of blood vessels in the kidney called the vasa recta. The central core model of the urine concentrating mechanism, which was initially proposed by Stephenson [4] *This research was supported in part by the National Aeronautics and Space Administration Grant NAGW-2865. 0895-7177/99/$ - see front matter @ 1999 Elsevier Science Ltd. All rights reserved. PII:SO895-7177(99)00055-2
Typeset by AM-T@
N. E. KOTTLERet
64
d.
in 1973, greatly simplified the modeling of the vasa recta while retaining all of the important physiological aspects of kidney function. This model quickly became the standard. In the past In 1987, Stephenson
20 years there have been many enhancements to the original model [ll-171.
et al. [16] considered for the first time electrolyte transport and introduced electric potential as a driving force for passive solute movement.
Other enhancements have been primarily due to
a greater understanding of renal physiology and to advances in computer technology that have allowed us to solve the model equations numerically [18]. As advanced as the model has become, however, there has been a disregard of the proximal tubule. In this paper, we extend Stephenson’s central core model of the nephron to include the proximal tubule. The physiologically based model contains a set of equations and parameters for all five nephron tubules, the medullary interstitium, and constants for the cortical interstitium. Because the model is an extension of that presented by Stephenson et al., parameter data is taken directly from Stephenson’s papers [16,17]. In addition, since the parameters for the proximal tubule were unavailable, they are estimated based on a nonlinear least-squares method with the cost function derived from the known behavior of the proximal tubule. The organization of the paper is as follows. We begin in Section 2 with a brief description of renal physiology. This summary of the renal physiology is important for the understanding of the mathematical model developed in subsequent sections. In Section 3, we first develop a set of differential and algebraic equations governing the flow of solvent and solutes in an arbitrary section of the nephron, then we customize these equations to each nephron tubule and to the interstitium. Section 4 presents our efforts on improving the model of the proximal tubule. The parameters used in the proximal tubule model are estimated based on a nonlinear least-squares method. The cost function used in the least-squares method is derived from the physiologically known behavior of the proximal tubule. Section 5 summarizes the steady-state nephron model of the kidney and discusses the numerical methods used to solve the nephron model and presents its numerical solution. Finally, Section 6 contains our concluding remarks including a comp~~on of our solution with the known behavior of the kidney.
2. THE KIDNEY The functional unit of the kidney, the nephron, is composed of a glomerulus and five connected tubules, the proximal tubule (PT), the descending limb of the loop of Henle (DHL), the ascending limb of the loop of Henle (AHL), the distal tubule (DT), and the collecting duct (CD). This general structure of the nephron is depicted in Figure 1. The tubules form a continuous system which modifies the filtrate and transports it to a cavity in the kidney called the renal pelvis. This filtrate which passes through the glomerulus is referred to as glomerular ultrafiltrate (GU). The combined effects of each tubule enable the kidney to form a concentrated waste product while maintaining homeostasis in the body-the primary function of the kidney. Each tubule has its own role integral to the success of this process. 2.1. Interstitium There are two interstitia within which the nephron lies, the cortical interstitium (CI), and the medullary interstitium (MI). While the MI has a set of equations modeling its compIex behavior, the CI is simply assumed to be a well-mixed bath of constant volume. The PT and DT lie in the CI; the DHL, AHL, and CD lie in the MI. 2.2. Proximal
Tubule
The PT acts as a large sieve, returning the majority of what it is delivered from the glomerulus back to the general circulation. This is accomplished by a sodium potassium pump which actively transports sodium and potassium ions into the CI. The active transport favors a passive movement of chloride ions (via the potential gradient) and water (via the osmotic gradient) into the CI.
SoIute and Water T’rsnsport
Figure 1. General structure
Because
the net flux of PIat- is matched
by water,
65
af the nephron.
sodium
absorption
is isotonic
along the length
of the PT, Movement of urea is primarily due to the effects of its concentration gradient and solvent drag. At the end of the proximal tubule, approximately 70% of the Hz0 and Na+, 85% of the K+, 42% of the Cl-, and 50% of the urea that entered the proximal tubule have been returned to the general circulation via the cortical interstitium. The goal of forming a concentrated waste product (i.e., urine) has already begun in the proximal tubule. The function af the remaining tubules of the nephron is to enhance this concentrating process. 2.3.
The Descending Limb of Henle’s Loop
Active transport of Na+ in the AHL causes the formation of a sodium in the MI. As the GU traverses down the DHL, the exterior environment more concentrated. In fact, by the bottom of the DHL, in the deepest concentration levels in the MI reach three to four times what they are The walls of the DHL are relatively impermeable to the passive diffusion very smalf amounts of Na+ and Cl- enter the DHL to equilibrate this
~oneentrati~~ gradient becomes increasingly region of the meduila, in the upper medulla. of NaCl. Hence, only concentration gradient.
The DHL, however, is very permeable to water. Thus water is drawn out of the DHL by osmosis and enters the capillary network to return to the vascular system. There is no active transport in the DHL; and solvent drag, concentration gradients, and electrical gradients lead only to minor movements of Nai, Cl-, and IS+. Simil~ly~ movement of urea is affected by solvent drag and the urea concentration gradient. In essence, the descending limb of Henle’s Ioop accomplishes two main tasks: the extrusion of an additional 50% of the water, and hence, the equilibration of the concentration gradient between the GU and the meduliary interstitium. Thus, by the end of the descending limb the tubular fluid is isosmotic with the exterior enviranment, and a highly concentrated filtrate is delivered to the ascending limb of Henle’s loop. 2.4. The Ascending Limb of Henle’s Loop The system most commonly regarded nephron, counter-current multiplication,
as the backbone of the urine ~on&entrat~ng ability of the resides in the loop of Henle. Saft is actively pumped
66
N. E. KOTTLER
et al.
out of the AHL creating a concentration gradient in the MI. As the filtrate in the DHL passes through the concentrated MI, water flows out of the filtrate by osmosis. Hence, the DHL delivers a solute concentrated filtrate to the AHL, allowing for the extrusion of even more salt, and so on. Protecting the concentration gradient in the MI, the walls of the AHL are impermeable to water. Therefore flow remains constant in the AHL. In addition, because there is no fiux of water out of the AHL, there is no solvent drag. Thus, urea is moved only by the urea concentration
gradient.
In essence, the counter-current mechanism traps large amounts of sodium in the medullary interstitium, enabling large amounts of water to be reabsorbed in the descending limb. It also lowers the osmolality
of the filtrate
as it moves up the ascending
limb toward
the distal convoluted
tubule. 2.5. The Distal Convoluted As a result
of the counter-current
Tubule multiplier
system,
the osmolality
of the GU is slightly
lower
than that of the CI, causing water to leave the DT by osmosis. The two forces moving urea are the urea concentration gradient and solvent drag. The reabsorption of Na+ and K+ is controlled by the hormone aldosterone. Aldosterone is secreted by the adrenal cortex in response to a high K+ concentration in the blood. An increase in aldosterone secretion yields an increase in Naf reabsorption (out of the DT), and an increase in K+ secretion (into the distal tubule). In fact, without any aldosterone secretion, 100% of the K+ is reabsorbed into the vascular system; whereas with aldosterone secretion at its highest, 100% of the Naf is reabsorbed. Thus, aldosterone can play a powerful role in maintaining homeostasis. For the purposes of this paper, an average aldosterone secretion, which would occur under hom~static conditions, is assumed. For a more complete and accurate model, a control system representing aldosterone secretion must be incorporated. 2.6.
The Collecting Duct
As a result of the counter-current multiplier system, the exterior environment becomes increasingly hypertonic as the filtrate traverses down the CD. In addition, the counter-current multiplier system has caused the fluid entering the CD to be hypotonic (despite the efforts of the DT). Therefore, counter-current multiplication creates a large concentration gradient across the wall of the CD. Salt does not enter the CD to relieve the concentration gradient because the walls of the CD are impermeable to salt. Furthermore, there is no active transport of any solute in the CD. Hence, strong osmotic effects cause water to be drawn into the MI, where it enters the capillary network and is returned to the general circulation. The rate at which water is reabsorbed is controlled by ADH, antidiuretic hormone. ADH, produced by the neurons in the hypothalamus and secreted by the posterior pituitary gland, increases the permeability of the collecting duct to water. For example, in a state of severe dehydration ADH is secreted, and up to 99.8% of the original GU is reabsorbed (the remaining 0.2% is the obligatory water loss). In a well hydrated person, ADH would not be secreted, and 99.2% of the original GU is reabsorbed. This control translates to urine volumes ranging from as small as 0.4L (dehydrated person), to 1.5 L (well hydrated person). In this paper, an average ADH concentration is assumed. Therefore, while this model may be accurate under homeostatic conditions, it will not accurately represent any deviations. As an expansion of this model, an ADH control system should be modeled. At this point, the processing by the kidney is complete. The urine flows from the collecting duct via the ureters to the bladder. The entire process is summarized in Tables 1 and 2. 3.
MATHEMATICAL
DESCFUPTION
The solutes we will model are Na’, Cl-, K+, and urea. Figure 2 shows a section of a genera1 tubule of the nephron. The mass balance equation for solute k in the differential volume element
67
Solute and Water Transport Table 1. Relative flow rates of water and relative amounts of sodium, urea, and potassium in different segments of the nephron of a rat kidney.
I
Nephron Segments Gtomerulus
Percentage Remaining Potassium Sodium Urea
Water 190
1ilQ
100
IQ6
End of Proximal Tubule End of Descending Limb
23.0 17.0
23.0 47.1
44.7 605.3
End of Ascending Limb
17.0
9.5
133.0
3.7
2.9 0.1
77.1 26.3
26.7 11.9
End of Distal Tubule End of Collecting Duct
6.0 0.24
16.4 86.9
Table 2. Concentrations of sodium, urea and potassium in different segments of the nephron of a rat kidney.
I
Concentration in mmol,/L
Nephron Segments
Sodium
Urea
Potassium
Gfomerulus End of Proximal Tubule
140.0 140.0
3.0 5.0
6.3 4.5
End of Descending Limb
294.4
31.8
10.3
End of Ascending Limb End of Distal Tubule
101.0 89.7
35.4 87.8
3.4 8.3
End of Collecting Duct
176.2
243.3
71.8
J J
fP f+
i
t
1
)
!
t,
*
!
x,=2
i f !
a
X,=L,
Figure 2. A general flow tubule. of tubule i shown in Figure
2 is given by
where AZ; is the length of the dierentiaf vdume element, Ai is the cross-sectional area of the differentid volume element, C:kfzi, t> is the ~~n~~nt~at~~nof solute k in tubule i, Fikfz+, t) is
N.
68
the axial solute flow of solute tubule
destroyed in tubule The corresponding
i, Jik is the net outward
k in tubule
i, and Sik is the average
E. KOTTLER et al.
net rate per unit volume
transmural
at which solute
flux of solute
k in
k is being created
or
i. mass balance d [AziAi]
for volume
(i.e., water)
flow in tubule
i is given by
= Fiv (xi, t) - Fiv (xi + A~ir t) - Axi Jiv,
at where Fiv is the axial volume
equation
i, and Jiv is the transmural
flow in tubule
t) with its Taylor series expansion,
In both cases, replacing the term F~,(xi+Az:i, and letting Axi ---f 0, yields, respectively,
a[Ad& (Xi,t)] = _ aFik(xi7t, _ Jik + at and
volume
&ik
aFiv (xi,t) ax
Equation (3) gives a partial differential equation arbitrary flow tubule. Equation (4) gives a partial in an arbitrary flow tubule.
dividing
i.
by Axi,
(3)
dX
aAi -=at
flux in tubule
_ J,
ZV.
for mass balance of an arbitrary solute in an differential equation for mass balance of water
For the purposes of this paper, we are interested only in the steady state problem, and thus, the left-hand sides of both equations (3) and (4) become zero. In addition, we make the simplifying assumption that no solute is being produced or destroyed in any tubule (i.e., sik = 0). Hence, the mass balance equations simplify to
(Xi,t> = _ J,
a&k
tk
ax and
8FiV (xilt) = _ ax
The transmural pressure gradients written as
J,
zv.
(6)
volume flux, Jiv, is determined by a combination that exist between the tubule and the interstitium.
Jiv = hiv
2
RT (elk - Cik) gik f Pi - PI
k=l
of the concentration and In general Jiv can be
1,
(7)
where hiv is the hydraulic permeability coefficient per unit length for tubule i, R is the universal gas constant, T is the absolute temperature, CIk is the concentration of solute k in the interstitium surrounding tubule i, C& is the concentration of solute k in tubule i, oik is the Staverman reflection coefficient of the wall of the ith tubule to solute k, Pi is the pressure in tubule i, and PI is the pressure in the interstitium surrounding tubule i. The transmural solute flux, Jik, is determined by a combination of solvent drag, Fick’s Law, electronic attraction, and active transport. Hence, Jik can be written as
Jik = J~v (1 -
aik>
CIk
i- Cik n L
+
hik (CIk - Cik) +
‘Ik; ‘jkg
(‘Pi -
e,,] +
(8)
J$,
where hik is the passive permeability of tubule i to solute k, zk is the valence of solute k, F is Faraday’s constant, 8i is the potential in tubule i, *I is the potential in the interstitium surrounding tubule i, and J,“kis the active transport of solute k from tubule i.
69
Solute and Water Transport
The active transport, J$, is a metabolically Menten kinetics. that is.
driven transport that is consistent with Michaelis-
where o;k is the m~mum rate of active transport, and &k is the ~on~entratian at which 2% reaches onehalf its m~mum value (Le., i%&is the Michaelis constant). The Nernst-Plank equation relates solute flow, Fit,, to v&me flow, &v, concentration, C&-l and potential !3Zrias follows:
where L&k is the diffusion coefficient of solute k in tubule a, and T&kis the mobility of solute k in tubule i. Furthermore, diffusivity and mobility are related by
Assuming Poiseuille flow, pressure drop along tubule i is given by the equation
where Pi is the hydrostatic pressure in tubule i; Rg = 8qfrf is the hydrodyn~~~ flow resiistance of tubule i; where 77is the viscosity of the solvent (water) and T is the radius of the cylindrical tubule. Since Na+, K+, and Cl- are ions, electroneutrality must be maintained. This condition yields the final equations, 3 c k=l
and
=
0
(13)
3 c k=l
Using
zkC$k
equation
(f3), the concentration
Zks&k =
0.
(14)
of Cf-, a negatively charged ion, can be written in
terms of the ~on~ent~at~on~ of Na+ and K+, the positively charged ions. Hence,
where Cis is the concentration of chloride in tubule i; C+r is the concentration of sodium in tubule i, and Ci3 is the concentration of potassium in tubule i. Finally, substituting equation (8) into equation (I$) and solving for 9i - 91 yieids
i&.-g,
=
5 Zk f&V (1 - %k) ((elk + cik) 12) f hik (elk 5 +ik ((Clk +Cik) 12)(F/W
Cik) + J$]
_k=l
k=l
Our general system can be summarized as follows: in each tubule, i, we have eleven dependent variables: volume flow, Fgv; sodium flow, 2%; sodium concentrations Gil; chloride flow, Fi2; chloride concentration, i&; potassium Bow, S’i3; potassium concentration, CS; urea flow, F&; urea concentrations C&; pressure, Pi; and potential, 96. These variables are matched by nine differential equations: (5), (6), (lo), and (12); and twa algebraic equations: (15) and (IS).
70
N. E. KOTTLER
&al.
In all tubules except the medullary interstitium, we assume that the diffusivity and mobility terms are negligible as compared to the convective flow term. Thus, in these tubules, equation (10) simplifies to the algebraic relationship,
Fik = FivCik. Taking the derivative on both sides and solving for 2
(17) yields
1
-=-
(18)
.
Substituting equations (6) and (5) into (18) gives an equation for the change in concentration of each solute, k. Namely,
JiVCi/c- Jik Fiv ’
acik -= ax
(19)
Once the four solute concentrations are found, equation (17) gives values for the four solute flows, Fi1, Fis, Fi3, Fi4. Since the solute flows can be computed algebraically, the four equations represented by differential equation (5) do not need to be evaluated. Hence, the total number of differential equations in each tubule, except for the medullary interstitium, is reduced from nine to five. These are: -
the volume flow equation a&v -
ax
-
ax
(21)
JivCi3- Ji3 Fiv ’
(22)
the urea concentration equation aCi4 -=
ax
-
JivCil - Jil Fiv ’
the potassium concentration equation aCi3 -= ax
-
(20)
the sodium concentration equation acil -=
-
= -Jiv,
JivCi4
-
Fiv
Ji4
(23)
’
the pressure equation
(24) The remaining values to be computed are chloride concentration, C~Z and potential, +. The chloride concentration, CQ, can be found from equation (15); and Qi can be computed using the algebraic relationship expressed in equation (16). Now we can adapt our general system of equations to each tubule by determining the appropriate forms for Jiv and Jik. 3.1. The Proximal Tubule Model Using equation (7), we can write the equation for volume flux in the proximal tubule,
Jov = hov f: RT (Cck - COk) ffOk + PO- PC , k=l
I
(25)
71
Solute and Water Transport
where CCk is the concentration of solute k in the cortical interstitium, and PC is the pressure in the cortical interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium. cc1
-t Co1
Joi = Jov (1 - COk)
(Cc1
2
-co1>+
“I
+ 2 ‘01
g
(%I
-
Ic)]
+
J&;
(26)
cc3
+ ‘03 2
g
(go
-
qc)]
+
Jg3;
(27)
for potassium,
Jo3
=
Jov
(1 -
cc3
ao3)
+ co3 2
and for urea, Jo4
=
Jov
(1. -
(Cc4
+ ho4
cc4
co4)
-
+ co4 2
Co4)
f
G4
+ G4 2
g
(ipo
-
Bc)]
_
[
Note that the flux equation for urea does not contain an active transport actively transported in the proximal tubule.
term since urea is not
3.2. The Descending Limb of Henle’s Loop Model Using equation (7), we can write the equation for volume flux in the descending limb,
Jlv = hv
&WIX
-
Clk)fllk
+ pl
1 ,
-PI
k=l
(2%
where Cfk is the concentration of solute k in the medullary interstitium; and PI is the pressure in the medullary interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium, Jll
=
&V
(1 -
Cl1
elk)
+ Cl1 ,, L
(30)
for potassium, A3
=
JlV
(1 -
CI3
613)
+ Cl3 2
(31) (C13
+ h13
-
G3)
+
cI3
+2 c13
g
(i&l
_
~~)I
;
+2 (24
g
(\El
_
*,,I
.
[
and for urea, J14 = JIV (1 -
+ hi4
CI4
~14)
((714
-+- Cl4
2 -
64)
(32) +
a4
Note that none of the flux equations contains an active transport transport in the descending limb.
term since there is no active
N. E. KOTTLER et al.
72 3.3.
The Ascending
Limb
of Henle’s
Loop
Model
Since the ascending limb is not permeable to water, there is no volume flux into or out of the ascending limb. Hence, the equation for volume flux is simply, (33)
Jzv = 0. The equations for solute flux in the ascending limb are, for sodium,
+.
J21 = h21 (Cr1 - C21) +
J& ;
(34)
for potassium, J23 = h23 (C13 - C23) + ‘I3 + c23 g 2
(‘P2 - QI)] + J2a33;
and for urea, - C24) + c14 +2 c24 g
J24 = h24 (64
(q2 - I,)]
.
(36)
Note that the flux equations have no solvent drag term, because J2v = 0. Also, the flux equation for urea does not contain an active transport term since urea is not actively transported in the ascending limb. 3.4. The Distal
Tubule
Model
The equation for volume flux in the distal tubule can be expressed as
J3v = h3v
f:
RT (Ccrc -
C3k)
c3k
+
p3 - PC
k=l
1,
(37)
where 63k = 1, V k. Hence, water is a strong force equilibrating the concentration gradient which exists between the distal tubule and the exterior environment, the cortical interstitium. In a similar manner, we can use equation (8) to write the equations for solute flux. These are, for sodium, J31 = h31 (Cc1 - C31) + cc1 + c31 g 2
(‘P3 - SC)]
+ J&l;
for potassium, J33 = h33 (Cc3 - C33) + co
+ c33 g 2
(\k3- Qc)]
;
(39)
and for urea, J34 = 0.
(40)
Notice that none of the flux equations contains a solvent drag term. The reason for their absence is that the reflection coefficients for all of the solutes in the distal tubule are at their maximum value (i.e., U3k = 1, Vk). Also, there is no active transport of potassium or urea in the distal tubule; thus the absence of the Michaelis-Menton terms in those equations. Finally, the distal tubule is impermeable to urea. Combining this fact with the lack of urea transport and with the reflection of urea at its maximum, there is no net flux of urea into or out of the distal tubule.
73
Solute and Water Transport
3.5.
The Collecting
Duct
Model
The equation for volume flux in the collecting duct can be expressed as
J4v
= h4v
2
RT
(CI~
-
C4k) g4k + p4 -
k=l
PI
1,
(41)
where U& = 1, V k. Similarly, we can write the equations for solute flux. These are, for sodium,
J41 = h41
(Crl
J43 = h43
(C13 -
C41) + cIl
-
+2 c41 g
(Q4 -
*I)]
;
+2 c43 g
(Q4 -
QI)]
;
for potassium, + c13
C43)
(43)
and for urea, J44 = 0.
(44)
Similar to the distal tubule, the reflection coefficients in the collecting duct are at their maximum. Hence, solvent drag is not a factor as all of the solutes are reflected back into the tubule at the walls of the collecting duct. Also, there is no active transport of any solute in the collecting duct, including sodium. Finally, the collecting duct is relatively impermeable to urea. Thus, ail of the urea delivered to the collecting duct remains inside and is transported as waste out of the kidney. 3.6. The Medullary
Interstitium
Model
Diffusivity and mobility terms can no longer be considered negligible as compared to convective flow in the medullary interstitium. Therefore, we no longer have an algebraic relationship between solute flow and solute concentration. However, we know that in the medullary interstitium, the sum of the solute flows at any level is equal to the final excretion of that solute. In other words, we can use the following equation to compute solute flows in the medullary interstitium:
c
Fik(z)
= F4k(N),
i
where i is summed over all tubules which lie in the medullary interstitium (i.e., i = I, 1,2,4); Fik is the flow of solute k in tubule i; N is the final discretized position in any tubule; and Fdk(N) is the final flow of solute k in the collecting duct. We can use equation (17) to compute the solute flows in each tubule except for the interstitium. Incorporating these into equation (45) yields solute flow equations for sodium, FII(~)
+ F~v(~)Cll(z)
+ F2v(z)C21(z)
+ F4v(z)C41(2)
-
F4v(N)C41(N)
= o;
(46)
+ Flv(z)Cl3(z)
+ FZV(z)C23(z)
+ F4v(x)C43(5)
-
F4v(N)C43(N)
= 0;
(47)
+ F1v(x)C14(z)
+ FZV(z)C24(5)
+ F4v(z)C44(2)
-
F4v(N)C44(N)
= 0.
(48)
for potassium, FIN
and for urea, F14(2)
We have a similar condition for solvent flow in the medullary interstitium. Namely, the sum of the volume flows at any level of the medullary interstitium is equal to the final flow in the
N. E. KOTTLER
74
et al.
collecting duct. This condition translates to the following algebraic equation for solvent flow in the meduliary interstitium: E&V(X)
=
F4V(N).
(49)
i
Equivalently, FlV(Z) + Applying the electroneutrality
FlV(X)
+ J72v(s)
+ F4v@)
-
iJ4v(N)
=
0.
(56)
condition (14) to the medullary interstitium, we have our final
equation, QJrk = 6,
e
(511
k=l
where Jlk: = - E
Jirc
(52)
and Jik is the net transmural flux out of each tubuk which lies in the medullary interstitium (i.e., i = 1,2,4). Hence, ZI (-JII
- J21 - J41) +
22
c-J12
- J2z - J42) + 23 t-513
- 523 - J43f = 0.
(53)
However, the flux of chloride out of tubule i was never explicitly calculated because chloride concentration and hence, flow, could be algebraically computed in all other tubules by the electroneutrality condition, (13). Therefore, we need an equation governing the net flux of chloride out of tubule i. This equation comes from equation (5), where Fi2 = FivCi2 in all tubules except for the interstitium. Thus we have WvG2 3:
=3:_ J, 22.
(54)
Differentiating and using the product rule, we obtain
Equation (55) holds in all tubules except for the interstitium. Hence, we can compute the flux of chloride out of the descending limb, Ji2, the ascending limb, J22, and the collecting duct, 542. Using these net chloride fluxes, we can complete equation (53). The remaining four equations for pressure and solute concentrations are as previously stated. However, it should be once again noted that the entire Nernst-Plank equation must be used to compute solute concentrations in the medullary interstitium. 3.7. The Cortical
Interstitium
Constants
Keeping in line with Stephenson’s model, we will assume that the cortical interstitium is a well-mixed bath of constant volume. Thus, concentration of each solute, (Ccl, CCZ, Cca, and Cc*); pressure, (PC); and potential, (tkc), are all constant in the cortical interstitium. Actual values for these constants can be found in Table 7. 4.
PARAMETERS
AND
BOUNDARY
CONDITIONS
In order to complete our model, initial conditions for the five independent variables (volume flow; sodium, potassium, and urea concentrations; and pressure) as well as parameter values are required. When available, these initial conditions and parameters were taken from Stephenson’s papers [16,17]. However, because Stephenson did not model the PT, several modifications and parameter identifications were necessary. Specifically, Stephenson’s initial conditions are those for the DHL. In order to get our initial conditions to match Stephenson’s at the PT-DHL boundary,
Solute and Water Transport
75
we extrapolated his values based on the physiology of the PT (see Tables 1 and 2). For example, Stephenson’s input volume flow into the descending limb is lO.OnL/min.
Table 1 indicates that
this value is 23% of the initial flow in the proximal tubule (in the literature, values range from 20% to 35%). Therefore, a flow of approximately 33 nL/min is used as the input volume flow for this problem. We also note that the initial pressure is a parameter (its calculation will be discussed subsequently).
Table 6 presents the initial conditions for each differential equation. The system
parameters and their values are listed in Table 8. In order to remedy problems associated with varying tubule length and poor scaling (volume flows are on the order of 10-i’ (L/s) while pressure values are on the order of 10 (mmHg)), the system parameters are normalized with respect to tubule length and entering volume flow as follows (note: Z’ denotes the normalized parameter x). For solute permeabilities,
hi/c&
hk= piv(o)
(56)
’
where hik is the nonnormalized permeability of tubule i to solute k; Li is the length of tubule i; and Fiv(O) M 5.5 x lop7 mL/sec is the entering flow rate in the proximal tubule. For hydraulic permeability, hi,
=
$-$
(57)
EV
where hiv is the nonnormalized hydraulic permeability of tubule i. For maximum active transport rates, I aik
=
aikLi ~~~(0)
(58)
’ .
where aik is the nonnormalized maximum active transport rate of solute k in tubule i. For flow resistance, Rip = 7.5062 X 10m4qLiFiv(0),
(59)
where RiF is the nonnormalized flow resistance of tubule i; and v is the viscosity of water in poise (note: 7.5062 x 10m4 is a conversion factor for converting poise to mmHg s). It is necessary to normalize neither the reflection coefficients, gik, which are dimensionless, nor the Michaelis constants, kik, which are independent of flow and tubule length. Nonnormalized parameter values for the DHL, AHL, DT, CD, and MI were taken directly from Stephenson [16,17]. However, no micropuncture data was available for PT parameters. Hence, an ‘optimal’ set of parameters for the PT was determined independently of the other four tubules and MI. This optimization was done so that the PT model behaved as closely as possible to the known behavior of the PT. The parameter estimation problem was formulated as follows: let Q = (o&i,koi,abs, ko3, hb,, h&, hbar hb4,uo4, PO(O)) be the PT parameter set. In order to limit the number of parameters being estimated, we do not estimate the reflection coefficients for Na+, Cl-, or K+; ph ysiological evidence shows these parameters are close to one [14]. In addition, the initial pressure is listed as a parameter. This inclusion is due to the fact that initial pressure varies greatly in the proximal tubule (range 20-40 mmHg); however, the final pressure should be very near 10.75mmHg for a smooth transition to the DHL. Thus, final pressure is taken into account in our cost function and we allow initial pressure to be varied. The parameter estimation problem is to find the set of parameters 4 that minimizes
J (q) = wi fl: [Co1 (Zj) - 0.14012 +
w2
[Fov
(ZN)
-
0.312
j=l + w3 [CO3 (ZN)
-
o.oo45]2 + w4 [Co4 (ZJV) - 0.00512
+ ws [PO(ZN) - 10.7512,
(60)
76
N. E. KOTTLER et al. Table 3. Weights for cost function J. Value
Weight constant sodium concentration
weight, wi
2.0
final volume flow weight, ws final potassium
1.0
cont. weight, w3
1000.0 1000.0
final urea cont. weight, w4 final pressure weight, ws
0.001
Table 4. Parameter bounds for the proximal tubule model. Parameter active transport
Lower Bound
of sodium, a&
._
Upper Bound
0.0
10.0
Michaelis constant for sodium, kor
0.0
10.0
active transport
0.0
10.0
0.0
10.0
of potassium,
a&s
Michaelis constant for potassium, hydraulic permeability,
h&,
sodium & chloride permeability, potassium permeability, urea permeability,
ko3
hb, , hb,
hbs
hb4
urea reflection coefficient,
of34
initial pressure, PO(O),
0.0
1.0
0.0
100.0
0.0
100.0
0.0
100.0
0.0
1.0
20.0
40.0
Table 5. Initial and final parameter values for the proximal tubule model. Parameter active transport
Initial Guess
of sodium, a&
Final Value
0.10
0.41
Michaelis constant for sodium, kol
0.15
0.15
active transport
0.10
0.27
of potassium,
Michaelis constant
for potassium,
hydraulic permeability, permeability,
urea permeability,
h&s
hb4
urea reflection coefficient, initial pressure, PO(O) value of the cost function numbers of iterations
ko3
hb,
sodium & chloride permeability, potassium
a&
004
h& , h&
0.15
0.15
0.001
0.04
1.5
1.4
1.5
2.0
1.5
10.0
0.50 35.0 0.345525
0.20 25.35 0.2968383-
04
22
where N is the final mesh point and Cor(zj), Fcv(srv), Cc3(5~), COG, and Po(z:N) satisfy the system differential equations (20)-(24) d escribing the proximal tubule. The summation term in equation (60) requires that sodium concentration remain constant at O.l40mmol/mL; final volume flow be 30% of its initial value; final potassium and urea concentrations be, respectively, 0.0045 and 0.005 mmol/L; and final pressure be 10.75 mmHg. The weights, 2uj are chosen so that a I% error in any of the five terms adds an equal amount to the cost function, J. The appropriate weights are listed in Table 3. The minimization of the cost function J is done using a finite difference quasi-Newton method with an active set strategy since the problem is one with simple bounds. This scheme is implemented via the IMSL routine DBCONF 1191.For a detailed description of the numerical methods, see [ZO]. The parameters in the vector 4 are physical quantities that are bounded below by zero (except PO) and above by different values as listed in Table 4. The optimized parameter values (in the column “Final Value”) are given in Table 5 along with the initial value and final value of the cost function J. The initial guesses of these parameters
Solute
are taken
to be average
determined
values
by Stephenson
and Water
of the similar
parameters
[16]. We note that
77
Transport
in other
the optimized
segments
parameter
values
of the cost function by an order of magnitude of 4. These final values for the PT parameters are used as approximations in our larger system.
5. THE
NEPHRON
The complete
steady-state
MODEL
AND
34 ordinary differential equations and algebraic equations (DAE) is summarized below.
(0
For each of the five nephron
tubules,
of the kidney
constraints.
there
as
the value
SOLUTION
consists
The system
are five ordinary
reduced
for those PT parameters
ITS NUMERICAL
model of a single nephron
of the nephron
of a coupled
of differential
differential
set of
algebraic
equations,
equa-
of volume (i.e., water) flow, sodium concentrations (20)-(24), d escribing the distributions tion, potassium concentration, urea concentration, and the pressure. In these equations, the expressions for the transmural volume and solute fluxes, J~v and Jik, are given in Sections 3.1-3.5. the concentrations of sodium, potassium, and urea are solved (ii) In the medullary interstitium, from the Nernst-Plank equation (10). In addition, the volume flow and solute (sodium, potassium, and urea) flows satisfy the algebraic constraints (46)-(49). The pressure still satisfies the pressure equation (24) and the transmural solute flux in the medularry interstitium, JIk, satisfies the electroneutr~ity condition (51) or equivalently (53). The total system of 34 differential algebraic equations (29 differenti~ equations and five algebraic equations) is solved numerically by a quasi-Newton method, implemented in the IMSL routine, DNEQBF [19]. For the analysis of the numerical method implemented here and other possible numerical methods for DAE, we refer the reader to Brenan et al. [21] (and the numerous references therein). First, each tubule is divided into N segments. Because our system is normalized with respect to tubule length, each discretized segment has the same length (l/N cm). Second, the system of differential equations is replaced by the corresponding set of finite difference equations. For example, using the Crank-Nicholson scheme,
becomes
Fov(j + 1) - Fov(d = -Jo&i + 1) - Jov(~)
l/N
(6%
2
Equivalently,
IJo4.i + 1) + JovfAl = o
Fov(j + 1) - Fov(jj +
2N
1
(63)
where Few(j) is the volume flow in the proximal tubule at the point, x = j/N; and Jov(j) is the net volume flux out of the proximal tubule at the point, x = j/N. It is well known that the above discretization method is unconditionally stable and has second-order accuracy [22]. Similarly, discretized equations can be found for the remaining four differential equations in each nephron tubule and the nine differential and algebraic equations in the medullary interstitium. Hence, the general discretized system can be summarized as follows: f(x)
= 0,
(64)
where f is a vector representing the set of 34N discretized equations; and x is a vector representing the 34N discrete values of flows, concentrations, pressures, and potential. A solution to this problem
is found
by solving
the optimization
problem
m,‘n ifT(s)f(x).
(65)
In order to sofve the minimization numerically, the quadratic function, F(z) mated by the first two terms of its Taytor expansion
:= fTf,
is approxi-
where (67) Furthermore, the IMSL routine which was imp~emented~ DNEQBF, uses a trust region for the Taylor expansion of F(s), Hence, the problem becomes
where ~5is updated at each iterate. That is, we are going to trust the quadratic model (66) in a small, variable neighborhood of the point zC. In a qua&Newton
method, the update on x is given by
where Q is a doubfe-dogleg step; the Jacobian, J(x), is approximated by z1 finite difference; and S(r,) is the BFGS appro~mat~on to the Hessian matrix. For a detailed description of the numerical methods, see [ZO]. The initial conditions for each differential equation are given in Tabie 6, the parameters in Table 8, and the CI constants in Table 7, With these initial conditions and parameters, the system is ready to be solved using the quasi-Newton method outlined above. The tubule length is discretized into 20 intervals {i.e., N = 20). Hence, a system of 680 equations and 680 unknowns must be solved. Table 6. Initial conditions pressure.
for flow; sodium, patassium, and urea concentrations;
Initial Value
Independent Variable volume fhv,
1.0 (normalized)
F&(O)
sodium concentration, potassium
concentration,
urea concentration, P-=-=%
0.140 mmof/mL
C&(O)
0.0063 mmoI/mL
Cm(O)
0.0026 mmol/mL
C&4(0)
25.35 mmHg
fofOf Table 7. Cortical
interstitium constants. Value
Cortical Variable cortical sodium concentration,
Cc1
cortical potassium concentration, cortical urea concentration,
Cc4
ccs
0.140 mmol/mL 0.0045 mmol/mL 0.0026 mmol/mL
cortical pressure, PC
0.0 mmHg
cortical potential,
OmV
Qc
1
and
Solute and Water Transport Table 8. Parameter stitium.
values for each tubule of the nephron and the medullary
DHL
, ai1 I
79
0.41
%2 I %3
0.27
G,
0.15
G2
AHL
DT
0.0321
0.0127
0.0642
0.0127
0.15
0.075
0.15
0.075
CD
k/3
0.15
hi”
0.04
0.000814
0.0
0.000107
0.000357
‘&
1.4
0.058
0.284
0.0136
0.0217
hi2
1.4
0.058
0.048
0.0136
0.0217
x3
2.0
0.09
0.855
0.0041
0.0326
10.0
0.054
0.0389
0.0
0.0
h:4 Oil
1.00
0.96
1.0
1.0
1.0
Qi2
1.0
0.96
1.0
1.0
1.0
oi3
1.0
0.96
1.0
1.0
1.0
oi4
0.20
0.95
1.0
1.0
1.0
inter-
MI
%
0.0221
%3
0.0133
0;4
0.0353
uil
0.000827
ui3
0.001321
WI4 pow
0.0 25.35 VOLUME FLOW ALONG LENGTH OF NEPHRON
POSITION
Figure 3.
Due to the large size of the problem, an accurate initial guess is important. With a poor initial guess, the program will converge to a physiologically inaccurate solution. Because much care was taken to have this model be an extension of the model presented by Stephenson [ll], Stephenson’s solutions, averaged over each tubule, were used as a flat initial guess for this system. All computations were performed on an IBM RS/SOOO at the Center for Research in Scientific Computation using software packages written in Fortran and double precision arithmetic with 16 significant digits. A typical run time to solve the large optimization problem (65) is about two hours. Figures 3-7 present the solution in graphical format. All variables are graphed in the direction of flow.
N. E. KOTTLERet al.
CONCENTRATION
PT
1
DHL
OF SODIUM ALONG LENGTH OF NEPHRON
2
AHL 3 DT POSITION
4
CD
5
MI
Figure 4. CONCENTRATION
PT
1
DHL
OF POTASSlUM ALONG LENGTH OF NEPHRON
2
AHplos;nofT
4
CD
5
MI
6
Figure 5.
6. CONCLUSIONS Comparing the results of the model to the known behavior of the nephron, we notice many consistencies, along with a few inconsistencies. Sodium reabsorption in the proximal tubule is isotonic (i.e., sodium concentration is nearly constant along the proximal tubule); flow in the proximal tubule reduces to appro~mately 19% of its original value (slightly smaller than the accepted flow values which range from 20% to 35%); and solute concentrations and pressure reduce accurately. In general, the proximal tubule continues to behave accurately when incorporated into the rest of the model. Volume flow decreases and solute concentrations increase in the DHL, reflecting the passive extrusion of water into the medullary interstitium (counter-current multiplication). Sodium and chloride concentrations fall in AHL as expected, due to the active extrusion of salt, while volume flow remains constant, due to the impermeability of the AHL to water. Because the direction of flow in the AHL is toward the cortical interstitium, opposite that in other tubules, volume flow in the AHL is negative. In the DT and CD, volume flow reduces slightly and there is a net increase in the concentration of each solute. Pressure decreases accurately in
Solute and Water Transport
CONCENTRATION
81
OF UREA ALONG LENGTH OF NEPHRON
POSITION
Figure 6. PRESSURE
1
DHL
2
ALONG LENGTH OF NEPHRON
AHL 3 DT POSITION
4
CD
5
MI
6
Figure 7.
the DT and CD, as well. Flow in the MI is negative, reflecting that the direction of MI flow is toward the CI. Notice also the sodium concentration gradient which forms in the MI. This gradient is the driving force of the counter-current multiplier system, and its existence is essential in any accurate renal model. Furthermore, urea concentration increases almost continuously throughout the nephron. This behavior appropriately reflects one of the ultimate goals of the kidney, to create a waste highly concentrated in urea. Unfortunately, the physiological properties used to measure the accuracy of the results are only guidelines. With access to more exact flow values, the solution could be compared quantitatively. The results were also compared with Stephenson’s [17]. Again, because the two systems were normalized with respect to different flow values, his solutions could only represent directions to follow, rather than exact values that must be achieved. In general, the results compared favorably with Stephenson’s, with the following exceptions. There is a slight oscillatory behavior in the PT volume flow and in the sodium and potassium concentrations (we believe that these mild oscillations are numerical artifacts); Stephenson’s sodium concentration gradientin the MI
N. E. KOTTLER et al.
a2
is slightly more pronounced; we experienced a flat pressure profile and flat urea concentration gradients in the MI. The lack of a urea gradient, however, was also experienced by Stephenson [I?‘]. Therefore, this problem accurate
as a primary
improvement,
must be found.
initial
a more robust
Furthermore,
guess for this problem,
detailed
optimization
routine,
capable
of solving
renal flow data, which would present
must also be obtained.
Other
extensions
a more
to this research
include solving the parameter identification problem, assuming first constant, then nonconstant parameter values along each tubule; and implementing control to mimic the effects of both aldosterone and ADH. Hopefully, research in this area will be continued. what will someday ments
be a complete
this goal will be achieved
and accurate
At this time, we can only catch a glimpse
renal model.
However,
with continued
and will open the door to a more complete
of
improve-
understanding
of
renal physiology.
REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.
16. 17. 18. 19. 20. 21. 22.
J.J. Jacquez, B. Carnahan and P. Abbrecht, A model of the renal cortex and medulla, Mathematical Biosciences, An International Journal 1, 227-261 (1967). E. Koushanpour, R.R. Tarica and W.F. Stevens, Mathematical simulation of normal nephron function in rat and man, Journal of Theoretical Biology 31, 177-214 (1971). J.L. Stephenson, Transient behavior of the single loop cycling countercurrent multiplier, Bulletin of Mathematical Biology 35, 183-195 (1973). J.L. Stephenson, Concentrating engines and the kidney I. Central core model of the renal medulla, Biophysical Journal 13, 512-545 (1973). J.J. Jacquez and E. Daniels, Solute concentration in the kidney-I. A model of the renal medulla and its limit cases, Mathematical Biosciences, An International Journal 32, 307-335 (1976). D. Foster, J.J. Jacquez and E. Daniels, Solute concentration in the kidney-II. Input-output studies on a central core model, Mathematical Biosciences an International JountaZ32, 337-360 (1976). G.L. Barrett and J.S. Packer, Dynamic simulation of the renal medulla, Medical and Biological Engineering and Computing 21, 324-332 (983). A.S. Wexler, R.E. Kalaba and D.J. Marsh, Three-dimensional anatomy and renal concentrating mechanism. I. Modeling results, American Journal of Physiology 260, F368-F383 (1991). A.S. Wexler, R.E. Kalaba and D.J. Marsh, Three-dimensional anatomy and renal concentrating mechanism. II. Sensitivity results, American Jou7714l of Physiology 260, F384-F394 (1991). F.C. Hoppensteaclt and C.S. Peskin, Mathematics in Medicine and Life Sciences, Springer-Verlag, New York, (1992). J.L. Stephenson, Analysis of the transient behavior of kidney models, Bulletin of Mathematical Biology 40, 211-221 (1978). J.L. Stephenson, J.B. Garner and KS. Crump, Transient behavior of the single loop solute cycling model of the renal medulla, Bulletin of Mathematical Biology 40, 274-300 (1978). R. Mejia and J.L. Stephenson, Numerical solution of multinephron kidney equations, Jownal of Computational Physics 32, 235-246 (1979). J.B. Garner, R.B. Kellogg and J.L. Stephenson, Mathematical analysis of a model for the renal concentrating mechanism, Mathematical Biosciences, An International Journal 65, 125-150 (1983). R. Mejia and J.L. Stephenson, Solution of a multinephron, multisolute model of the mammalian kidney by Newton and continuation methods, Mathematical Biosciences, An International Journal 68, 279-298 (1984). J.L. Stephenson, Y. Zhang, A. Eftekhari and R. Tewarson, Electrolyte transport in a central core model of the renal medulla, American Journal of Physiology 253, F982-F997 (1987). J.L. Stephenson, Y. Zhang and R. Tewarson, Electrolyte, urea, and water transport in a two-nephron central core model of the renal medulla, American Journal of Physiology 257, F399-F413 (1989). D.E. Wessell, Parameter estimation for the proximal tubule in a single nephron steady state model of the kidney, Master’s Thesis, North Carolina State University, Raleigh, NC, (1993). IMSL, User’s Manual, IMSL Math/Libmry, FORTRAN Subroutines for Mathematical Applications, IMSL, Inc., Houston, TX, (1989). J.E. Dennis, Jr. and R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice-Hall, Englewood Cliffs, NJ, (1983). K.E. Brenan, S.L. Campbell and L.R. Petsold, Numerical Solution of InitiaG Value Problems in Dinerential-Algebraic Equations, SIAM Classics Series in Applied Mathematics, Philadelphia, PA, (1996). R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, New York, (1984).