(‘r,mp,,r.Biol. Med. PergamonPress1978. Vol 8,pp.5744. Printed mGreat Britain.

AN EFFICIENT NUMERICAL METHOD FOR SOLVING THE DIFFERENTIAL EQUATION OF RENAL COUNTERFLOW SYSTEMS* PARVE FARAHZAD~ Department

and REGINALD P. TEWARSON

of Applied Mathematics and Statistics, State University New York, Stony Brook, NY 11790. U.S.A.

(Received 15 December 1976; irl revisedform

of

I I April 1977)

Abstract-It is well known that the crucial step in the mathematical modeling of the kidney involves the determination of the roots of large systems of sparse non-linear algebraic equations. We utilize the information about the flow directions and the physiological connectivity of the various tubes in the kidney to develop a sparse matrix version of Newton’s method for the solution of these equations. The storage requirement for our method is independent of the number of tubes in the model. and its rate of convergence is that of Newton’s method. Non-linear ing storage

equations Minimizing

Sparse matrices computer time

Newton’s method Kidney Counter current multiplier

modeling MinimizMammalian kidney

INTRODUCTION

The mathematical models of renal counterflow systems contribute towards a better understanding of the renal function. These models involve systems of non-linear differential equations which, in general, can only be solved numerically [l, 23. The most suitable methods for handling such equations are known to be the finite difference schemes, which lead to non-linear systems of algebraic equations [3]. For a given model, the size of such a system is proportional to the number of discretization steps, and to the number of tubes and solutes in the model. In general, to obtain accurate solutions, small chop sizes (discretization steps) are needed, and this leads to a large number of equations. When solving such large systems, many practical difficulties are encountered. Attempts to overcome these difficulties led to the development of several methods C&6]. In this paper we will describe a new method which is a significant improvement in terms of storage and running time over the ones given in [46]. The motivation for developing the present method came from the need to handle comprehensive, multinephron models of the concentrating mechanism of the mammalian kidney without the use of auxiliary slow storage (discs or tapes). It has been pointed out recently that the functional unit of the kidney is not a single nephron but a nephrovascular unit consisting of a large number of tightly coupled nephrons [7]. The mathematical models describing such concentrating mechanisms involve a large number of interacting tubes which lead to non-linear algebraic equations of very large sizes. Fortunately, these equations are very sparse (there are very few unknowns in each equation), and therefore, sparse matrix techniques [8] can be developed for the efficient solution of these equations. In this paper we give a sparse matrix version of Newton’s method for solving these equations. The storage requirement of this method is independent of the number of tubes. Furthermore, its overall solution time is much less than for other methods. To illustrate our method, we describe a six tube vasa recta medullary model of the kidney in the next section. In the third Section we show how the physiological connectivity of the tubes in the kidney can be utilized to obtain a certain ordering of the equations and the unknowns. This ordering is used in the algorithm described in the fourth Section. *This research was supported in part by the National Institute of Arthritis, Metabolism and Diseases, NIH, Bethesda. Maryland. Grant No. RMOI AM 17593. t Present address: Maths Department, University of Vermont, Burlington, VT 05401. U.S.A. 57

Digestive

PARVIZ FARAHZAD and REGINALD P. TEWARSON

58

CORTEX

MEDULLA

Fig.

,

counter-flow system. (b) The axial 1. (a) A six tube vasa recta model of the medullary flows (-+) and the transmural fluxes w in the six tube vasa recta model.

In the final section, we give some computational with those published earlier. THE

results comparing the present method

MATirEMATICAL

MODEL

In the six tubes vasa recta model of the medullary counter-flow system (Fig. I), fluid enters tube 1 (descending Henle’s limb-DHL), flows through tube 2 (the ascending Henle’s lim&AHL), tube 3 (the distal nephron--l>N), and then to the tube 4 (the collecting duct-CD) where it emerges as concentrated urine. Blood enters tube 5 (the descending vasa recta-DVR) and leaves from tube 6 (the ascending vasa recta--AI/R), (see Fig. 1). We assume that tube 3 (DN) exchanges only with the Cortical Interstitium (C1ka bath of constant concentrations, while the other tubes exchange solute and water with tube 6 (AT/R). This model can be described by the following differential equations: kFitx)

&tci/cFi)

+

