Experimental evaluation of lime spray drying for SO₂ absorption

1.   Introduction

End-stage renal disease is a debilitating disease affecting approximately 4.9–9.7 million patients worldwide, with only 2.6 million patients receiving treatment ^[1,2]. The best solution currently available is kidney transplantation but not all patients are eligible for transplantation and there is limited organ availability. In turn, most patients are recommended dialysis as a bridge-to-transplantation, which seeks to reduce the renal load by artificially removing toxins and injecting solutes back to the blood ^[2]. However, dialysis covers only a small part of the kidney function, i.e., it mimics the kidney's glomerulus filtration function and, as such, only removes small and middle molecule uremic toxins ^[3,4]. The physiological function of the proximal tubule, namely, to actively remove the larger protein-bound uremic toxins, is not recapitulated in dialysis. These toxins then accumulate, resulting in anemia, insulin resistance, epileptic seizures, and renal failure, among many other bodily dysfunctions ^[5,6].

Microstructural (membrane porosity, composition and thickness) and macrostructural (overall membrane configuration, flow rate, blood thinners, membrane surface area) properties of the dialyser are being consistently analysed to improve the dialyser performance. One of the proposed microstructural solutions is to introduce wavy undulations to reduce deposition of solutes on the straight fibers ^[7]. Wavy fibers were simulated to observe the stress distribution and it was shown that crimped fibers experience less stress than straight ones ^[8]. This is corroborated by the ex vivo study performed by Leypoldt et al., where they show that mass transfer for low and middle weight molecules is increased by utilizing dialysers with wavy membranes ^[9].

Another way to improve the removal of free and protein-bound uremic toxins is by enhancing classic hemodialysis (referred to herein as 'dialysis' for simplicity, though the authors recognize other forms of dialysis exist) by culturing a living cell layer on top of the hollow fiber membranes ^{[10,11,12,13,14,15]}. This cell layer then contributes to the active removal of protein-bound toxins through dedicated transporters (e.g., organic anion transporter 1 (OAT1)) on their basolateral and apical surface. For example, Jansen et al. cultured conditionally immortalized proximal tubule epithelial cells with the overexpression of OAT1 (ciPTEC-OAT1) in vitro on dialysis hollow-fiber membranes. Their results show that the ciPTEC-OAT1 cells formed a functional barrier with the ability to clear indoxyl sulfate and kynurenic acid, two important protein-bound uremic toxins, and resorb albumin from the flow chamber ^[16].

It is also known, however, that such cell layers are fragile, and wall-shear stress at flow conditions common in commercial dialysis devices may damage or destroy the cell monolayer ^[17,18]. The use of a hollow fiber incorporating wavy undulations to reduce shear stress within the troughs may help to prevent such effects, allowing for the application of a cell monolayer in commercial dialysis fibers. Additionally, Hu et al. have shown that cell growth is promoted on wavy undulations ^[19], and Sheng et. al. found benefits to renal cell function ^[20] giving further justification for their use. Still, optimizing such a device is far from trivial due to the non-intuitive spatiotemporal dynamics of the toxin removal process. Mathematical modeling can address this challenge by quantitatively investigating the relation between the behavior of the cell monolayer under various flow conditions and the resulting toxin clearance.

In this present work, we introduce a model for such wall-shear stress effects that may simulate the complex interaction of wall-shear stress on ciPTEC cells cultured on a wavy hollow fiber membrane, and the resulting impact on toxin clearance. The article is outlined as follows. We begin by introducing the mathematical model and important notational aspects. We proceed to formally analyze the introduced model at the continuous level, establishing its well-posedness. We then seek to quantify how the wavy geometry influences toxin clearance, and in particular what is the optimal configuration of waves in terms of their number and size. We expect to see regions of increased wall shear transport near the sides of the waves, where the shear-stress is higher, with the troughs preventing the destruction of cells by preventing the wall shear stress from reaching excessively high levels. Finally, we conclude by summarizing the important findings and discussing directions for renal replacement therapies.

2.   Materials and methods

2.1.   Basic notation and mathematical preliminaries

Throughout this work, we denote an inner product of two scalar functions $f, \, g$ and vector functions $\boldsymbol{f}, \, \boldsymbol{g}$ defined on a domain $\Omega$ in the following way:

