Mathematical Modeling of the Dynamic Exchange of Solutes in a Prototype Hemodialyzer

Hemodialysis (HD) is one type of procedure for eliminating toxic chemicals and infusing bicarbonate in patients with end-stage renal disease (ESRD). Research and development in the hemodialyzer industry have, hitherto, depended mostly on empirical evidence to optimize HD therapy. This is often costly and involves numerous clinical trials. Developing a comprehensive time-dependent mathematical model to examine the dynamic exchange of solutes (HCO-3 and pCO2), blood pH and H+ ions in a prototype hollow-fiber hemodialyzer is essential in optimizing future design and improvement. A comprehensive mathematical model which is represented by a coupled set of transport equations and delineates the blood and dialyzate compartments of the hemodialyzer, and includes bicarbonate-buffering reaction in the blood channel and bicarbonate replenishment mechanism in the dialysate, is used to describe the time-dependent bulk concentration and exit concentration of solutes, blood pH and H+ ions in the hollow-fiber prototype hemodialyzer. A numerical simulation of the model is used to test several time-dependent bulks and exit concentration profiles of solutes in the blood and dialyzate. Results obtained from the numerical solution of the model show the bulk and exit concentrations of solute at various distances along the blood and dialyzate channels at different times. This modeling exercise will also allow us in our next study to examine some physical mechanisms of the hemodialyzer.


Introduction
One of the leading goals of HD therapy, aside the elimination of electrolytes, How to cite this paper: Boamah, E.K. 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 [1] [2] [3]. 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. [4], 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. [5] 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. [6], 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. [7] 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 [8]- [13].
The only model based on the description of dialyzate bicarbonate concentration is the research by Heineken et al. [14]. They rearranged Sargent and Gotch's [15] 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.

E. K. Boamah
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.  The nonlinear reactive term which deals with buffering of the blood,

Model
The replenishment term which replenishes bicarbonate concentration in the dialysate.

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.
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 r T called the transmittance coefficient is the fraction of solutes that penetrate the membrane pores. Similarly, the convection-diffusion-reaction equa- 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

Nondimensionalization
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  ; ; are the dimensionless parameters. Similarly, on the dialyzate side, the non-dimensional parameters and the mass transport model is given as ; ;

Chemical Reactions
The most important buffering reaction in the blood is the inter-conversion of CO 2 and − 3 HCO . These undergo the following reversible reaction (catalyzed in the presence of carbonic anhydrase (C.A) or uncatalyzed)

Dialyzate Replenishment
The body HCO from a bath [17]. The addition of − 3 HCO 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 − 3 HCO 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 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).

Method
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

The Bulk Concentration
The bulk concentrations 1_

Clinical Data
The following tables provide initial data for pre-hemodialysis concentrations for − 3 HCO and CO 2 ( 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 − 3 HCO and CO 2 were obtained from literature [18] and [19]. Also permeabilities of − 3 HCO and CO 2 , fiber radius, membrane thickness, tube length for the model were based on a high-flux Fresenius type of hollow fibers from Gambro.

Results
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 − 3 HCO concentrations in the blood and dialysate were 19 mEq/l and 35 mEq/l respectively, and the respective inlet pCO 2 concentrations in the blood and dialyzate were 40 mmHg and 60 mmHg.   HCO 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 suit-able bulk exit concentration of  in the blood and dialyzer channels respectively at increasing selected times.

Discussions
As time progresses, the bulk concentrations of pCO 2 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 pCO 2 replenishment mechanism of the dialysate. The graph in Figure 3 also shows that the optimal concentration of pCO 2 is reached at the exit of the blood channel. It is crucial to obtain suitable bulk exit concentration of pCO 2 because this is the concentration that enters the body of the ESRD patient. In Figure 4 however, we again see hypothetical representations of dialysate pCO 2 over the hollow-fiber length as time increases. The decline in dialysate pCO 2 concentration is most rapid in the far end of the tube. The more rapid decline is due to the initially low pCO 2 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     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.

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 [20].
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 [21], 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.

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.