Jio(ci,

+

Jik(Ci,

C6r xl = O,

c6,

x,

=

O,

(2.la) (2.k)

where the subscripts i, k and v refer to the ith tube, the k’h solute and the volume. F is the axial volume flow, J denotes the outward transmural flux, C stands for the concentration, and X, which varies from 0 to 1, is the distance measured from the cortico-medullary junction down the medulla. For tube 3 (the DN), x is measured along the tube such that s = 0 and x = 1 correspond, respectively, to its junctions with tube 2 (AHL) and tube 4 (CD). Ci is the vector of concentrations whose elements

Method for solving equation of renal counteflow

systems

59

are Cik where k varies over all solutes. The J’s are known functions of their arguments and are assumed to satisfy:

Js,~= - ;

J,,,

(2.2)

where /I varies over all solutes and volume flows. The entering flows and concentrations in tubes 1 and 5, and the constant values of the concentration in CZ are given. In addition, the volume flows and the concentrations are assumed to match at the tube junctions. To obtain a finite difference analogue of (2.1), we first assume that each tube is divided into N equal parts of size h(h = l/N) and then perform all integrations by the trapezoidal rule. If we use a non-uniform chop size or other quadrature formulas for the integrations, then our sparse matrix method is still applicable with suitable modifications. The discrete forms of (2.la) and (2.lb) can be written as:

am -

Fi(j

-

1) + i(J,(j)

+ Jiu(j - 1)) = 0,

(2.3a)

and C,(,j)Fi(j) - Cik(j - l)Fi(j - 1) + i(Jia(j) +

Jdj - 1))= 0,

(2.3b)

where j varies from 1 to N. The difference equations (2.3) constitute a system of non-linear algebraic equations in terms of Cik(j) and F,(j). Because equations (2.3a) give easy access to Fio’yS as functions of the concentrations Ciko’ys,it has become a practice in solving such systems to delete F,(i)‘s from the equations, and to work with a reduced system of equations [2-71. Although the reduced system is somewhat smaller in size, the resulting equations are not as sparse as before. In the following sections we will give a sparse matrix method for solving the system of equations (2.3) simultaneously. The method requires approximately (P./S)’ storage locations, while the best known method for the reduced system [9] requires approximately T (N.S)’ storage locations, where T N and S are, respectively. the number of tubes, space chops and solutes in the model. STRUCTURE

OF THE

JACOBIAN

The so-called “occurrence matrix” of (2.3) after a suitable permutation, is shown in Fig. 2. The rows and columns of this matrix correspond to the equations and variables. Blanks indicate the absence of variables in the relevant equations. The particular structure shown in Fig. 2 can be obtained by a proper numbering of the equations and the unknowns. In this section we show how the physiological connectivity of the tubes in the kidney can be utilized to obtain such a numbering. From now on we will refer to FiO’)and CikO’)as the principal variables of the equations (2.3a) and (2.3b), respectively. The tube segments are numbered in the direction of flow so that the first and the last chops correspond to the top of tube 1 (AH-L), and tube 6 (AI/R), respectively (see Fig. 3). With regard to the above enumeration, the numbering of the equations is such that equations r(S + 1) + 1, r(S + 1) + 2,. . . , and (r -I- 1) (S + 1) are the ones obtained from the discretization of (2.1) at the rth chop. The unknowns are numbered in such a way that the pth unknown corresponds to the principal variable of the plh equation. In the remainder of this paper, the non-linear system (2.3) will be denoted by 4(Y) = 0,

(3.1)

where 4 and y are vectors whose components are the equations and variables of (2.3), respectively. We assume that the numbering is done in a manner which is consistent

60

PARVIZFARAHZADand REGINALDP. TEWARSON

Fig. 2. Structure of the “occurrence matrix”

for the six tube

model, after suitable permutation.