An inner product taken over a portion of a boundary (or portion of a boundary) $\Gamma$ is defined similarly:

Square integrable functions in $\Omega$ are denoted by the space $L^2(\Omega)$ with the norm $\sqrt{(f, \, f)}$. The Sobolev space of functions whose square is integrable in $\Omega$ together with the derivatives up to the order $s$ ($s \in \mathbb{N}$) are denoted by $H^s(\Omega)$. The definition can be extended to $s$ real non-negative numbers ^[21]. The corresponding spaces for vector-valued functions will be denoted by bold-face fonts. For a $H^1(\Omega)$ function $f$, the norm reads

If $f$ is a function in $H^p(\Omega)$ $p\geq 1$, and $\Gamma \subset \partial \Omega$ is a portion of the boundary of $\Omega$ with positive measure, we recall that there exists a bounded trace operator $\gamma$ extracting the value of the function on $\Gamma$ in the space by $H^{p-1/2} (\Gamma)$ ^[22]. The trace inequality states in particular that there exists a constant $\beta_t$ such that:

for all $\phi$ in $H^1(\Omega).$ With a popular notation, the space of $H^p$-functions with null trace on $\Gamma$ will be denoted by $H^p_{\Gamma}(\Omega)$. One may refer to ^[23] for more information on these topics.

2.2.   Problem definition

2.2.1.   Membrane transport

We consider a cylindrical domain $\Omega$ such that $\Omega = \Omega_{ \boldsymbol{b}} \cup \Omega_{ \boldsymbol{d}}$, $\Omega_{ \boldsymbol{b}} \cap \Omega_{ \boldsymbol{d}} = \emptyset$. A visual depiction of the model is in Figure 1. The flow of dialysate occurs in $\Omega_{ \boldsymbol{d}}$, and the flow of blood in $\Omega_{ \boldsymbol{b}}$. Both flows are modeled with a linear viscosity (denoted by $\nu_{ \boldsymbol{d}}$ and $\nu_{ \boldsymbol{b}}$ respectively). Let ${{\mathbb T}}$ be a generic tensor and $\mathbb{S}({{\mathbb T}}) \equiv 1/2({{\mathbb T}}+{{\mathbb T}}^T)$ its symmetric part. Then, we introduce the tensor:

where ${{\mathbb I}}$ denotes the identity tensor. The subdomains are separated by a membrane $\Gamma_m$ that is assumed to be not permeable to the fluids. The case where the membrane is semi-permeable is also important and a possible subject of future work. Assuming a steady flow, fluid velocities $\boldsymbol{u}_{ \boldsymbol{d}}$ and $\boldsymbol{u}_{ \boldsymbol{b}}$, and pressures $p_{ \boldsymbol{d}}$ and $p_{ \boldsymbol{b}}$ are described by the steady incompressible Navier-Stokes equations:

where $\rho_{ \boldsymbol{d}}$ and $\rho_{ \boldsymbol{b}}$ are the densities of dialysate and blood respectively, and $\boldsymbol{f}_{ \boldsymbol{d}}$, $\boldsymbol{f}_{ \boldsymbol{b}}$ possible body forces (e.g., the gravity). The equations are completed with the following boundary conditions:

where $\boldsymbol{n}_{ \boldsymbol{b}}$ ($\boldsymbol{n}_{ \boldsymbol{d}}$) is as usual the outward-pointing normal to $\Omega_{ \boldsymbol{b}}$ ($\Omega_{ \boldsymbol{d}}$). Notice that $\boldsymbol{n}_{ \boldsymbol{b}} = - \boldsymbol{n}_{ \boldsymbol{d}}$ on $\Gamma_m$. Also, $\boldsymbol{u}_{ \boldsymbol{b}, IN}$ is an inflow function in the blood domain, while $\boldsymbol{u}_{ \boldsymbol{d}, IN}$ a similar inflow function in the dialysate domain.

