DOI: 10.2478/amcs-2014-0062
A MATHEMATICAL MODEL FOR FLUID–GLUCOSE–ALBUMIN TRANSPORT IN PERITONEAL DIALYSIS
ROMANCHERNIHA∗,∗∗∗, JOANNASTACHOWSKA-PI ˛ETKA∗∗, JACEKWANIEWSKI∗∗
∗Institute of Mathematics
National Academy of Sciences of Ukraine, Tereshchenkivs’ka Street 3, 01601 Kyiv, Ukraine e-mail:cherniha@gmail.com
∗∗Nał˛ecz Institute of Biocybernetics and Biomedical Engineering Polish Academy of Sciences, Ks. Trojdena 4, 02-109 Warsaw, Poland
e-mail:{jstachowska,jacekwan}@ibib.waw.pl
∗∗∗School of Mathematical Sciences
University of Nottingham, University Park, Nottingham, NG7 2RD, UK e-mail:Roman.Cherniha@nottingham.ac.uk
A mathematical model for fluid and solute transport in peritoneal dialysis is constructed. The model is based on a three- component nonlinear system of two-dimensional partial differential equations for fluid, glucose and albumin transport with the relevant boundary and initial conditions. Our aim is to model ultrafiltration of water combined with inflow of glucose to the tissue and removal of albumin from the body during dialysis, by finding the spatial distributions of glucose and albumin concentrations as well as hydrostatic pressure. The model is developed in one spatial dimension approximation, and a governing equation for each of the variables is derived from physical principles. Under some assumptions the model can be simplified to obtain exact formulae for spatially non-uniform steady-state solutions. As a result, the exact formulae for fluid fluxes from blood to the tissue and across the tissue are constructed, together with two linear autonomous ODEs for glucose and albumin concentrations in the tissue. The obtained analytical results are checked for their applicability for the description of fluid-glucose-albumin transport during peritoneal dialysis.
Keywords: fluid transport, transport in peritoneal dialysis, nonlinear partial differential equation, ordinary differential equ- ation, steady-state solution.
1. Introduction
Peritoneal dialysis is a life saving treatment for chronic patients with end-stage renal disease (Gokal and Nolph 1994). The peritoneal cavity, an empty space that separates bowels, abdominal muscles and other organs in the abdominal cavity, is applied as a container for the dialysis fluid, which is infused there through a permanent catheter and left in the cavity for a few hours.
During this time small metabolites (urea, creatinine) and other uremic toxins diffuse from blood that perfuses the tissue layers close to the peritoneal cavity to the dialysis fluid, and finally are removed together with the drained fluid. The treatment cycle (infusion, dwell, drainage) is repeated several times every day. The peritoneal transport occurs between the dialysis fluid in the peritoneal cavity
and blood passing down the capillaries in the tissue surrounding the peritoneal cavity. The capillaries are distributed within the tissue at different distances from the tissue surface in contact with the dialysis fluid. The solutes, which are transported between blood and the dialysis fluid, have to cross two transport barriers: the capillary wall and a tissue layer (Flessner 2006).
Typically, many solutes are transported from blood to dialysate, but some solutes such as, for example, an osmotic agent (typically glucose), which is present in high concentration in dialysis fluid, are transported in the opposite direction, i.e., to the blood. This kind of transport system appears also in other medical treatments, such as local delivery of anticancer medications, and in some experimental and natural physiological phenomena (see below). Typically a distributed approach is applied to take
into account spatial properties of these systems. The first applications of the distributed model were limited to the diffusive transport of gases between blood and artificial gas pockets within the body (Piiper et al., 1962), between subcutaneous pockets and blood (Van Liew, 1968; Collins, 1981), and transport of heat and solutes between blood and the tissue (Perl, 1963; 1962).
Applications of the distributed approach to the modeling of diffusive transport of small solutes include the description of transport from the cerebrospinal fluid to the brain (Patlak 1975), delivery of drugs to the human bladder during intravesical chemotherapy, and drug delivery from the skin surface to the dermis in normal and cancer tissue (Gupta et al., 1995; Wientjes et al., 1993; 1991). Finally, the distributed approach was also proposed for a theoretical description of fluid and solute transport in solid tumors (Baxter and Jain, 1989; 1990;
1991). The mathematical descriptions of these systems were obtained using partial differential equations based on the simplification that capillaries are homogeneously distributed within the tissue. Experimental evidence confirmed the good applicability of such models (see, for example, the papers of Waniewski et al. (1996a; 1996b;
2007), Smit et al. (2004a), Flessner (2006), Parikova et al. (2006), Stachowska-Pietka et al. (2012), Guest et al.
(2012) and the references therein).
An important objective of peritoneal dialysis is to remove excess water from the patient (Gokal and Nolph, 1994). Typical values of water ultrafiltration measured during peritoneal dialysis are 10–20 mL/min (Heimbürger et al., 1992; Waniewski et al., 1996a, 1996b; Smit et al., 2004a; 2004b). This is achieved by inducing osmotic pressure in the dialysis fluid by adding a solute (called the osmotic agent) in high concentration. The most frequently used osmotic agent is glucose. This medical application of high osmotic pressure is unique for peritoneal dialysis.
The flow of water from blood across the tissue to the dialysis fluid in the peritoneal cavity carries solutes of different size, including large proteins, and adds a convective component to their diffusive transport.
The mathematical description of fluid and solute transport between blood and dialysis fluid in the peritoneal cavity has not been fully formulated yet, in spite of the well-known basic physical laws for such transport. The complexity of peritoneal fluid transport modelling comes mainly from the fact that, whereas diffusive transport of small solutes is linear, the process of water removal during peritoneal dialysis by osmosis is nonlinear. A first formulation of the general distributed model for combined solute and fluid transport was proposed by Flessner et al. (1984) and applied later for the description of the peritoneal transport of small molecules (Flessner et al., 1985).
The next attempt to model fluid and solute transport did not result in a satisfactory description. It was assumed
in that model that the mesothelium is a very efficient osmotic barrier for glucose with the same transport characteristics as the endothelium (Seams et al., 1990).
The assumption resulted in negative interstitial hydrostatic pressures during osmotically driven ultrafiltration from blood to the peritoneal cavity during peritoneal dialysis (Seams et al., 1990). This contradicted the experimental evidence on positive interstitial hydrostatic pressure during the ultrafiltration period of peritoneal dialysis (Flessner 1994). Moreover, the mesothelium, being a highly permeable layer, cannot provide enough resistance to small solute transport to be an osmotic barrier for solutes such as glucose (Flessner, 1994; 2006; 2009;
Czyzewska et al., 2000).
Recent mathematical, theoretical and numerical studies have introduced new concepts on peritoneal transport and yielded a better description of particular processes such as pure water transport, a combined osmotic fluid flow and small solute transport, and transport of water and proteins (Flessner, 2001; Cherniha and Waniewski, 2005; Stachowska-Pietka et al., 2006;
2007; Cherniha et al., 2007; Waniewski et al., 2007;
2009). A recent study by Stachowska-Pietka et al. (2012) addresses again the combined transport of fluid (water) and several small solutes.
However, the problem of a combined description of osmotic ultrafiltration to the peritoneal cavity, absorption of the osmotic agent from the peritoneal cavity and leak of macromolecules (e.g., albumin) from blood to the peritoneal cavity has not been addressed yet. Therefore, we present here an extended model for these phenomena and investigate its mathematical structure. In particular, we investigate fundamental questions concerning the role of various transport components, such as osmotic and oncotic gradients and the hydrostatic pressure gradient.
It should be stressed that the oncotic gradient leading to the leak of macromolecules from blood to the peritoneal cavity has the opposite sign to the osmotic gradient; hence, their combination may lead to new effects, which cannot arise in the case of the simplified models mentioned above. The model addresses also the disputed problem of the values of the Staverman reflection coefficients for glucose and albumin in the tissue, which cannot be directly measured and need to be indirectly derived from clinical data (Katchalsky and Curran, 1965;
Waniewski, 2013).
The paper is organized as follows. In Section 2, a mathematical model of glucose and albumin transport in peritoneal dialysis is constructed. In Section 3, non-uniform steady-state solutions of the model are constructed and their properties are investigated.
Moreover, these solutions are tested using real parameters that represent clinical treatments of peritoneal dialysis.
The results are compared with those derived by numerical simulations for simplified models (Cherniha et al.,
2007; Waniewski et al., 2007). Finally, we present some conclusions and discussion in the last section.
2. Mathematical model
Here we present a new model of fluid and solute transport in peritoneal dialysis. The model is developed in one spatial dimension with x = 0 representing the boundary of the peritoneal cavity and x = L representing the end of the tissue surrounding the peritoneal cavity (see the work of Stachowska-Pietka et al. (2012) for a discussion of the assumptions involved in this approach).
The mathematical description of transport processes within the tissue consists of the local balance of the fluid volume and solute mass. For incompressible fluid, the change of volume may occur only due to elasticity of the tissue. The fractional fluid void volume, i.e., the volume occupied by the fluid in the interstitium (the rest of the tissue being cells and macromolecules forming the solid structure of the interstitium) expressed per unit volume of the whole tissue, is denoted by ν(t, x), and its time evolution is described as
∂ν
∂t = −∂jU
∂x + qU − ql, (1) where jU(t, x) is the volumetric fluid flux across the tissue (ultrafiltration), qU(t, x) is the density of the volumetric fluid flux from blood capillaries to the tissue, and qlis the density of the volumetric fluid flux from the tissue to the lymphatic vessels (hereafter we assume that it is a known positive constant: nevertheless, it can also be a function of hydrostatic pressure (Stachowska-Pietka et al., 2006;
2012).
Similarly to many distributed models, our model involves the spreading of the source within the whole tissue as an approximation to the discrete structure of blood and lymphatic capillaries.
The independent variables are time t and distance x within the tissue from the tissue surface in contact with the dialysis fluid in the direction perpendicular to this surface (flat geometry of the tissue is here assumed with a finite width, see below). The solutes, glucose and albumin, are distributed only within the interstitial fluid (or part of it, see below), and their concentrations in this fluid are denoted by CG(t, x) and CA(t, x), respectively.
The equation that describes the local changes of glucose amount in the tissue, νCG, is
∂(νCG)
∂t = −∂jG
∂x + qG, (2) where jG(t, x) is the glucose flux through the tissue, and qG(t, x) is the density of the glucose flux from blood.
The cellular uptake of absorbed glucose is not taken into account in Eqn. (2) because this process leads to a small correction in bulk absorption of glucose to the capillaries.
So, we neglect these intracellular changes that have been observed experimentally (Zakaria et al., 2000).
Similarly, the equation that describes the local changes of albumin amount in the tissue, ανCA, is
∂(ανCA)
∂t = −∂jA
∂x + qA, (3) where jA(t, x) is the albumin flux through the tissue, qA(t, x) is the density of the albumin flux from blood. The coefficient α < 1 takes into account that only a part of the fractional fluid void volume ν that is available for fluid is accessible for albumin because of its large molecular size (Flessner, 2001; Stachowska-Pietka et al., 2007). In other words, the inclusion of the term αν in (3) implies that CA(t, x) is the concentration of albumin in that part of the interstitium across which the albumin molecules can pass.
In the general case, Eqn. (3) involves a new fluid void volume function νA(t, x), which depends on the hydrostatic pressure similarly to the function ν (see below) and satisfies the inequality νA < ν. Hereafter we set νA= αν for simplicity.
The flows of fluid and solutes through the tissue are described according to linear non-equilibrium thermodynamics. Osmotic pressure of glucose and oncotic pressure of albumin are described by van’t Hoff law, i.e., assuming that the corresponding pressures are proportional to the relevant concentrations.
The fluid flux across the tissue is generated by hydrostatic, osmotic, and oncotic (i.e., osmotic pressure of large proteins) pressure gradients according to the extended Darcy law (Katchalsky and Curran 1965):
jU = −νK∂P
∂x + σT GνKRT∂CG
∂x + σT AνKRT∂CA
∂x , (4)
where K is the hydraulic conductivity of the tissue that is assumed constant for simplicity (K may also depend on the pressure P ), R is the gas constant, T is absolute temperature, and σT G and σT A are the Staverman reflection coefficients for glucose and albumin in the tissue, respectively. The Staverman reflection coefficient σ is a thermodynamic parameter and describes the effectiveness of osmotic pressure in selectively permeable membrane: if σ = 0, then no osmotic pressure can be induced by this solute across the membrane, and if σ = 1, the maximal theoretically possible osmotic effect is induced (ideal semi-permeable membrane).
The fluid flows along the hydrostatic pressure gradient (i.e., from higher to lower pressure) and up the concentration gradient of the osmotic agent (i.e., from solution of lower concentration to solution of higher concentration, osmosis), and therefore there is a difference in signs at the gradients in Eqn. (4). Intermediate values of σ represent non-ideal semipermeable membranes.
Katchalsky and Curran (1965) thoroughly address the problem of Staverman reflection coefficients.
The density of the fluid flux from blood to the tissue is generated, according to the Starling law, by the hydrostatic, osmotic and oncotic pressure differences between blood and the tissue:
qU = Lpa(PB− P ) − LpaσGRT (CGB− CG)
−LpaσART (CAB− CA), (5) where P (t, x) is the hydrostatic pressure, Lpa is the hydraulic conductance of the capillary wall, PB is the hydrostatic pressure of blood, CGB and CAB are the glucose and albumin concentrations in blood, and σG
and σAare Staverman reflection coefficients for glucose and albumin in the capillary wall, respectively. Contrary to other parameters, there is an unsolved problem of the values of σG and σA. In particular, the values of σG were found low (about 0.005 − 0.03) in many experiments in contrast to some newer experimental data that suggest the values close to 0.5 (see the discussion of this controversy in the works of Waniewski et al.
(2009) and Stachowska-Pietka et al. (2012)). We also assume that blood concentrations of glucose and albumin are constant according to clinical and experimental data that demonstrate only negligible variation of these concentrations during peritoneal dwell of dialysis fluid (Heimburger et al. 1992). This observation is related to the quasi-continuous mode of continuous ambulatory peritoneal dialysis with fluid exchanges every few hours and was applied in most previous theoretical and numerical studies on peritoneal dialysis.
The glucose flux across the tissue is composed of a diffusive component (proportional to the glucose concentration gradient) and a convective component (proportional to glucose concentration and the fluid flux)
jG = −νDG∂CG
∂x + ST GCGjU, (6) where DGis the diffusivity of glucose in tissue, ST Gis the sieving coefficients of glucose in the tissue. According to non-equilibrium thermodynamics, ST G = 1 − σT Gfor a homogenous membrane (Katchalsky and Curran 1965).
The density of the glucose flux between blood and tissue describes the number of moles of glucose per unit total volume of tissue per unit time that move between blood and the tissue. It has a diffusive component (proportional to the difference of glucose concentration in blood, CGB, and glucose concentration in the tissue, CG), a convective component (proportional to the density of the fluid flow from blood to the tissue, qU) and a component that represents lymphatic absorption of solutes (proportional to the density of the volumetric lymph flux, ql):
qG= pGa(CGB− CG) + SGqUCG− qlCG. (7)
where pGa is the diffusive permeability of the capillary wall for glucose.
In a similar way, the albumin flux across the tissue, jA(t, x), and the density of the albumin flux to the tissue, qA(t, x), can be described as
jA= −ανDA∂CA
∂x + ST ACAjU, (8) qA= pAa(CAB− CA) + SAqUCA− qlCA, (9) where ST A = 1 − σT A is the sieving coefficient of albumin in the tissue, SA = 1 − σA is the sieving coefficient of glucose and albumin in the capillary wall, DA is the diffusivity of albumin in the tissue, and pAa is the diffusive permeability of the capillary wall for albumin. Typical values of the model parameters are listed in Table 1.
Equations (1)–(3), together with Eqns. (4)–(9) for flows, form a system of three nonlinear partial differential equations with four variables: ν, P, CA, and CG. Therefore, an additional, constitutive, equation is necessary, and this is the equation describing how fractional fluid void volume, ν, depends on interstitial pressure, P . This dependence can be established using data from experimental studies (Stachowska-Pietka et al., 2006). It turns out that
ν = F (P ), (10)
where F is a monotonically non-decreasing bounded function with the limits F → νmin if P → Pmin
and F → νmax if P → Pmax (particularly, one may take Pmin = −∞, Pmax = ∞). Here νmin < 1 and νmax < 1 are empirically measured constants. An example of the function F based on experimental data taken from the work of Zakaria et al. (1999) was presented by Stachowska-Pietka et al. (2006) and Cherniha et al.
(2007).
Boundary conditions for a tissue layer of width L impermeable at x = L and in contact with the dialysis fluid at x = 0 are
x = 0 : P = PD, CG= CGD, CA= CAD, (11) ex = L : ∂P
∂x = 0, ∂CG
∂x = 0, ∂CA
∂x = 0. (12) Generally speaking, intraperitoneal pressure PD, glucose CGD and albumin CAD concentrations in the peritoneal cavity may depend on time. However, experimental data and theoretical studies suggest that they change at a low rate compared to the rate of transport processes in the tissue (Stachowska-Pietka et al., 2006; 2007; Waniewski, 2007). Therefore, we may assume that PD, CGD and CAD are constant for some time period and assess the steady-state solution for these particular boundary conditions that may be considered an approximated quasi
steady-state solution for the full model of peritoneal dialysis with time-dependent boundary conditions. This approximation was applied to previous models with variable boundary conditions for small solutes (as glucose) and water transport (but without proteins, as albumin) (see Dedrick, 1981; Flessner, et al., 1984; 1985;
Seames et al., 1990; Stachowska-Pietka et al., 2006; 2007;
Waniewski, 2001; 2002; Waniewski et al., 2009).
The initial conditions describe an equilibrium within the tissue without any contact with the dialysis fluid:
t = 0 : P = P∗, CG= CG∗, CA= CA∗, (13) where P∗, CG∗, and CA∗ are some non-negative values, which will be estimated below.
Note that Eqns. (1)–(10) can be united into three nonlinear partial differential equations (PDEs) for hydrostatic pressure P (t, x), glucose concentration CG(t, x) and albumin concentration CA(t, x). Thus, these three PDEs, together with the boundary and initial conditions (11)–(13), form a nonlinear boundary-value problem. Possible values of the parameters arising in this problem are presented in Table 1 (see the relevant comments in Section 4).
The fluid flux jU(t, x) at x = 0 describes the net ultrafiltration flow, i.e., the exchange of fluid between the tissue and the peritoneal cavity across the peritoneal surface, and therefore directly the efficiency of water removal during peritoneal dialysis. The assessment of ultrafiltration flow is important from a practical point of view because low values of this flow in some patients indicate that some problems with osmotic fluid removal have occurred, which may finally result in a failure of the therapy (Parikova, 2006).
3. Steady-state solutions of the model and their applications
The time needed to approach the steady state is of the order of minutes for small solutes, as glucose, but increases for larger solutes, and for albumin may be expected to be a few hours (Waniewski, 2007; 2008;
Stachowska-Pietka et al., 2006). Thus, the transport of small solutes is close to the steady state even during a single exchange of the dialysis fluid. We should, however, take also into account that patients are in continuously repeated treatment with several hours of fluid dwell in the peritoneal cavity and 10–20 minutes of fluid exchange, and there are typically four exchanges every day (Stachowska-Pietka et al., 2006). Thus, the transport system for large molecules after many exchanges is also close to the steady state. Hence, solutions for the steady state of the system may be considered good approximations for real conditions in the tissue in this clinical setting.
First of all, we consider the special case of tissue in its physiological state without dialysis, and, therefore, no transport to the peritoneal cavity occurs. In this case the boundary conditions at x = 0 given by Eqn. (11) are replaced by zero Neumann conditions, and the steady-state solution can be easily found because it does not depend on x. In fact, by solving the algebraic equations
qU− ql= 0, qG= 0, qA= 0, (14) one easily obtains the spatially uniform steady-state solution
CG∗ = pGa
pGa + σGqlCGB, CA∗ = pAa
pAa + σAqlCAB, P∗= PB− ql
Lpa
− qlRT
σ2GCGB
pGa + σGql + σ2ACAB
pAa + σAql
.
(15)
In the case ql = 0, i.e., the zero flux from the tissue to the lymphatic vessels, the formulae (15) produce
CG∗ = CGB, CA∗ = CAB, P∗= PB, (16) otherwise
CG∗ < CGB, CA∗ < CAB, P∗< PB. (17) This uniform solution describes the system in equilibrium if no dialysis is performed, and therefore we may use the values P∗, CG∗, and CA∗ calculated above as the initial profile for simulation of the transport processes after the initiation of dialysis (see the formulae (13)).
To find spatially non-uniform steady-state solutions, we reduce Eqns. (1)–(3) to an equivalent form by introducing scaled non-dimensional independent and dependent variables (except for ν, which is a non-dimensional variable),
x∗= x
L, t∗=KPDt
L2 , (18)
p(t∗, x∗) = P PD, u(t∗, x∗) = CG− CGB
CGD− GGB, w(t∗, x∗) = CA
CGD− GGB.
(19)
Here we consider only the case with CGD > CGB, and therefore variables u and w are well defined and positive, if the appropriate initial conditions, Eqns. (13), (16) and (17), are assumed. For alternative methods of scaling, see the works of Waniewski (2001; 2002).
Thus, taking into account Eqns. (4), (6), and (8), and after rather simple calculations, one obtains Eqns. (1)–(3) in the form (hereafter the upper index∗ is omitted)
∂ν
∂t = ∂
∂x
ν∂p
∂x
− t0σ1 ∂
∂x
ν∂u
∂x
− t0σ2 ∂
∂x
ν∂w
∂x
+ t0(qU − ql), (20)
∂
∂t
ν(u + u0)
= d1t0 ∂
∂x
ν∂u
∂x
+ ST G ∂
∂x
uν∂p
∂x
− ST Gt0σ1 ∂
∂x
uν∂u
∂x
− ST Gt0σ2 ∂
∂x
uν∂w
∂x
− ST Gu0 ∂
∂x
ν∂p
∂x
− ST Gu0t0σ1 ∂
∂x
ν∂u
∂x
+ ST Gu0t0σ2 ∂
∂x
ν∂w
∂x
+ t0SG(u + u0)qU− t0b1u − t0u0ql,
(21)
∂(ανw)
∂t
= d2t0 ∂
∂x
ν∂w
∂x
+ ST A ∂
∂x
wν∂p
∂x
− ST At0σ1 ∂
∂x
wν∂u
∂x
− ST At0σ2 ∂
∂x
wν∂w
∂x
+ t0SAwqU − t0b2w − t0w0(b2− ql),
(22)
where
qU = β
1
t0(p0− p) + σGσ1
σT G u +σAσ2
σT Aw∗
,
σ1= σT GKRTCGD− GGB
L2 , d1=DG
L2, σ2= σT AKRTCGD− GGB
L2 , d2= αDA
L2 , b1= pGa + ql, b2= pAa + ql, β = LpaL2
K , u0= CGB
CGD− GGB, w0= CAB
CGD− GGB, p0= PB
PD, t0= L2 KPD.
(23)
Since Eqns. (21)–(22) have essentially a different structure, we modify them using Eqn. (20) and the notation w∗ = w − w0; hence, the system of governing
equations consisting of Eqn. (20) and the equations
∂(νu)
∂t + σT Gu0∂ν
∂t
= d1t0 ∂
∂x
ν∂u
∂x
+ ST G ∂
∂x
uν∂p
∂x
− ST Gt0σ1 ∂
∂x
uν∂u
∂x
− ST Gt0σ2 ∂
∂x
uν∂w
∂x
+ t0(SGu + u0(SG− ST G))qU − t0b1u
− t0σT Gu0ql,
(24)
∂(ανw)
∂t − ST Aw0∂ν
∂t
= d2t0 ∂
∂x
ν∂w
∂x
+ ST A ∂
∂x
w∗ν∂p
∂x
− ST At0σ1 ∂
∂x
w∗ν∂u
∂x
− ST At0σ2 ∂
∂x
w∗ν∂w
∂x
+ t0(SAw − ST Aw0))qU− t0b2w∗− t0σT Aw0ql
(25) is obtained.
We want to find the steady-state solutions of Eqns. (20), (24)–(25) satisfying the boundary conditions (11)–(12). They take the form
x = 0 : p = 1, u = 1, w = CAD
CGD− GGB, (26) x = 1 : ∂p
∂x = 0, ∂u
∂x = 0, ∂w
∂x = 0. (27) for the non-dimensional variables.
Note that, to find the steady-state solutions, Eqns. (20), (24)–(25) can be reduced to the system of ordinary differential equations (ODEs)
1 t0
d dx
νdp dx
− σ1 d dx
νdu dx
− σ2 d dx
νdw dx
= ql− qU, (28)
d1 d dx
νdu dx
+ST G
t0
d dx
uνdp dx
− ST Gσ1 d dx
uνdu dx
− ST Gσ2 d dx
uνdw dx
+ (SGu + u0(SG− ST G))qU− b1u − σT Gu0ql= 0, (29) d2 d
dx
νdw dx
+ST A
t0
d dx
w∗νdp dx
− ST Aσ1 d dx
w∗νdu dx
− ST Aσ2 d dx
w∗νdw dx
+ (SAw − ST Aw0)qU− b2w∗− σT Aw0ql= 0.
(30)
The non-linear system of ODEs (28)–(30) is still very complex and cannot be integrated in the case of arbitrary coefficients. Thus, we look for correctly specified coefficients, for which this system can be simplified. It can be noted that the relations
SA= ST A, SG= ST G (31) lead to an essential (this means that automatically σG = σT G, σA = σT A) simplification of this system. This assumption is introduced for mathematical reasons only:
a specific symmetry of the equations allows much easier rigorous analysis. On the other hand, it is shown in the next section that even in this special case the solutions of the model are qualitatively/quantitatively similar to those obtained via other simplified models, which do not use this assumption.
So, using the assumption (31), expressions for qU
from (23) and jU from (4), rewritten in non-dimensional variables
jU = Lν
−1 t0
∂p
∂x+ σ1∂u
∂x + σ2∂w
∂x
, (32) we obtain the relation
jU = Kν LpaL
dqU
dx (33)
allowing jU to be found provided the function qU is known. Using the formulae (31)–(33), the nonlinear ODE system (28)–(30) can be simplified to the form
1 β
d dx
νdqU
dx
− qU+ ql= 0, (34)
d1 d dx
νdu
dx
−SG
β d dx
νudqU
dx
+ SGuqU − b1u − σGu0ql= 0,
(35)
d2 d dx
νdw
dx
−SA
β d dx
νwdqU
dx
+ SAwqU− b2(w − w0) − w0ql= 0.
(36)
The linear semi-coupled system of ODEs (33)–(34) can be extracted to find the functions qU and jU provided the function ν is known. However, ν depends on pressure, which is also an unknown function, and therefore we need to assume additional restrictions on the function F from the formula (10).
Let us consider the first case, in which we assume that F is a constant function. This assumption was applied in many studies, especially for the description of small solute transport (Dedrick, 1981; Flessner, 1984;
Waniewski, 2001; 2002). In this case,
ν(x) = νm, (37)
where νmis a positive constant. Substituting (37) into the system (33)–(34), its general solution can be found:
qU = C1e−λx+ C2eλx+ ql, (38)
jU =L
λ(−C1e−λx+ C2eλx), λ =
LpaL2
Kνm . (39) The arbitrary constants C1and C2can be specified using the boundary conditions (26)–(27) since the functions qU
and jU are expressed via p, u, w and its first-order derivatives (see the formulae (23) and (32)). Making rather simple calculations, we obtain
C1= (q0− ql) e2λ
1 + e2λ, C2= (q0− ql) 1
1 + e2λ, (40) where
q0= β
1
t0(p0− 1) + σ1+ σ2CAD− CAB CGD− GGB
. (41) Having the explicit formulae for qU and jU, Eqns. (35) and (36) can be reduced to two linear autonomous ODEs:
d1νmd2u dx2+SG
λ (C1e−λx−C2eλx)du
dx−κ1u = u01 (42) and
d2νmd2w dx2 +SA
λ (C1e−λx− C2eλx)dw
− κ2(w − wdx0) = w01 (43)
with unknown functions u(x) and w(x). Hereafter the notation
κ1= pGa + σGql, u01= σGu0ql,
κ2= pAa + σAql, w01= σAw0ql (44) is used. Note the similarities in the structure of Eqns. (42) and (43) However, to the best of our knowledge, the general solutions of the ODEs (42) and (43) are unknown.
On the other hand, since the unknown functions u(x) and w(x) should satisfy the boundary conditions (26)–(27), the corresponding linear problems can be numerically solved using, e.g., the Maple program package. Finally, using two expressions for qU from (23) and (38), we obtain the function
p(x) = p0+ t0σ1u + t0σ2(w − w0)
−t0
β
C1e−λx+ C2eλx+ ql
. (45) In the next section, the numerical non-uniform steady-state solutions based on realistic parameter values arising in the formulae derived above will be presented for this case, i.e., with the restrictions (31) and (37).
Let us now consider the second type of restriction on the function ν. Instead of the rather restrictive assumption (37), we examine the case when the function ν is non-constant and satisfies the general conditions described after the formula (10). According to the experimental data, the hydrostatic pressure during peritoneal dialysis is a decreasing function with respect to the distance x from the peritoneal cavity (Flessner, 1994;
Zakaria et al. 1999; 2000). Hence, function F (p(x)) is decreasing (with respect to x) provided p(x) is a spatially non-uniform steady-state solution. The simplest case of such a pattern occurs when ν is the linear, monotonically decreasing function of x:
ν(p(x)) ≡ ν(x) = νmax− (νmax− νmin)x, (46) where x ∈ [0, 1]. Substituting (46) into (34), we obtain the linear ODE
(νmax− (νmax− νmin)x)d2qU
dx2
− (νmax− νmin)dqU
dx − β(qU− ql) = 0. (47) It can be shown by substitution
y2= 4δ∗(ν∗− x), qV = qU − ql, ν∗= νmax
νmax− νmin > 1,
δ∗= β
(νmax− νmin)> 0
(48)
that the linear ODE (Eqn. (47)) reduces to the modified Bessel equation of the zero order (see, e.g., Polyanin and Zaitsev, 2003)
y2d2qV
dy2 + ydqV
dy − y2qV = 0. (49) The general solution of Eqn. (49) is well known.
Hence, using the formulae (48), we obtain the solution of Eqn. (47):
qU = C1I0(2
δ∗(ν∗− x)) + C2K0(2
δ∗(ν∗− x)) + ql, (50) where I0and K0are modified Bessel functions of the first and third kinds, respectively.
Substituting the obtained function qU into Eqn. (33) and using the well-known relations between the Bessel functions (Bateman, 1974), we find the function
jU = −L
ν∗− x δ∗
C1I1(2
δ∗(ν∗− x))
− C2K1(2
δ∗(ν∗− x)) ,
(51)
where I1and K1are modified Bessel functions. Note that, similarly to the previous case, the constants C1and C2can be determined from the boundary conditions. Omitting rather simple calculations, we present only the result:
C1=(q0− ql)K1(2
δ∗(ν∗− 1))
Δ , (52)
C2=(q0− ql)I1(2
δ∗(ν∗− 1))
Δ , (53)
where
Δ ≡ I0(2
δ∗ν∗)K1(2
δ∗(ν∗− 1)) + K0(2
δ∗ν∗)I1(2
δ∗(ν∗− 1)) and q0is defined by (41).
Thus, we have found the explicit formulae for qU
and jU. Having the formulae (50)–(51), the system of ODEs (35)–(36) can be reduced to two linear autonomous ODEs with the unknown functions u(x) and w(x). These equations have the form
d1(νmax− νmin)
(ν∗− x)d2u dx2 −du
dx
+√SG
δ∗
d dx
√ν∗− xΔ1(x)u
+
SGΔ0(x) − κ1
u − u01= 0
(54)
and
d2(νmax− νmin)
(ν∗− x)d2w dx2 −dw
dx
+√SA
δ∗
d dx
√ν∗− xΔ1(x)(w − w0)
+
SAΔ0(x) − κ2
(w − w0) − w01= 0, (55)
where
Δ1(x) ≡ C1I1(2
δ∗(ν∗− x))
− C2K1(2
δ∗(ν∗− x)) and
Δ0(x) ≡ C1I0(2
δ∗(ν∗− x)) + C2K0(2
δ∗(ν∗− x)).
Although both the equations are linear second order ODEs with the same structure, we could not find their general solutions because of their awkwardness. Thus, we solve them numerically together with the boundary conditions (26)–(27) using the Maple program package.
In the next section, realistic values of the parameters for the formulae (50)–(55) will be selected and applied in numerical simulations to calculate the non-uniform steady-state solutions.
4. Numerical results and their biomedical interpretation
Here we present numerical results based on the formulae derived in Section 3. Our aim is to check whether they are applicable for describing the fluid-glucose-albumin transport in peritoneal dialysis. The parameter values in these formulae were derived from experimental and clinical data and have been applied in previous mathematical studies (Van Liew, 1968; Imholz et al., 1998; Zakaria et al., 1999; Flessner, 2001; Smit et al.
2004a; 2004b; Waniewski, 2001; Stachowska-Pietka et al., 2006; 2007; Cherniha et al., 2007; Waniewski et al., 2009 ). Most of the parameters, especially those for water and glucose, were derived from experimental data or obtained by fitting the distributed model to clinical data, and are discussed in detail in a recent paper of Stachowska-Pietka et al. (2012).
Some phenomena, such as vasodilatation and change of the parameters of interstitial transport with the change of tissue hydration, were not included in our model because our objective was the mathematical analysis of the model, so its structure had to remain not too complicated. Nevertheless, the model covers all basic transport phenomena and provides a good background for further modifications. However, only numerical studies are available for such extended models (see Smit et al., 2004a; 2004b; Stachowska-Pietka et al., 2006). Furthermore, those parameters without firmly established experimental values (i.e., Staverman reflection coefficients), were varied to check their impact on the results of modelling . The diffusivity of albumin in the interstitium is not well known, but it is much lower than interstitial diffusivity of glucose (see the works of Waniewski (2001) and Stachowska-Pietka et al. (2012) for more details). The values of parameters and absolute constants applied in numerical simulations are listed in Table 1.
Let us consider the first case of a constant fractional fluid void volume, i.e., with the restrictions (31) and (37).
We remind the reader that the assumption that F is a constant was applied in many studies and this implies that ν(x) is also a constant. It seems to be reasonable to set νm = (νmax+ νmin)/2 = 0.26, i.e., we assume that the fractional fluid void volume at the steady-state stage of the peritoneal transport is an intermediate value between its maximum and minimum. In order to compare the numerical results obtained here with those for osmotic peritoneal transport obtained in earlier models, in which albumin transport was not considered, we neglect the oncotic pressure as a driving fluid force across the tissue, i.e., we set the Staverman reflection coefficients for albumin σT A = σA = 0. This means that the fluid flux across the tissue, jU, and the fluid flux from blood to the tissue, qU (see the formulae (4) and (5)), do not
depend on the albumin concentrations. The Staverman reflection coefficients for glucose in the tissue and in the capillary wall are equal to σT G= σG= 0.001. Hereafter, the values of other parameters and absolute constants are taken from Table 1.
Figure 1 presents the spatial distributions of the steady-state density of the fluid flux from blood to the tissue qU and the fluid flux across the tissue jU, calculated using the formulae (38)–(41). The negative sign of jU indicates the net fluid flux occurs across the tissue towards the peritoneal cavity. Therefore, it corresponds to the water removal by ultrafiltration. The monotonically decreasing (with the distance from the peritoneal surface) function qU(x) and the monotonically increasing function jU(x) are in agreement with experimental data and previously obtained numerical results for the models that took into account only the glucose transport (Cherniha et al., 2007; Waniewski et al., 2007). It should be stressed that in those previous models albumin transport was not considered and the restrictions (31) and (37) were not used.
Using the value of the fluid flux jUat the point x = 0, one may calculate the reverse water flow (i.e., out of the tissue to the cavity). The total fluid outflow from the tissue to the cavity (ultrafiltration), calculated assuming that the surface area of the contact between the dialysis fluid and the peritoneum is equal to 5 · 103cm2 (this surface area measured in 10 patients on peritoneal dialysis was found to be within the range from 0.41 to 0.76 m2 (Chagnac et al., 2002)), is about 0.90 mL/min. Note that a similar value was obtained previously by Cherniha et al., (2007) using numerical simulations. Moreover, it comes from the formula (4) for x = 0 that the ultrafiltration increases with growing σT G. For example, if one sets σT G = 0.01 into Eqn. (4) then the total fluid outflow from the tissue to the cavity is 5.2 mL/min, which is very close to the value obtained by Cherniha et al. (2007) for the same parameters.
Figure 2 presents the spatial distributions of the glucose concentration in the tissue for σT G = 0.001 and σT G = 0.01 (see Fig. 2, the top panel). The interstitial glucose concentration CG decreases rapidly with the distance from the peritoneal surface to the constant steady-state value of CG∗ (see the formula (15)) in the deeper tissue layer independently of the σT Gvalues and is practically CG∗ for x > 0.3 (see Fig. 2, the bottom panel, where both curves coincide). Thus, the width of the tissue layer with the increased glucose concentration (that is around 0.3 cm) does not depend on σT G. This remains in agreement with the previous results obtained by Cherniha et al. (2007).
We may conclude that, although the restriction in the form of the assumption (37) is rather artificial from physiological point of view, the analytical formulae derived in Section 3 lead to results which are similar to
those obtained earlier with numerical simulations of pure glucose and water peritoneal transport (Cherniha et al., 2007), where this assumption was not used.
Let us now consider the second case, which is more realistic, i.e., hereafter the restrictions (31) and (46) apply.
Remark 1. In the case σT A = 0.0 the results obtained via the formulae (50)–(53) and the ODEs (54)–(55) practically coincide with those presented above (see Fig. 1).
Now we assume that the Staverman reflection coefficient for albumin is non-zero and equal to σA = σT A = 0.5, i.e., the maximum value of σT A (see Table 1) is taken. In other words, we assume that the oncotic pressure plays an important role, contrary to the previous case. Here the fluid flux across the tissue jU and the fluid flux from blood to the tissue qU (see the formulae (4) and (5)) depend on the interstitial concentrations of glucose and albumin. We performed many calculations using the formulae (50)–(53) and ODEs (54)–(55) for a wide range of values of the parameter σT G, from very small (0.001) to large (0.03) values. Of course, other parameters can also be varied; however, we restricted ourselves to this parameter because it is included in assumption (31).
The results obtained for σT G= 0.001, σT G= 0.002 and σT G= 0.01 are presented in Figs. 3 and 4. It is quite interesting that the profiles for functions qU(x) and jU(x) shown in Fig. 3 are very similar to those in Fig. 1, although the relevant formulae are essentially different (the reader may compare (50)–(53) with (38)–(41)) and σT A = 0.5.
Moreover, the form of these profiles are the same for a wide range of the values of σT G.
Using the value of the fluid flux jUat the point x = 0, one may again calculate the ultrafiltration flow to the peritoneal cavity that can be obtained under the assumed restrictions on the model parameters. In the case σT G = 0.01 the total fluid outflow from the tissue to the cavity is approximately equal to 4.8 mL/min, while it is very small (0.06 mL/min) for σT G = 0.001. To obtain values of the ultrafiltration that correspond to those measured during peritoneal dialysis, we need to set σT G≥ 0.02. For example, setting σT G= 0.02 and σT G= 0.03, we obtain the ultrafiltration rates of 10 mL/min and 15 mL/min, which are close to those measured in clinical conditions for similar boundary concentration of glucose (Smit et al., 2004a; 2004b; Waniewski et al., 1996a; 2009; Stachowska et al., 2012). The initial values of ultrafiltration remain in agreement with those obtained by our group in clinical studies (Waniewski et al., 1996b; Waniewski, 2007). The initial rates of ultrafiltration for 3.86% glucose solution were found to be of 15 mL/min, which was much higher than those for 2.27 and 1.36% glucose solutions, 8 and 6 ml/min, respectively (Waniewski et al. 1996a; 1996b).
Similar values in the range of 14–18 mL/min were measured during the initial minute and much lower values
of 4-8 mL/min at the end of a 4 h dwell study with 3.86%
glucose solution (Smit et al. 2004a; 2004b).
Table 1. Parameters of the model used for numerical analysis of peritoneal transport.
Parameter name Parameter symbol,
value and unit
Minimal fractional void volume νmin= 0.17 Maximal fractional void volume νmax= 0.35 Staverman reflection coefficient 0 < σT G ≤ 0.03 for glucose
Sieving coefficient of glucose ST G= 1 − σT G
in the tissue
Staverman reflection coefficient 0.05 ≤ σT A≤ 0.5 for albumin
Sieving coefficient of albumin ST A= 1 − σT A
in the tissue
Hydraulic permeability of the tissue K = 5.14 · 10−5 cm2/ min /mm Hg Gas constant times temperature RT = 18 · 103
mm Hg· mL/ mmol Width of the tissue layer L = 1.0 cm Hydraulic permeability of LPa = 7.3 · 10−5 capillary wall times density 1/ min / mm Hg of capillary surface area
Volumetric fluid flux to ql= 0.26 · 10−4
lymphatic vessels min−1
Diffusivity of glucose in the tissue DG= 12.11 · 10−5
divided byνmin cm2/ min
Diffusivity of albumin in the tissue DA= 0.2 · 10−5
divided byνmin cm2/ min
Permeability of the capillary wall pGa =
for glucose times density 3.4 · 10−2min−1 of capillary surface area
Permeability of the capillary wall pAa = for albumin times density 6 · 10−5min−1 of the capillary surface area
Glucose concentration in blood CGB= 6 · 10−3 mmol / mL Albumin concentration in blood CAB= 0.6 · 10−3
mmol / mL Glucose concentration in dialysate CGD= 180 · 10−3
mmol / mL Albumin concentration in dialysate CAD= 0 Hydrostatic pressure of blood PB= 15 mm Hg Hydrostatic pressure of dialysate PD= 12 mm Hg Non-dimensional parameter α = 0.8
The spatial distributions of the glucose and albumin concentrations for different values of σT G are displayed in Fig. 4. Note that the glucose concentration in the tissue, CG, again decreases rapidly with the distance from the peritoneal cavity to the constant steady-state CG∗ in the deeper tissue layer. The glucose concentration is practically equal to CG∗ for any x > 0.1 cm if σT G is
