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.