In the dialysis process, the solute $c$ is carried out of the blood across $\Gamma_m$ into the dialysate. The filtration rate is described by the flux $q$ of the solute across the membrane. We assume therefore that the solutes $c_{ \boldsymbol{b}}$ and $c_{ \boldsymbol{d}}$ in their respective domains $\Omega_{ \boldsymbol{b}}$, $\Omega_{ \boldsymbol{d}}$ obey an advection-diffusion problem where the advection field is given by the dialysate and blood velocities $\boldsymbol{u}_{ \boldsymbol{b}}$ and $\boldsymbol{u}_{ \boldsymbol{d}}$ respectively. The diffusivity coefficients will be denoted by $\mu_{ \boldsymbol{b}}$ and $\mu_{ \boldsymbol{d}}$ respectively and will be assumed to be constant too. In general, also these coefficients may depend non-linearly on the concentrations $c_{ \boldsymbol{b}}$ and $c_{ \boldsymbol{d}}$ respectively. This option will be considered elsewhere. The solutes obey therefore the following equations:

Here, $f_{ \boldsymbol{b}}$ and ${f}_{ \boldsymbol{d}}$ refer to appropriate body terms in the blood and dialysate domains as indicated by the subscripts. The interface condition (2.8c) is defined by a function $\widetilde{q}(c_{ \boldsymbol{d}}, c_{ \boldsymbol{b}})$ of the concentrations across the membrane. This function is a major ingredient of our model, so we discuss it in the next section.

2.2.2.   Modelling cell-layer transport

The most important aspect of the model (2.6, 2.8) is the definition of the transport across the membrane. The simplest approach would be a diffusive model

where $K$ represents the porosity of the membrane. Such a model was used, for example, for solute transport across blood vessel walls in ^[24]. Here we are interested in examining the effects of cultured kidney cell-layer aided dialysis, which feature more sophisticated transport dynamics and offer a promising improvement over existing membrane-only techniques ^{[13,16,25,26]}. In particular, the model (2.9) for $K$ constant does not account for the role of wall shear stress (which may be hereafter abbreviated as WSS) in the blood, defined as

In reality, WSS has a complex effect on solute transport. It has been experimentally shown that shear stress improves renal epithelial characteristics, such as epithelial polarization and transport function ^{[17,27,28,29]}. In contrast, high shear stress levels (>0.2–0.6 Pa) reduce the expression of epithelial characteristics, including ZO-1, E-cadherin and tubular cilia ^[17] as well as the cell viability ^[17,18]. We therefore propose the following modification of (2.9) incorporating this phenomenon. The wall shear stress-induced active transport model is defined by the following porosity

where $\| \cdot \|$ is understood to mean the standard Euclidean norm, $\widetilde{\tau}$ is a threshold value for the Euclidean norm of the WSS and $\widetilde{K}$ and $\xi$ are constants depending on the membrane. The model (2.11) responds to the necessity to acquire an optimal shear stress for maximum toxin clearance and it is the novel contribution of this work. We note that this model assumes adequate shear stress for the transporter expression, and for $\| \boldsymbol{\tau} \| = 0$ (corresponding to the situation $\boldsymbol{u}_{ \boldsymbol{b}} = 0$) we will have no transport. Further experimental validation is necessary to clarify this phenomenon and improve upon this modeling assumption.

3.   Results

3.1.   Well-posedness analysis of the WSS-dependent model

Before we consider the numerical approximation of our model, we first establish the well-posedness of the WSS-dependent model (2.8, 2.11) model, assuming that the fluid model (2.6) is well posed. Also, we will assume the fluid flow to be regular enough, and specifically that both the velocities $\boldsymbol{u}_{ \boldsymbol{b}}$ and $\boldsymbol{u}_{ \boldsymbol{d}}$ are nonzero and $\textbf{H}^2$ regular in their respective domains. Regularity conditions on the data and the domains that obtain this regularity can be found in ^[30]. For ease of notation, denote the function spaces:

as well as the bilinear forms:

The variational form of the problem then reads: find $c_{ \boldsymbol{d}}$ in $V_{ \boldsymbol{d}} \oplus \mathcal{L}_{ \boldsymbol{d}} (\widetilde{c}_{ \boldsymbol{d}})$, $c_{ \boldsymbol{b}}$ in $V_{ \boldsymbol{b}} \oplus \mathcal{L}_{ \boldsymbol{b}} (\widetilde{c}_{ \boldsymbol{b}})$ such that:

Here, $\mathcal{L}_{ \boldsymbol{d}} (\widetilde{c}_{ \boldsymbol{d}})$ and $\mathcal{L}_{ \boldsymbol{b}} (\widetilde{c}_{ \boldsymbol{b}})$ denote appropriate lifting operators of the inlet boundary conditions in the domains $\Omega_{ \boldsymbol{d}}$ and $\Omega_{ \boldsymbol{b}}$, respectively. In the subsequent analysis, to reduce technicalities, we will assume to work with homogeneous inlet conditions.