with that given in the previous paragraph. Thus the “occurrence matrix” of (3.1) is in the form of Fig. 2. Non-linear algebraic equations are usually solved by an iterative scheme. Among these algorithms, Newton’s method was found most suited to our problem. The use of other methods, e.g. non-linear Gauss Seidel, the Jacobi, and some storage saving projection methods [lo], leads to very slow convergence rates. For a description of these and other non-linear iterative methods, the reader is referred to [l 11. The Newton method applied to (3.1) is equivalent to solving the following sequence of linear equations: 4’(Y”)ZP = 4(YP)?

(3.2)

where Yp is the pth approximation to the solution of (3.1); +(yp) and #(yp) are, respectively, the non-linear system and its associated Jacobian evaluated at yp; and the solution Zp of (3.2) is the correction vector to yp. The next approximation yp+’ is given by: 1VP+1

=

YP_ZP.

(3.3)

Examination of equations (2.3) shows that if .Ii, and Jik are smooth functions of their arguments, then for all small values of h, all the diagonal elements of 4’(Y) should turn out to be excellent choices for pivoting. This information can be utilized to reduce the number of logical operations needed in solving (3.2). For example, in the Gaussian elimination method, the pivots can be taken on the diagonal. We recall that Fig. 2 shows the structure of the occurrence matrix for equation (3.1) and that the blanks in the matrix indicate that the associated equations do not contain the relevant variables. This implies that the corresponding elements of &(Y”) are zero. The information about the location of zeroes and the domination of the diagonal ele-

Fig. 3. Numbering of the tube-chops in the direction

of flow.

Method

for solving

equation

of renal counterflow

Fig. 4. Shaded blocks corresponding to the elements of A which typical Gaussian elimination step, when pivots are taken

hl

systems

will be modified on the diagonal.

during

a

ments can be utilized to give efficient sparse matrix methods for solving (3.2). In the next section we will describe one such method. A SPARSE

MATRIX

METHOD

We now describe how (3.2) can be solved by Gaussian elimination. we drop the iteration variable p and write (3.2) as AZ = b,

For convenience (4.1)

where A = @(y”) and b = 4(yp). The computations involved in solving (4.1) can be arranged in such a manner that the storage requirement at each elimination step is very small. In the previous section we mentioned that the diagonal elements of A are good candidates for pivoting. We will use this pivoting strategy in the elimination process. The cases in which these pivots turn out to be small will be discussed later. In view of the above-mentioned fact and Fig. 2, it is evident that if the pivots are taken on the diagonal starting from the top left hand corner of A, then Gaussian elimination does not create any additional non-zero elements. It is also clear from Fig. 2 that at each elimination step, only a small portion of A will be modified. We have illustrated this for a typical elimination step by the shaded regions in Fig. 4. At each stage, we only need to store those elements which will be modified. This suggests an efficient sparse matrix method for solving (4.1). Since our main objective for finding the solution of (4.1) is to update yp, the following algorithm produces the updated yP+ ’ directly.

Algorithm

Let S, T and N, respectively, denote the number of solutes, tubes and space chops. Choose W and d as two working areas of sizes [S(N + 2) + l] x [S(N + 3) + 21, and [S(N + 2) + l] x 1, respectively. Let W and d be partitioned as shown in Fig. 5, and do the following steps. 1. Generate the bottom right corner block of size S(N + 1) x S(N + I), and the top right corner block of size 2(S + 1) x S(N + 1) of A. Store the first block in W,, and the second one in WI, and Wz2, respectively. 2. Generate the first 2(S + 1) and the last S(N + 1) elements of b. Store the first block in d1 and d2 and the second block in d,. 3. Let the subiteration variable I = 1. 4. Generate the blocks of A corresponding to columns (A - l)(S + 1) + 1 to A(S + 1) androws@. - l)(S + 1) + lto(1 + l)(S + l),andN(S + l)(T - 1) + 1 to(S + l)N(T 1) + NS. Store the first block in WI 1 and Wllr and the second block in W, 1. 5. Use the Forward course of the Gaussian Elimination, with pivots on the diagonal.

62

PARVIZ FARAHZAD and REGINALDP. TEWARSQN

WI2 w22

w32

Fig. 5. Partitioning of the working area for the sparse

matrix method.

to reduce the first (S + 1) columns of lV, to an upper triangular form. Do the same operations on the right hand side “d.” 6. If 1 2 N(T - l), go to 9; otherwise continue. 7. Replace W,, and d, by W,, and d,, respectively. 8. Generate the block, which refers to (A + 2) (S + 1) + 1 to (A + 3) (S + 1) + 1 rows, and the last S(N + 1) columns of A, and store it in W,,. Compute the block corresponding to (A + 2)(S + 1) + 1 to (A + 3) (S + 1) elements of b; store it in d,. Increment ,? by 1 and go to 4. 9. Solve the system of linear equations, W,, ZA= d,,

(4.2)

by some inversion technique, and store the solution za, and the inverse W3;’ in d, and W32, respectively. Note that saving of this inverse is only necessary when the fixed Newton iterations are to be used. 10. Using (3.3) and the correction vector zA, which is available in d,, update the last S(N + 1) elements of y. Use this updated y in the following calculations, and for v = 1,2,.... , (S - 1) (T - 1) - 1, and (S - 1) (T - I), do the following steps. 11. Generate the diagonal block, which corresponds to (v - 1) (S + 1) + 1 to v(S + 1) rows and columns of A, and store it in WI,. Compute the block referring to (v - 1) (S + 1) + 1 to v(S + 1) elements of b and store it in dl. 12. Solve the system of linear equation W11&=

d1,

(4.3)

by Gaussian elimination, with pivots on the diagonal. 13. By using (3.3) and the correction vector z,,, update the (v - 1) (S + 1) + 1 to v(S + 1) elements of y. Use this updated y in the following calculations. 14. If v 2 (S - l)(T - l), stop; otherwise, increment v by 1 and go to 11. This completes one iteration of Newton’s method. For fixed Newton iterations, start with Step 2. Restrict all operations on the W matrix to w1, W,, and W13, leaving the remainder of the matrix intact, and replace Step 9 by: 9’. Compute z1 by: zA = W,, d,.

(4.4)

This completes one iteration of the fixed Newton’s method. Note that the usual fixed Newton method requires the complete storage of [+‘I-‘, the inverse of the Jacobian, but in our case we have reduced the storage requirement substantially. In situations where a particular diagonal element Uii of A is too small to be used for pivoting, the variable Zi in (4.1) is set equal to zero, and the reduced system is solved. This corresponds to a modified version of Newton’s method, where at each iteration only a subset of the elements of yp are being updated. Our experience has

Method for solving equation of renal counterflow systems

63

confirmed the suitability of such diagonals for pivoting. In all cases that we considered. no difficulties were encountered from such a strategy. COMPUTATIONAL

RESULTS

For numerical experiments, we selected a two solute kidney model of the form described in the second Section. The non-zero transmural fluxes were chosen as follows (s and u denote, respectively, salt and urea): J1, = 10(% + c6u - cl, - clu). Jzs = 1.8

0 < .X