Nazanin Maani1, Nitash Balsara2, Steven W. Hetts3, Vitaliy L. Rayz1
1Weldon School of Biomedical Engineering, Purdue University
2Department of Chemical and Biomolecular Engineering, University of California, Berkeley
3Radiology and Biomedical Imaging, University of California San Francisco

Abstract

A group of drugs used in Intra-Arterial Chemotherapy (IAC) have intrinsic ionic properties, which can be used for filtering excessive drugs from blood in order to reduce systemic toxicity. The ion-exchange mechanism is utilized in an endovascular Chemofilter device which can be deployed during the IAC for capturing ionic drugs after they have had their effect on the tumor. In this study, the concentrated solution theory is used to account for the effect of electrochemical forces on the drug transport and adsorption by introducing an effective diffusion coefficient in the advection-diffusion-reaction equation. Consequently, a multi-physics model coupling hemodynamic and electrochemical forces is developed and applied to simulations of the transport and binding of Doxorubicine in the Chemofilter device. A comparison of drug adsorption predicted by the computations to that measured in animal studies demonstrated the benefits of using the concentrated solution theory over the Nernst-Plank relations for modeling drug binding.

Introduction

Drug transport is one of the key aspects in the design of endovascular devices for drug filtration or elution. Drugs interact with blood and with medical devices through different mechanochemical processes, which depend on their properties. A large group of chemotherapeutic drugs have intrinsic ionic properties; therefore, the ion-exchange mechanisms can be utilized to enhance their effectiveness by reducing adverse side effects. Doxorubicin (Dox) is one of the most commonly used chemotherapeutic drugs used in the Intra-Arterial Chemotherapy (IAC) for treatment of solid tumors, e.g. the primary liver cancer. In the IAC procedure, a cocktail of chemotherapeutic drugs, including Dox, is injected into the arterial blood flow that feeds the tumor. Even though the IAC provides a targeted delivery, less than 50% of drug absorbs to the tumor, while the excessive drug remains in the circulation causing side effects such as an irreversible heart failure 1,2, 3. Dox is positively charged and, therefore, can be filtered from blood by binding to an ionic resin via ion-exchange mechanism in order to reduce systemic toxicity caused by the IAC. A catheter-based Chemofilter device was proposed for eliminating Dox from the venous flow after it has had its effect on the malignant cells of the tumor 1, 4-10. The Chemofilter deployed during the IAC procedure in veins draining the tumor would adsorb the drug to its surface coated with ionic resin, thus allowing for safer and more efficient treatment.
Significant advances in computational methods have contributed to development of mathematical models which elucidate the underlying mechanisms of convection, diffusion, and reaction. The motivation for this study was to develop a multi-physics modeling approach to analyze transport phenomena in concentrated and dilute electrochemical solutions, as a part of the Chemofilter design project. We present a new method for numerical modeling of drug transport and binding in an electrochemical system as well as the application of this method to simulation of transport and binding of Dox to the Chemofilter device. The modeling results are supported by comparison to measurements reported from animal studies.
A computational model capturing the transport of Dox in the flow and its binding to the ionic surface of the device has to couple the hemodynamic and electrochemical forces. In this study, the mathematical relationships are developed for electrochemical interaction of ionic particles with an ion-exchange surface resulting in a flux of the drug towards the surface. While traditionally the electrochemical force and chemical reaction are represented by a source term in the Navier-Stokes and advection-diffusion-reaction equations, herein the electrochemical body force is embedded in the diffusive term. Consequently, the passive diffusion coefficient is replaced by an effective diffusion coefficient, thus allowing to avoid the source term in the coupled transport equations. Electrochemical systems are typically modelled with the well-established Nernst-Plank equation; however, its application is limited to the dilute solutions. There are few studies of the concentrated solutions, such as that of Dox in blood. In this work, the contribution of electrochemical forces to the coupled Navier-Stokes and Advection-diffusion equations is modelled by using the concentrated solution theory.
In the Theory section of the paper, the relations for the effective diffusion coefficient are developed based on two alternative approaches: the dilute solution theory, and the concentrated solution theory. In the subsequent CFD Modeling section, the obtained relationships are used for numerical simulations of Dox transported through the Chemofilter device. The computational results obtained with each approach are compared to available experimental data.

Theory