To carry out the analysis of the coupled problems (3.4), (3.5), we define the functional space: $\boldsymbol{V} : = V_{ \boldsymbol{d}} \, \times V_{ \boldsymbol{b}}$, endowed with the norm:

Adding together (3.4), (3.5), we obtain a new formulation of the problem in terms of the unknown vector $\bf{C}: = \lbrack c_{ \boldsymbol{d}}, \, c_{ \boldsymbol{b}}\rbrack^T$, forcing vector $\boldsymbol{F}: = \lbrack f_{ \boldsymbol{d}}, \, f_{ \boldsymbol{b}}\rbrack^T$, and test vector $\boldsymbol{\Phi} : = \lbrack \phi_{ \boldsymbol{d}}, \, \phi_{ \boldsymbol{b}}\rbrack^T$:

Find $\bf{C}$ in $\boldsymbol{V}$ such that:

for all $\boldsymbol{\Phi}$ in $\boldsymbol{V}$, where

Lemma 1. For $\widetilde{\tau} > 0$, $\boldsymbol{u}_{ \boldsymbol{b}} \in H^2$, $\widetilde{K} > 0$ the function $(2.11)$ is positive and essentially bounded over $\lbrack 0, \, \infty)$.

Proof. The positivity follows from the fact that $\| \boldsymbol{\tau}\| > 0$, and $\widetilde{K} > 0$. Then, it is enough to observe that the function

for $g(x, \xi, \widetilde{\tau}) \equiv \exp\left(-\xi(x-\widetilde{\tau})\right)$ has a (finite) maximum in $x^* > 0$, the root of the equation $g(x, \xi, \widetilde{\tau}) = \xi x -1$. The positivity of $x^*$ is promptly verified by noting that $g(0, \xi, \widetilde{\tau}) > 0$ while the right hand side is negative in $0$, and $g$ is monotonically decreasing, while $\xi x -1$ is monotonically increasing. In fact, we can state the bound

Lemma 2. If the data and the domains in $(2.6)$, $(2.7)$ are regular enough so that $\boldsymbol{u}_{ \boldsymbol{b}}$ and $\boldsymbol{u}_{ \boldsymbol{d}}$ are $H^2$-regular, then the bilinear forms $(3.2)$, $(3.3)$ are both continuous and coercive. Also, If $f_{ \boldsymbol{d}} \in L^2(\Omega_{ \boldsymbol{d}})$ and $f_{ \boldsymbol{b}} \in L^2(\Omega_{ \boldsymbol{b}})$, then the right hand sides of $(3.2)$, $(3.3)$ are continuous.

Proof. The continuity of the bilinear forms and of the functionals on the right hand side is a trivial consequence of the Cauchy-Schwarz and the trace inequality. For the bilinear form, we use the bound found in Lemma 1. In particular,

where $\beta_e$ comes from the embedding and $\beta_t \equiv \max(\beta_{t, { \boldsymbol{d}}}, \beta_{t, { \boldsymbol{b}}})$ stems from (2.4), applied to the two domains, respectively. Together with Lemma 1, this establishes the continuity of the left-hand side of (3.6). The $H^{1/2}$ integrability of the terms on $\Gamma_m$ ensure the bound in (3.9) by standard Sobolev embeddings (see e.g., ^[31]).

For the coercivity, notice that, since $\nabla\cdot \boldsymbol{u}_{ \boldsymbol{d}} = 0$ and $\boldsymbol{u}_{ \boldsymbol{d}} |_{\Gamma_M \cup \Gamma_{ \boldsymbol{d}, IN}} = 0$, we have

Assuming that on the outflow $\boldsymbol{u}_{ \boldsymbol{d}} \cdot \boldsymbol{n}_{ \boldsymbol{d}} \geq 0$, it follows that $(\boldsymbol{u}_{ \boldsymbol{d}}\cdot \nabla c_{ \boldsymbol{d}}, \, c_{ \boldsymbol{d}}) \geq 0.$ By standard arguments, the coercivity of $a_{ \boldsymbol{d}}(\cdot, \cdot)$ follows. Similar arguments apply to $a_{ \boldsymbol{b}}(\cdot, \cdot)$.

