One of the leading goals of HD therapy, aside the elimination of electrolytes, toxic chemicals, and water, is correcting metabolic acidosis by the infusion of bicarbonate as a buffer from the dialyzate into the bloodstream in a prototype hemodialyzer. Normalization of metabolic acidosis sometimes results in, for example, metabolic alkalosis due to over compensation of the acidosis, causing dialysis symptoms such as mental confusion and muscle cramps   . Thus, it is of vital importance to obtain the closest if not the exact correction of metabolic acidosis during HD as time progresses.
Research and development in the hemodialyzer technology and HD therapy have hitherto depended mostly on empirical evidence. This is costly and often involves numerous clinical trials. Gotch et al. , conducted clinical experiments on eight ESRD patients. In this experiment, the researchers used 36 mEq/l of bicarbonate in the dialyzate, however, results were not good in terms of acid-base balance. In a similar experiment performed over 12 weeks, Ward et al.  were faced with the problem of how to control bicarbonate during HD therapy and that resulted in abnormally high pH value which puts the patient at risk of post-dialytic alkalosis. To address the issue of metabolic acidosis, La Greca et al. , used 31 mEq/l of bicarbonate and 5 mEq/l of acetate in the dialyzate during HD therapy. Here, the abnormally low pH values as observed in the ESRD patients were not adequately corrected. In an attempt to address the discrepancies in the experimental results, Gotch et al., considered a black box input-output ordinary differential equation model to describe the mass balance of hydrogen ions during HD. Black box model in this case means that the internal variables of the hollow fibers were not taken into consideration. The model by Gotch et al., failed to predict the acid-base state during a single dialysis session because several processes controlling solute concentrations were not taken into account. Due to variations that were still observed in the experimental results of HD procedures, there have been increasing need to consider mathematical models. Monti et al.  also used a black box input-output ordinary differential equation model to describe the dynamic exchange process of bicarbonate during and after HD. This black-box model was oversimplified because the internal structure and variables of the prototype hemodialyzer unit were not taken into account. Thus, the model did not allow for the prediction of internal variables such as the surface area of the dialyzer, dialyzate flow rates, the thickness and permeability of the membrane  - .
The only model based on the description of dialyzate bicarbonate concentration is the research by Heineken et al. . They rearranged Sargent and Gotch’s  modified mass balance equation to determine the dialyzate bicarbonate concentrations for buffering hydrogen ions during HD. Although Heineken et al., showed the benefits of long-term dialyzate bicarbonate concentration, their black box model did not take into account the exchange process inside the hemodialyzer unit. Thus, the model failed to predict bicarbonate dynamics over the course of a single HD procedure. Hemodialysis procedure is still problematic and poses challenging questions of great importance.
In order to go beyond black-box input-output ordinary differential equation models, and efficiently address the dynamic exchange of solutes in the prototype hemodialyzer, a comprehensive partial differential equation model incorporating internal properties, the dimensions and the surface area of the hollow fibers used during hemodialysis, the nature of solute transfer across the hemodialyzer semi permeable membrane, the hemodialysis flow rates and the duration of the administration of hemodialysis were taken into account. The knowledge gained from this study will eventually lead to a means of predicting the underlying mechanisms of solute transfer in a prototype hemodialyzer, thereby minimizing the need to perform costly and time-consuming clinical trials.
2. Model Description
Concentration of species in the blood,
Concentration of species in the dialyzate,
Velocity of blood in the axial direction,
Velocity of dialyzate in the axial direction,
Velocity of blood in the radial direction,
Velocity of dialyzate in the radial direction,
Wall velocity in the membrane-blood channel,
Wall velocity in the membrane-dialyzate channel,
Radius of the blood channel,
Radius of the dialyzate channel,
Sum of and ,
Blood-membrane Sherwood number,
Dialyzate-membrane Sherwood number,
Length Pecle’t number (blood side),
Radial Pecle’t number (blood side),
Length Pecle’t number (dialyzate side),
Radial Pecle’t number (dialyzate side),
Coaxial entrance length,
Radius of the annulus,
Membrane permeability of species,
Maximum axial velocity of blood flow,
Maximum axial velocity of dialyzate flow,
Dimensionless time scale,
Initial concentration of species in the blood,
Initial concentration of species in the dialyzate,
Pressure at the entrance of the annular,
Pressure at the end of the annular,
Diffusion coefficient of species in the radial direction,
Axial grid size,
Radial grid size,
Number of tubes,
Blood flow rate,
Dialyzate flow rate,
Ultrafiltration flow rate,
Mass flux in the blood side channel,
Mass flux in the dialyzate side channel,
The nonlinear reactive term which deals with buffering of the blood,
The replenishment term which replenishes bicarbonate concentration in the dialysate.
2.2. Model Assumptions
We will use the following fundamental assumptions to formulate the models.
· We consider dilute, aqueous and Newtonian flow mechanism in the hemodialyzer.
· We consider permeability, diffusivity and density as constants.
· Axial diffusion was considered negligible.
· Flows were considered laminar.
· The flow mechanisms for the blood and the dialyzate are countercurrent.
· There are no angular gradients (axis-symmetry approximation).
The rate of change of species concentration in the blood and dialyzate channels are given as
The initial and boundary conditions for Equations (2.1) to (2.3) are:
where called the transmittance coefficient is the fraction of solutes that penetrate the membrane pores. Similarly, the convection-diffusion-reaction equation on the dialyzate side is given as
with respect to the axial and the radial velocity profiles,
where is the viscosity of dialyzate. The initial and boundary conditions for (2.8) to (2.11) are given as
We nondimensionalize the models and drop asterisks for simplicity to obtain the following nondimensional model and parameters. The dimensionless mass transport model in the blood channel is given as
The initial and boundary conditions for (2.16) to (2.18) are given as
and are the dimensionless parameters. Similarly, on the dialyzate side, the non-dimensional parameters and the mass transport model is given as
and are the dimensionless parameters
2.4. Chemical Reactions
The most important buffering reaction in the blood is the inter-conversion of CO2 and . These undergo the following reversible reaction (catalyzed in the presence of carbonic anhydrase (C.A) or uncatalyzed)
where and are the forward and reversible reaction rate constants respectively. Chemically, (Reaction 2.31) forms a major buffer system of the blood. The net rate of reaction of species by chemical (reaction 2.31) per unit volume can be expressed as
The interrelationship between pH, , and pCO2 is then expressed by the Hendernson-Hasselbalch equation as
where is the partial pressure of carbon dioxide and is the solubility constant in mmol/l . The concentration of H+ is obtained from the pH of the blood by the equation,
By (Reaction 2.31),
We simplify and nondimensionalize (2.35) and (2.36) to obtain the following equations:
where and .
2.5. Dialyzate Replenishment
The body stores are replenished by adding from a bath . The addition of is partially consumed by titration with the body buffers and by organic acid production. Bicarbonate remaining is then added to the pool and that causes to rise in the blood returning to the hemodialysis membrane. In our model, we used a differential rate equation to determine the total amount of used for replenishment in the dialyzate during the time of treatment. A suitable nonlinear equation is given as
In nondimensional form, we obtain the equation,
In all, the problem consists of four nonlinear dimensionless equations. Two equations were each obtained from the blood and dialyzate channels (Figure 1).
3. Numerical Formulation
The method of lines (MOL) is employed to reduce the model equations to a system of ordinary differential equations. Basically, the system of equations obtained from our model equations as a result of the application of the MOL using central difference, backward difference and up-winding schemes is given as
Figure 1. This plot shows the transient boundary and bulk concentration profiles of bicarbonate in the blood as a function of the tube length after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins. The inlet bicarbonate concentration was 19 mEq/l and the blood and dialyzate flow rates were 300 ml/min and 500 ml/min respectively.
and are matrices of dimensions and respectively. and are respectively the blood side, the dialyzate side coupling matrices, the dialyzate side and the blood side coupling matrices, each of dimension . The vectors u, v and w are each of dimension .
3.2. The Bulk Concentration
The bulk concentrations and which is the most meaningful concentration measured at the outlet of the blood channel. These are given as
3.3. Clinical Data
The following tables provide initial data for pre-hemodialysis concentrations for and CO2 (Table 1) and values for constants and operational hemodialysis parameters (Table 2) for running the numerical simulations. The listed initial data for pre-hemodialysis for and CO2 were obtained from literature  and . Also permeabilities of and CO2, fiber radius, membrane thickness, tube length for the model were based on a high-flux Fresenius type of hollow fibers from Gambro.
Table 1. Initial hemodialysis concentrations used in the numerical simulations.
Table 2. Constants and operational hemodialysis parameters used in the numerical simulations.
4. Results and Discussions
Following the discussion on the non-steady state model system, we present the results of the numerical simulations. The transient concentration profiles for different time steps along the tube length are shown in Figure 1 to Figure 6 below. The inlet concentrations in the blood and dialysate were 19 mEq/l and 35 mEq/l respectively, and the respective inlet pCO2 concentrations in the blood and dialyzate were 40 mmHg and 60 mmHg.
Figure 2. The plot shows the time evolution of boundary and bulk concentration profiles of bicarbonate in the dialyzate as a function of the tube length after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins.
Figure 3. This plot shows the transient boundary and bulk concentration profiles of partial pressure of carbon dioxide in the blood as a function of the tube length after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins.
As time progresses, the bulk concentrations of in mEq/l increases in the blood side channel (Figure 1) and stays approximately constant in the dialysate side channel which is expected due to the replenishment mechanism of in the dialysate. The graph in Figure 1 also shows that the optimal concentration is reached at the exit of the blood channel. It is crucial to obtain suitable bulk exit concentration of because this is the concentration that enters the body of the ESRD patient. In Figure 2 however, we see hypothetical representations of dialyzate over the hollow-fiber length as time increases. The decline in dialysate concentration is most rapid in the far end of the tube. The more rapid decline is due to the initially low concentration gradient between the blood and the dialyzate.
As time progresses, the bulk concentrations of pCO2 in mmHg increases in the blood side channel (Figure 3) and stays approximately constant (Figure 4) in the dialysate side channel which is expected due to the pCO2 replenishment mechanism of the dialysate. The graph in Figure 3 also shows that the optimal concentration of pCO2 is reached at the exit of the blood channel. It is crucial to obtain suitable bulk exit concentration of pCO2 because this is the concentration that enters the body of the ESRD patient. In Figure 4 however, we again see hypothetical representations of dialysate pCO2 over the hollow-fiber length as time increases. The decline in dialysate pCO2 concentration is most rapid in the far end of the tube. The more rapid decline is due to the initially low pCO2 concentration gradient between the blood and the dialyzate.
In Figure 5, the pH increases in the blood channel as time progresses.
The normal pH of the body lies in the range of 7.35 to 7.45. As seen from the graph in Figure 5, the pH almost normalizes at the exit of the blood channel which is crucial for the ESRD patient. As increases in the blood, hydrogen ions are reduced as seen in Figure 6.
Figure 4. This plot shows the transient boundary and bulk concentration profiles of partial pressure of carbon dioxide in the dialyzate as a function of the tube length after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins.
Figure 5. The plot shows the transient boundary and bulk concentration profiles of pH value in the blood as a function of the tube length after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins.
Figure 6. The plots in this figure show the boundary and bulk concentration profiles respectively of H+ concentrations in the blood channel after time steps: 1.8 s, 3.0 s, 4.5 s, 7.2 s, 13.7 s, 23.8 s, 40.1 s, 1.7 mins, 2 mins, 6.8 mins, 32.7 mins and 240 mins.
Figure 6 shows that hydrogen ion concentration is at a minimum at the exit of the blood side channel. From the graphs (Figure 1 to Figure 6), as time increases, there are no significant changes in concentrations a few minutes after the dialysis procedure. The reaction process is very fast and almost reaches steady state between the time intervals from 1.8 s to 13.7 s. We noticed that the uppermost curves for , and pH value in the blood channel for the interval from 23.8 s to 240 mins are all overlapping each other. In a similar situation, the lowest curve for , in dialyzate and H+ ions in the blood overlap each other at the given times. At this point, we are not sure of the cause.
In this paper, we have used mathematical models to study the dynamic exchange of solutes during HD. This model considers replenishing bicarbonate concentration in the dialyzate compartment and bicarbonate buffering mechanism in the blood compartment. In general, we have simulated the profiles of solute concentrations inside the prototype hemodialyzer and have used known techniques to develop numerical algorithms for the non-steady state model of the prototype hemodialyzer.
We have developed mathematical models that are more comprehensive than the models of Gotch et al. and Monti et al. The models of Gotch and coworkers and Monti and coworkers describe black box models. Such models do not represent the structure of HD, neither do they predict variations in solute concentrations internal to the prototype hemodialyzer.
In our HD model, the dynamic exchange mechanism of the solutes (bicarbonate, carbon dioxide, and hydrogen ions) in the hemodialyzer has been described by partial differential equations using concepts from computational fluid dynamics. Using the model simulations and graphical representations (Figure 1 to Figure 6), we have been able to examine the underlying mechanisms, and how the effects of the changes in various parameters affect the dynamic exchange of solutes in the hemodialyzer.
4.3. Model Limitations
The presented model is based on a number of assumptions that might limit its direct applicability to HD. However, results obtained do reveal interesting underlying mechanisms that contribute to the current body of knowledge on HD therapy. In addition, it provides an opportunity for further expansion and development. Although the most important blood buffering process is the inter-conversion of bicarbonate-carbon dioxide, it might make a difference to incorporate non-bicarbonate buffers such as albumin, ammonia, phosphate and hemoglobin into the model .
During HD therapy, a decrease in the overall permeability which is caused by the adhesion of high protein adsorption on the membrane surface or the occlusion of higher molecular weight solutes needs to be investigated because as the adhesions increase, membrane resistance increases and, thus, permeability decreases. The adhesions may also impede adequate flow of solutes during HD. In order to observe the overall effect of adhesions on the membrane surfaces during the exchange mechanism of solutes, the assumption of constant permeability needs further investigation. My future work will include a modified permeability coefficient.
In the model system, the effects of electrolytes such as potassium, sodium and chloride in the blood during HD were ignored. In reality, these solutes might impact the results of the solutions to the models. These components will be investigated in my future work. Also, the charge on the semi-permeable membrane was assumed to have no influence on the dynamic exchange. In reality, Membrane charges and the mass transfer of ionic solutes may be influenced by the transmembrane gradient of electrical potential due to Gibbs-Donnan effect , and may cause adsorption of solutes. In my future work, I will extend the model to incorporate the effect of membrane and solute charges. Despite these limitations, the model does provide a very good starting point for further investigations.
5. Concluding Remarks
We have developed algorithms and solved the important case of time dependence. In the case of the time dependent problem, the full system was solved and numerical results obtained using MATLAB. The time dependent model system and simulation results offer the flexibility to vary and control the dialyzate bicarbonate concentration during HD. The transfer of bicarbonate during HD also depends on the type of the hemodialyzer, flow rates, duration of HD, and dialyzate bicarbonate concentration. We observed that the amount of bicarbonate transferred, for example, might be best achieved by adjustments in the dialyzate bicarbonate concentration.
In summary, we have described the underlying mechanisms of the dynamic exchange of solutes during HD therapy. Numerical simulations performed give valuable insight to concentration distributions of the solutes (bicarbonate, partial pressure of carbon dioxide and hydrogen ions) and the outlet pH values. This study has helped to provide valuable insights into the flow dynamics and the mass transfer behavior of the solutes and the pH values in the prototype hemodialyzer during HD therapy. We expect that the simulation systems coupled with our model refinements will be useful for fully understanding the nature of solute transfer during HD.
The insight gained from this simulation will eventually pave the way to providing suggestions for clinical HD therapy and optimizing future design and improvement of the hemodialyzer. The model can be used to determine the individual dialyzate bicarbonate concentration for the ESRD patients, to choose the type of hemodialyzer for an ESRD patient and to determine the outlet blood bicarbonate concentration.
The problem of convection, diffusion, and reaction of the solutes between the blood and the dialyzate during HD therapy is complex and as a result, further work will be needed to identify predictive tools necessary for hemodialysis adequacy. Also, experimental studies might be needed to measure how solute adhesions on the semi-permeable membrane impede adequate flow, and treatment duration for optimal HD therapy.
This work was as a result of my Ph. D dissertation, consequently, the author acknowledges Dr. Daniel Bentil (my Ph. D dissertation advisor) of the University of Vermont for his helpful discussions during my Ph. D dissertation write-up. The author also acknowledges with gratitude the financial support provided by the University of Vermont during my Ph. D work.
 Ahmad, S., Pagel, M., Vizzo, J. and Scribner, B.H. (1980) Effect of the Normalization of Acid-Base Balance on Postdialysis Plasma Bicarbonate. Transactions—American Society for Artificial Internal Organs, 26, 318-340.
 Raja, R.M., Kramer, M.S., Rosenbaum, J.L., Fernandes, M. and Barber, K. (1982) Hemodialysis with Varying Dialysate Bicarbonate Concentration. Transactions—American Society for Artificial Internal Organs, 28, 514-521.
 La Greca, G., Fabris, A., Feriani, M., Chiaramonte, S. and Ronco, C. (1989) Acid-Base Homeostasis in Clinical Dialysis. 3rd Edition, In: Maher, J.F., Ed., Replacement of Renal Function by Dialysis, Kluwer Academic Publishers, Boston, 808-826.
 Monti, J.-P., Sarrazin, M.Y., Baz, M., Murisasco, A., Crevat, A.D. and Elsen, R. (1990) Kinetic Modeling of Intradialytic and Interdialytic pH Shifts during and after Acetate and Bicarbonate Hemodialysis. Artificial Organs, 13, 799-802.
 Ronco, C., Brendolan, A., Crepaldi, C., Rodighiero, M. and Scabard, M. (2002) Blood and Dialyzate Flow Distributions in Hollow-Fiber Hemodialyzers Analyzed by Computerized Helical Scanning Technique. Journal of the American Society of Nephrology, 13, S53-S61.
 Eloot, S., De Wachter, D., Van Tricht, I. and Verdonck, P. (2002) Computational Flow Modeling in Hollow-Fiber Dialyzers Artificial Organs, 26, 590-599.
 Chelamcharla, M., Leypoldt, J.K. and Cheung, A.K. (2005) Dialyzer Membranes as Determinants of the Adequacy of Dialysis. Seminars in Nephrology, 25, 81-89.
 Sato, Y., Mineshima, M., Ishimori, I., Kaneko, I., Akiba, T. and Teraoka, S. (2003) Effect of Hollow Fiber Length on Solute Removal and Quantification of Internal Filtration Rate by Doppler Ultrasound. The International Journal of Artificial Organs, 26, 129-134.
 Heineken, F.G., Brady-Smith, M., Haynie, J. and Van Stone, J.C. (1988) Prescribingdialyzate Bicarbonate Concentrations for Hemodialysis Patients. The International Journal of Artificial Organs, 11, 45-50.
 Atherton, L.J., Maddox, D.A., Gennari, F.J. and Deen, W.M. (1988) Analysis of PCO2 Variations in the Renal Cortex. I. Single Nephron. American Journal of Physiology, 255, F349-F360.
 Ronco, C., Bowry, S.K., Brendolan, A., Crepaldi, C., Soffiati, G., Fortunato, A., Bordoni, V., Granziero, A., Torsello, G. and La Greca, G. (2002) Hemodialyzer: From Macro-Design to Membrane Nanostructure; the Case of the FX-Class of Hemodialyzers. Kidney International, 61, S126-S142.
 Ronco, C., Orlandini, G., Brendolan, A., Lupi, A. and La Greca, G. (1998) Enhancement of Convective Transport by Internal Filtration in a Modified Experimental Hemodialyzer. Kidney International, 54, 979-985.
 Kilmartin, J.V. and Rossi-Bernardi, L. (1973) Interaction of Hemoglobin with Hydrogen Ions, Carbon Dioxide and Organic Phosphates. Physiological Reviews, 53, 836-890.
 Petitclerc, T. (1998) Estimation of Mass Transfer through a Hemodialyzer: Theoretical Approach and Clinical Applications. Artificial Organs, 22, 601-607.