The transport and binding of Dox in the Chemofilter device is affected by a complex interplay of hemodynamic and electrochemical forces. Electrochemical forces are dominant in the electric double layer (EDL), adjacent to the surface and defined as the region where binding of the particles is assumed to be instantaneous. In the EDL, the ion’s velocity field may be not divergence-free and is a function of the electrophoretic mobility of ions in the solution. The ion’s mobility is, in turn, a function of the electrostatic potential of the electrochemical system which is correlated to the migration of ions11. In this study, it was hypothesized that the diffusive and migration terms in the material balance equation could be combined by introducing an effective diffusion coefficient, as described in this section. In the dilute solution approximation, the induced migration of ions was accounted for by using the Nernst-Plank equations4, 11, 12, where the solution of Dox in plasma was approximated as 1) a binary solution, and 2) a non-binary solution. In the concentrated solution approximation, the closed form of the effective diffusion coefficient was derived for a binary solution of Dox in plasma. These alternative models for the effective diffusion coefficient are briefly described below.

Dilute Solution Theory (Nernst-Plank equation)

In an electrochemical system, the flux density of each dissolved species, \(N_{i}\) [mol/cm2], is due to three transport mechanisms – migration, diffusion, and advection:
\(N_{i}=-z_{i}u_{i}Fc_{i}\nabla\psi_{i}\ -D_{i}\nabla c_{i}+c_{i}\mathbf{v}\)(1)
where \(z_{i}\) is the number of proton charges carried by an ioni , \(u_{i}\) is the mobility of species i , \(F\) is the Faraday’s constant, \(c_{\text{i\ \ }}\)is the concentration of ioni , \(\psi\) is the electrostatic potential, \(D_{i}\) is the diffusion coefficient of species i , and \(\mathbf{v}\) is the advective velocity. The material balance for a minor component is formulated as:
\(\frac{\partial c_{i}}{\partial t}=-\ \nabla\bullet N_{i}+\mathcal{R}_{i}\)(2)
where \(\mathcal{R}_{i}\) is the reaction (source) term. In the binding of Dox to the ionic resin, the electrochemical force between\(\text{Dox}^{+}\) and \({\text{SO}_{3}}^{-}\)results in attraction of the particles towards the surface. Assuming a binary electrolyte, the Dox particles in the injected \(Dox-Cl\) solution dissociate (Eq. 3) to \(\text{Dox}^{+}\) cations in plasma (solvent) and bind to the surface (Eq. 4), where they form solid species which remain on the surface.
\(Dox-Cl\rightarrow\) \(\text{Dox}^{+}\ +\text{Cl}^{-}\)(dissociation in solution) (3)
\(\text{Dox}^{+}\ +{\text{SO}_{3}}^{-}\rightarrow Dox-\text{SO}_{3}\)(binding of Dox to the surface) (4)
By substituting the ion’s flux (Eq. 1) to the material balance (Eq. 2) and rearranging the terms, the balance equation for Dox cations reduces to:
\(\frac{\partial c}{\partial t}\mathbf{+v}\bullet\nabla c=D_{\text{eff}}\nabla^{2}c\mathcal{+R}\)(5)
where the concentration of the electrolyte is\(c=\ \frac{c_{+}}{\nu_{+}}=\frac{c_{-}}{\nu_{-}}\) and\(\nu_{+}\)and \(\nu_{-}\) are the number of moles of cations (+) and anions (-) that are produced from dissociation of one mole of the electrolyte (Appendix A). Eq. 5 is the advection-diffusion-reaction equation11, with the passive diffusion coefficient replaced by the effective diffusion coefficient:
\(D_{eff-db}=\ \frac{z_{+}u_{+}D_{-}\ -\ z_{-}u_{-}D_{+}}{z_{+}u_{+}\ -\ z_{-}u_{-}}\)(6)
where db stands for dilute-binary and the subscripts + and – stand for cation and anion, respectively. The migration term, the first term in Eq. 1, is present when an external potential is applied to the system. In the case of Dox binding to the ionic surface, no external potential is present. Instead, an induced potential term appears in the balance equation due to the electrophoretic mobility of ions. The mobility results in addition of a non-divergence-free term in the velocity, \(\mathbf{v}_{+}=\ -\ {z_{+}u}_{+}F\nabla\psi\), where\(\mathbf{v}=\ \mathbf{v}_{0}+\mathbf{v}_{+}\) (see details in Appendix). Therefore, a charge density is induced in the presence of a concentration gradient or a difference of the diffusion coefficients of the anion and cation (note that the diffusion coefficient of Dox is about one order of magnitude lower than that of the other ions such asCl- or Na+ ). This charge density creates a non-uniform potential which accelerates the ions with smaller diffusion coefficient towards the surface11.
For a non-binary electrolyte, the effective diffusion coefficient was derived for the case where the binding sites were occupied by other ions and the ionic bond on the surface should be overcome before the Dox particles could bind to the surface. This condition was investigated by Schlogl for an electrostatic system,13, 14 where the exchange of cations from the solution and the resin was considered as:
\({A_{s}}^{+}\ +{B_{r}}^{+}\rightarrow{A_{r}}^{+}+\ {B_{s}}^{+}\)(7)
where r denotes resin and s denotes solution. In Schloglet al., 13, 14 the Nernst-Einstein relationship was used to derive the flux of particles as a function of concentration gradient for an electrostatic system (Appendix A). In the current study, Schlogl’s model was expanded to model flow dynamics, which results in the material balance equation presented in Eq. 5 with an effective diffusion coefficient expressed as:
\(D_{eff-dnb}=\ -\frac{2D_{A}\ (C_{A}+C_{B})}{\left(\frac{D_{A}}{D_{B}}+1\right)C_{A}+{2C}_{B}}\)(8)
Where dnb stands for dilute-non-binary, and A and Bcorrespond to Dox+ andNa+ . In the new effective diffusion coefficient based on Schlogl binding conditions, \(D_{eff-dnb}\) is a variable which depends on the concentration of the ions.