Collecting all these results, we have the following Theorem.

Theorem. The coupled problem given by $(3.6)$ is well-posed.

Proof. The bilinear form $\boldsymbol{A}(\cdot, \cdot)$ is continuous as an immediate consequence of the previous Lemmas. For the coercivity, note that:

By the Lax-Milgram theorem, this establishes the well-posedness of the problem (3.6).

3.2.   Numerical experiments

In this section we investigate the relationship between geometry, wall-shear stress, and fiber performance. We performed the simulations using the finite element library FreeFem++^[32]. In order to reduce the computational burden, we considered the three-dimensional problem as an axisymmetric problem. A mesh convergence analysis was performed to ensure simulation reliability and numerical stability. We solved the steady Navier-Stokes problem with Taylor-Hood $P^2$-$P^1$ elements and the scalar convection-diffusion problem with $P^1$ linear elements.

3.2.1.   Geometry and problem setup

We consider an axisymmetric problem in 3D. In the $r$ direction, the diameter of the blood domain in the hollow fiber is 0.5 mm (0.25 mm in the two-dimensional representation) and 0.5 mm for the dialysate domain. We consider a length $L$ of 20 mm in the $z$ direction. The dialysate domain is placed above the blood domain; in 3D, this corresponds to the blood domain being surrounded completely by the dialysate. The interface between the regions is defined as:

We refer to $\eta$ as the wave magnitude and $\omega$ as the wave frequency. The 2D geometry shown in the $(r, z)$ coordinates is shown in Figure 2 (left).

For the physical problem setup, we consider a countercurrent flow. We prescribe a flow profile vanishing at walls, resulting in a quadratic profile along the annulus, and a flow rate $Q_{ \boldsymbol{b}}$ of 0.1 mL/min in the blood domain and $Q_{ \boldsymbol{d}}$ = 0.2 mL/min in the dialysate domain, corresponding to physical flow conditions in a hollow fiber ^[26]. The blood viscosities are {0.03 g/cm$\cdot$s and 0.01 g/cm$\cdot$s for the blood} and dialysate respectively, with both fluids having a density of 1 gm/cm$^3$. The blood viscosity was assumed somewhat lower than normal due to the common use of blood thinners among dialysis patients ^[33]. $\mu_{ \boldsymbol{d}}, \, \mu_{ \boldsymbol{b}}$ are given by the diffusion coefficient of urea at 37 $^o$C in water and blood, respectively ^[33]. At the blood inlet, we consider a concentration of solute $c_{ \boldsymbol{b}, \, in}$ of 1 g/mL, with a corresponding value $c_{ \boldsymbol{d}, \, in}$ of 0 g/mL at the dialysate inlet. The efficacy of a given flow configuration is then evaluated with the clearance rate $Cl_s$, defined as: $Cl_s = Q_{ \boldsymbol{b}}\, c_{ \boldsymbol{d}, \, out}/c_{ \boldsymbol{b}, \, in}$. Additional numerical tests (not shown) showed that this quantity is insensitive to $c_{ \boldsymbol{b}, \, in}$. Our convergence analysis demonstrated grid-independence of the solution for the mesh-size adopted, and we did not have convection-dominated instabilities, so no stabilization was needed. Different (non axial-symmetric) geometries may suffer from instabilities and special numerical techniques will be considered ^[34].

3.2.2.   Wave configuration and clearance

We now investigate the relationship between the wave configuration and the clearance. We depict different values of wave magnitude ($\eta$) in Figure 2 (right). We test the geometry for the wave frequency of $\omega = 1, \, 2, \, 4, \, 8$ and then, for each $\omega$, a range of different $\eta$. We seek to identify what is the relationship between wave frequency, wave size, and clearance, and in turn, determine if an optimal configuration, in terms of maximal clearance, may be found.