Concentrated Solution Theory

Even though Nernst-Plank equation is widely used in the modeling of electrochemical systems, its application is limited to dilute solutions and its results are based on several other assumptions, such as considering the ions as point charges in order to utilize the Nernst-Einstein relationship 11. Therefore, in this study, we applied the concentrated solution theory to derive a more generally valid model of the electrochemical binding of Dox to the ionic resin 11, 15.
In general, there are three main parameters that define the performance of an electrolyte in a concentrated electrochemical system and can be found experimentally; namely the diffusion coefficient, D , the ionic conductivity, κ , and the transference number,\(\ t_{+}^{0}\) 15, 16. In the concentrated solution theory, the migration and diffusion is expressed in terms of electrochemical potential as shown in Eq. 9. This equation is analogous to Eq. 1 for dilute solution11.
\(c_{i}\nabla\mu_{i}=RT\ \sum_{j}{\frac{c_{i}c_{j}}{c_{T}\mathfrak{D}_{\mathbf{\text{ij}}}}(\mathbf{v}_{j}-\mathbf{v}_{i})}\)\(i,\ j=\ +,\ -,\ 0\) (9)
Where \(\mathbf{v}\) is the velocity of the species and its subscripts are corresponding to the cations (\(+)\), the anions (\(-)\),the solution (0), and \(c_{T}\) is the total concentration (\(\sum_{i}c_{i}\)). In this equation, the interaction of species\(i,\ j\) is expressed in terms of Stefan-Maxwell diffusion coefficients,\(\ \mathfrak{D}_{\text{ij}}\), to quantify the relationship between the species velocity, \(\mathbf{v}_{i}\), and the electrochemical potential gradient,\(\nabla\mu_{i}\)11.
For a binary electrolyte, the flux of the cation and the anion is expressed as:
\({\mathbf{N}_{+}=c}_{+}\mathbf{v}_{+}=-\frac{\nu_{+}\mathfrak{D}}{\nu\text{RT}}\frac{c_{T}}{c_{0}}\nabla\mu_{e}+\frac{\mathbf{i}\text{~{}t}_{+}^{0}}{z_{+}F}+c_{+}\mathbf{v}_{0}\)(10)
\({\mathbf{N}_{-}=c}_{-}\mathbf{v}_{-}=-\frac{\nu_{-}\mathfrak{D}}{\nu\text{RT}}\frac{c_{T}}{c_{0}}\nabla\mu_{e}+\frac{\mathbf{i}\text{~{}t}_{-}^{0}}{z_{-}F}+c_{-}\mathbf{v}_{0}\)(11)
where \(\nu=\nu_{+}+\nu_{-}\) and\(\mu_{e}=\nu_{+}\mu_{+}+\nu_{-}\mu_{-}\ \)and the current density is defined as:
\(\mathbf{i}=-\kappa\nabla\psi-\frac{\kappa}{F}\left(\frac{s_{i}}{n\nu_{+}}+\frac{\text{~{}t}_{+}^{0}}{z_{+}\nu_{+}}-\frac{s_{0}c}{n\text{sc}_{0}}\right)\nabla\mu_{e}\ \)(12)
where \(s_{i}\) is the stoichiometric coefficient of species iand n is defined as: \(s_{+}z_{+}+s_{-}z_{-}=-n\). It was assumed that the concentration gradient of the anion,\(\ \text{Cl}^{-}\), is negligible, since blood plasma already contains a high concentration of \(\text{Cl}^{-}\), and \(n=-1\). Moreover, the first term on the right- hand-side of Eq. 12 was neglected, as there was no external electric potential in the system; thus, the induced current is defined in terms of \(\nabla\mu_{e}\). The gradient of concentration is related to the gradient of electrochemical potential as described by Newman:11
\(\frac{\mathfrak{D}}{\nu\text{RT}}\frac{c_{T}}{c_{0}}\ c\ \nabla\mu_{e}=D\left(1-\frac{\text{dln}c_{0}}{\text{dlnc}}\right)\nabla c\)(13)
Rearranging the equations, as detailed in the appendix, the material balance equation was finally reduced to:
\(\frac{\partial c}{\partial t}+\nabla\bullet\left(c\mathbf{v}_{0}\right)=\nabla\bullet\left[D_{eff-cb}\nabla c\right]\mathcal{+R}\)(14)
where
\(D_{eff-cb}=\ D\left(1-\frac{\text{dln}c_{0}}{\text{dlnc}}\right)\left(1-\frac{\kappa}{z_{+}F}\left(\frac{s_{+}}{n\nu_{+}}+\frac{\text{~{}t}_{+}^{0}}{z_{+}\nu_{+}}-\frac{s_{0}c}{n\text{sc}_{0}}\right)\left(\frac{{\text{νRT}t}_{+}^{0}}{\mathfrak{D}c}\frac{c_{0}}{c_{T}}\right)\right)\)(15)
Eq. 15 provides a closed form expression for the effective diffusion coefficient, \(D_{eff-cb}\), which includes the effect of induced migration of ions in the concentrated solution, and is expressed as a function of the passive and Stefan-Maxwell diffusion coefficients (\(\text{D\ and\ }\mathfrak{D}\)), transference number (\(\text{~{}t}_{+}^{0}\)), conductivity (\(\kappa\)), and concentration of cation and solvent, which are all measurable in experiments.