For $\widetilde{\tau}$, we use 0.4 Pa, a value obtained using the Poiseulle law for wall shear stress when $\eta = 0$ in the given configuration, and is the median of the literature values shown in ^[17,18]. We first perform a sensitivity analysis for a range of physiological values of $\widetilde{\tau}$, comparing a straight cylinder design against a representative case of wavy design where $\omega = 2$ and $\eta = 0.03125$ mm. In Figure 4 (left), we plot the clearance results against the different threshold parameters. We see that, when the threshold is extremely low, neither the curved geometry nor the straight geometry achieves good clearance. For intermediate values, we achieve good clearance for the wavy geometry while the straight geometry fails to give clearance. Finally, for higher values of the threshold, we see that the difference in performance is negligible; hence, the curved geometry can only potentially improve performance. A similar sensitivity analysis was then performed on the flow rate, with analogous results (Figure 3). For a lower flow rate, the overall WSS is lower, and similar clearances are achieved for both curved and straight fibers. However, at higher flow rates, the straight fiber fails whereas the curved fiber continues to produce toxin clearance. We note that the baseline flow rates of 0.2 mL/min and 0.1 mL/min shown here are in line with current hemodialysis machines ^[26].

After confirming the positive effect of waves in general when compared to the straight geometry, we proceed to examine the impact on clearance of different wave configurations. We show the clearance for the different values of $\omega$ with respect to wave size $\eta$ in Figure 4 (right). In general, a large number of small waves results in a higher clearance compared to a small number of larger waves. While there is much variation on clearance based on $\eta$ within each value of $\omega$, we generally observe that higher $\omega$ is associated with more clearance. For a given $\omega$, we observe a positive impact of $\eta$ on clearance up to a certain level, at which the clearance reaches a maximum. After this point, further increases to $\eta$ begin to hamper clearance. This is expected; when examining (2.11), we see that increasing $\| \boldsymbol{\tau}\|$, as increases to $\eta$ do in general, will have a positive effect on transport until $\widetilde{\tau}$ is exceeded too much. The numerical results confirm the findings in ^[17,18] which show that balancing these effects is not straightforward, and provide compelling evidence that one may find optimal designs in terms of wave size and frequency. In summary, we feel the results demonstrate the potential for computer-aided design of such components.

4.   Conclusions

In the present work, we have introduced a first step towards the modelling of the cell-aided dialysis process. We have developed a model incorporating the effect of wall-shear stress on cell-layer toxin transport; namely, depending on the shear stress values, the shear stress can aid toxin transport through enhancing the epithelial function, or hinder it through reducing cell viability. We established the well-posedness of the model, then proceeded to numerical studies, in which our results show that a hollow fiber incorporating a wavy design may provide an effective way to counter the negative effects for flow configurations which may generate too much wall-shear stress for the cell-layer to perform adequately. In other situations in which the cell-layer is more resistant to stress effects, the waves positively impact transport. We further demonstrated, via numerical simulation, optimal design considerations in terms of number and size of waves. In general, we found many small waves perform better than a few large waves; however, for all given wave frequencies, an optimal wave size was found. This highlights the important role such simulation studies may play in the design of such devices.

This is a work in its infancy, and there are many fundamental steps necessary for future research. In future work, we may examine the wall-shear stress effects together with more advanced toxin transport models, such as Michaelis-Menten models. From the fluid dynamics end, the dependence of the viscosities and diffusivities on shear rate should be investigated. We consider only a single toxin in the present work. Our ultimate goal of detailed dialysis simulation, however, will require us to consider other uremic solutes and toxins, which may differ from patient-to-patient and do not display uniform physical properties. The presented model does not consider the mechanical properties of the membrane, and in particular the effects of the membrane thickness or stiffness on the resulting transport rate. Such properties may, however, influence the transport and should be considered in future work. In particular, we may expect membrane deformation to have an influence on the WSS profile, which would naturally influence the results of the presented model. We must also incorporate experimental results to both calibrate and validate the modelling framework discussed here. In summary, our results demonstrate the potential for the use of computational models in the study and optimization of advanced renal therapies.

Acknowledgments

This work is supported by the partners of Regenerative Medicine Crossing Borders (RegMed XB), a public-private partnership that uses regenerative medicine strategies to cure common chronic diseases. This collaboration project is financed by the Dutch Ministry of Economic Affairs by means of the PPP allowance made available by the Top Sector Life Sciences & Health to stimulate public-private partnerships. Alessandro Veneziani acknowledges the support of US NSF Project DMS 2012686. The authors would like to acknowledge Martina Moschella of Politecnico di Milano for help with graphics.

Conflict of interest

The authors declare no competing interests.