CFD Modeling

The relationships derived above were used to conduct CFD simulations of Dox transport and binding in the Chemofilter for two alternative configurations of the device. The first configuration, named a Honeycomb Chemofilter (Fig. 1a), consists of parallel hexagonal channels arranged in three separate stages.9 Detailed CFD modeling of this three-staged design predicted its superior hemodynamic and drug adsorption performance.9 The other configuration, named a Strutted Chemofilter (Fig. 1b), was used in animal studies,4 and therefore CFD results for this configuration could be compared to the device filtration measured in vivo. The geometries for both configurations were generated in SolidWorks software and the CAD files were then imported to ANSYS ICEM for numerical mesh generation. The details of discretization, mesh sensitivity analysis, and numerical schemes used in the simulations were described in our previous publications.8, 9
Figure The configuration of (a) 3-stage twisted perforated honeycomb Chemofilter, and (b) single strutted Chemofilter
The coupled Navier-Stokes and Advection-Diffusion-Reaction equations were numerically solved in ANSYS Fluent, to calculate the flow and transport of Dox through the Chemofilter. The Chemofilter device is intended to be deployed in the hepatic veins, therefore the effect of the cardiac pulse can be neglected and the flow was modeled as steady. The flow was also modeled as laminar since the Reynolds numbers in these vessels are less than a hundred. The inlet velocity was set to 0.01 m/s to match the venous flow measured in porcine models, and the outflow boundary condition was assigned to the outlet of the model. In the adsorption of Dox to the Chemofilter, the chemical reaction on the surface was expressed as a sink term to model the elimination of the captured Dox from the system. During the IAC procedure, a steady dose of Dox is injected in the artery for about 10 minutes. In this time interval, it can be assumed that the percentage of Dox that intercalates in tumor cells, as well as the binding capacity of the Chemofilter, i.e. the number of available binding sites, remain constant. In other words, we assume that the mass fraction of Dox at the vein’s inlet is constant and the Chemofilter surface does not saturate. Therefore, the steady state conditions were assumed for modeling the Dox injection and transport.
To model Dox binding to the Chemofilter, the energy and species transport modules were activated in Fluent. The chemical reaction was modeled by solving the material balance equation for all species that were introduced in a mixture except the bulk fluid (blood). The chemical reactions were based on Arrhenius model (\(k=A_{r}\ {T\ }^{\beta_{r}}{e\ }^{{-E}_{r}/RT}\)), where \(k\) is the rate constant. A finite rate reaction was chosen for the flow and Arrhenius constants were \(A_{r}=1e20\), \(\beta_{r}=0\), and\(E_{r}=0\). Mass deposition source was activated to include the effect of surface mass transfer in the continuity equation. The density was calculated based on volume rated mixing law and the diffusivity of Dox in plasma (2.442x10-10 m2/s) was used in the simulations, since the majority of the Dox molecules reside in blood plasma.
Dox mass fraction at the inlet was set to 0.005, and the species site density (of the sulfonate group on the surface) was set to 10-8 kgmol/m2. It was assumed that the temperature did not change in this process, thus, the energy balance equation was not solved. The diffusion coefficients of chloride anion and sodium cation in blood were set to 2.032x10-9 , and 1.334x10-9 m2/s, respectively.11 In the simulations based on concentrated solution theory, the electrochemical properties of a high-molecular-weight non-structured polystyrene-b -poly(ethylene oxide) (SEO) copolymer electrolyte doped with a lithium salt was utilized16, due to the lack of experimental data for Dox performance in plasma. In the species transport module, the effective diffusion coefficient was incorporated into each numerical model with an external User Defined Function (UDF). The UDF was developed in C language and implemented via the species transport dialogue box in Fluent.

Results

In order to assess the influence of the diffusion coefficient on Chemofilter filtration, numerical simulations of Dox binding to the Honeycomb Chemofilter were conducted with three different diffusion coefficients (10-10, 10-9, and 10-8 m2/s). The simulations were then conducted with the estimated effective diffusion coefficient based on the binary concentrated solution approximation for both Strutted and Honeycomb configurations. Moreover, the simulations with the effective diffusion coefficient derived from dilute solution theory were performed to determine the difference between the results obtained using these two alternative theories. For the dilute solution approximation, the simulations were conducted for both binary and non-binary electrolytes. The mathematical relationships used in the above models are summarized in Table 1.
Table . Computational models considered for Dox binding to the Chemofilter surface

Model 1: Filtration based on Passive binding

Figure 2 shows the results of the first study (Table 1), where Dox transport was simulated with different values of the passive diffusivity of Dox particles in plasma. The concentration of Dox decreases as blood flows through each stage of the Honeycomb Chemofilter. The passive diffusion coefficient was set to 10-10, 10-9, and 10-8m2/s, corresponding to the Peclet number of 50000, 5000, and 500 for the average velocity of 0.01 m/s. The overall Dox mass fraction downstream of the device was predicted to be 0.00476, 0.00373, and 0.00167, corresponding to Dox concentration reduction of 4.7%, 25.4%, and 66.5% for the diffusion coefficients of 10-10, 10-9,10-8m2/s, respectively. The binding predicted in the simulation with diffusion coefficient of 10-8m2/s was the closest match with the binding measured in animal studies, giving the cue that the effective diffusion coefficient must be in the same order of magnitude as the coefficient used in this simulation.
Figure 2 The heat map of Dox concentration changes computed for the Honeycomb Chemofilter using a constant diffusion coefficient of a) 10-10 m2/s, b) 10-9 m2/s, and c) 10-8 m2/s. (In all cases, inlet velocity was 0.01 m/s and inlet Dox mass fraction was 0.005)

Model 2: Filtration based on Concentrated Solution Theory

Figures 3a and 3b show a qualitative comparison of the transport and binding of Dox for the Honeycomb and Strutted configurations of the Chemofilter, respectively. This model was based on the concentrated solution theory, which provides a general platform for modeling an electrochemical system. These results were obtained for the effective diffusion coefficient, Deff, based on the properties of the SEO polymer electrolyte16. The calculated effective diffusion coefficient for the SEO polymer electrolyte was about 65 times larger than the passive diffusion coefficient of Dox particles in plasma. The reduction of Dox concentration in blood for the Honeycomb Chemofilter was predicted to be 58.4%. The concentration reduction for the Strutted Chemofilter was predicted to be 43.28%. The computational model slightly underestimated the concentration reduction of 54.1±5%, that was measured in the animal studies. It should be noted that the Deff was calculated for SEO electrolyte and it is suspected that the effective diffusion coefficient of Dox in blood would be larger than that of SEO copolymer.