PDF
Electrochemical Treatment of Tumors
Introduction
The electrochemical treatment of tumors implies that diseased tissue is treated with direct current through the use of metallic electrodes inserted in the tumor. When tissue is electrolyzed, two competing reactions take place at the anode: oxygen evolution and chlorine production. The oxygen-evolution reaction also produces H+ ions, which lower the pH close to the anode. It should be stressed that chlorine production also leads to lowered pH through the hydrolysis of chlorine. One effect of low pH is the permanent destruction of hemoglobin in the tissue, which results in destruction of tumor tissue.
(1)
One challenge in developing this method of cancer treatment is in predicting the doses required for tumor destruction. One possible tool for dose planning is by modeling the reactions that take place close to the electrodes.
This example presents a first simple model for the development of dose-planning methods. More advanced models for dose planning, including secondary effects of chlorine, are found in Ref. 1, which also presents and solves models for the cathode.
Model Definition
This model uses the Tertiary Current Distribution, Nernst–Planck interface to predict the transport and reaction in the electrolysis of tumor tissue in a liver. A needle electrode is placed in the tumor, and transport is assumed to take place radially to and from this electrode. Because you can assume rotational symmetry, the computational domain reduces to a line (ra, rr) where ra is 1 mm and rr is 6 cm (see Figure 1).
The species you consider in the model are the protons, chloride, and sodium. At the surface of the anode you account for the chlorine and oxygen evolution reactions; see Equation 1.
Figure 1: Diagram of the cylindrical modeling domain inside a tumor.
This simplified model considers only a 1D model of the transport between two points, that is, between the two electrodes. The material balance for the species i is given by
where ci is the concentration (SI unit: mol/m3), Di give the diffusivities (SI unit: m2/s), zi equals the charge, umi represents the mobility (SI unit: (mol·m2)/(J·s)), and Ri is the production term for species i (SI unit: mol/(m3·s)), F denotes Faraday’s constant (SI unit: C/mol), and is the electrolyte potential (SI unit: V). The mobility, um i, can be expressed in terms of Di, R, and T as
The conservation of electric charge is obtained through the divergence of the current density:
At the electrode surface (r = ra) you use the Electrode Surface boundary node to specify the electrode reactions and the resulting fluxes for the ionic species that are included in the electrode reactions, H+ and Cl-. For the inert ionic species, Na+, the transport through the electrode surface equals zero. The expressions for molar fluxes at the boundary are based on the electrode reaction currents according to
where Ni is the flux, νij represents the stoichiometric coefficient for the ionic species i in reaction  j, nj is the number of electrons in reaction j, and jj is the current density of the reaction j.
You can express the current densities for the two reactions using the Electrode Reaction nodes. Introducing dimensionless pressure, P = p/pb, and concentration, C = c/cb, (where b denotes the reference value), the current density for the oxygen evolution is
where jI0, is the exchange current density (SI unit: A/m2) and ηI is the overpotential for the oxygen evolution reaction, defined as
where Eeq,I (SI unit: V) is the equilibrium potential for the oxygen evolution reaction.
Set the electrode potential to
where t denotes time.
The chlorine evolution reaction is similarly given by the expression
Using the input values nI = nII = 1, νH,I = -1, and νCl,II = 1, gives the fluxes at the electrode surface:
At the exterior boundary, assume the concentration is constant, ci = ci0, and ground the electrolyte potential.
The initial concentration is constant: ci = ci0.
Results and Discussion
The plot in Figure 2 shows the pH for different time steps. You can see that values below pH 2 are reached somewhere between 1800 and 2400 s. A closer examination reveals that it occurs after 2000 s. At this pH, tumor destruction starts to occur very rapidly according to the experimental and theoretical findings in Ref. 1.
Figure 2: pH-profiles at different time steps during the treatment.
The corresponding H+ profile in Figure 3 shows that the concentration maximum is not at the anode surface.
Figure 3: Proton concentration in the domain at different time steps.
This result arises because the current density is not constant over time. At high current densities, large amounts of protons are produced and this front moves inward in the domain as the current density is lowered.
The corresponding plot for chloride (Figure 4) shows a continuous decrease of chloride concentration close to the anode surface. This in turn decreases the production of chlorine and increases oxygen evolution.
Figure 4: Chloride concentration at different time steps.
The current-density plot in Figure 5 shows that the total current decreases rapidly as the concentration overvoltage for chlorine formation increases, due to lowered chloride concentration at the anode surface. The potential is then increased, which results in an increase in total current through increased oxygen evolution.
Figure 5: Total current density and current density for the competing reactions on the anode surface. Oxygen evolution is the lowest graph.
Reference
1. E. Nilsson, Modelling of the Electrochemical Treatment of Tumors, PhD. thesis, Dept. Chemical Engineering and Technology, Royal Inst. of Technology, Stockholm, Sweden, 2000.
Application Library path: Electrochemistry_Module/Electrochemical_Engineering/tumor
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  1D Axisymmetric.
2
In the Select Physics tree, select Electrochemistry > Tertiary Current Distribution, Nernst–Planck > Tertiary, Electroneutrality (tcd).
3
Click Add.
4
In the Number of species text field, type 3.
5
In the Concentrations (mol/m³) table, enter the following settings:
Na
H
Cl
6
Click  Study.
7
In the Select Study tree, select General Studies > Time Dependent.
8
Click  Done.
Global Definitions
Start by loading the model parameters and variables from text files.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file tumor_parameters.txt.
Definitions
Variables 1
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file tumor_variables.txt.
Geometry 1
Create the geometry as a single interval.
Interval 1 (i1)
1
In the Model Builder window, under Component 1 (comp1) right-click Geometry 1 and choose Interval.
2
In the Settings window for Interval, locate the Interval section.
3
In the table, enter the following settings:
Coordinates (m)
r_a
r_ext
4
Click  Build Selected.
5
Click the  Zoom Extents button in the Graphics toolbar.
Tertiary Current Distribution, Nernst–Planck (tcd)
Now start with the physics, begin with the electrolyte settings.
Species Charges 1
1
In the Model Builder window, under Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck (tcd) click Species Charges 1.
2
In the Settings window for Species Charges, locate the Charge section.
3
In the zNa text field, type z_Na.
4
In the zH text field, type z_H.
5
In the zCl text field, type z_Cl.
Electrolyte 1
1
In the Model Builder window, click Electrolyte 1.
2
In the Settings window for Electrolyte, locate the Diffusion section.
3
In the DNa text field, type D_Na.
4
In the DH text field, type D_H.
5
In the DCl text field, type D_Cl.
Initial Values 1
1
In the Model Builder window, click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the H text field, type H0.
4
In the Cl text field, type Cl0.
Electrolyte Potential 1
Now ground the potential of the exterior end, and then define the concentration at the same location.
1
In the Physics toolbar, click  Boundaries and choose Electrolyte Potential.
2
Select Boundary 2 only.
Concentration 1
1
In the Physics toolbar, click  Boundaries and choose Concentration.
2
Select Boundary 2 only.
3
In the Settings window for Concentration, locate the Concentration section.
4
Select the Species H checkbox.
5
In the c0,H text field, type H0.
6
Select the Species Cl checkbox.
7
In the c0,Cl text field, type Cl0.
Electrode Surface 1
Now set up the anode, and the anode reactions.
1
In the Physics toolbar, click  Boundaries and choose Electrode Surface.
2
Select Boundary 1 only.
3
In the Settings window for Electrode Surface, locate the Electrode Phase Potential Condition section.
4
In the ϕs,ext text field, type E_cell.
Electrode Reaction 1
1
In the Model Builder window, click Electrode Reaction 1.
2
In the Settings window for Electrode Reaction, locate the Stoichiometric Coefficients section.
3
In the νH text field, type -1.
4
Locate the Equilibrium Potential section. In the Eeq,ref(T) text field, type E_eqI.
5
Click to expand the Reference Concentrations section. In the table, enter the following settings:
Electrolyte species
Reference concentrations (mol/m^3)
H
H0
6
Locate the Electrode Kinetics section. In the i0,ref(T) text field, type j_I0.
Electrode Reaction 2
1
In the Physics toolbar, click  Attributes and choose Electrode Reaction.
2
In the Settings window for Electrode Reaction, locate the Stoichiometric Coefficients section.
3
In the νCl text field, type 1.
4
Locate the Equilibrium Potential section. In the Eeq,ref(T) text field, type E_eqII.
5
Locate the Reference Concentrations section. In the table, enter the following settings:
Electrolyte species
Reference concentrations (mol/m^3)
Cl
Cl0
6
Locate the Electrode Kinetics section. In the i0,ref(T) text field, type j_II0.
Global Definitions
Default Model Inputs
Set up the temperature value used in the entire model.
1
In the Model Builder window, under Global Definitions click Default Model Inputs.
2
In the Settings window for Default Model Inputs, locate the Browse Model Inputs section.
3
In the tree, select General > Temperature (K) - minput.T.
4
Find the Expression for remaining selection subsection. In the Temperature text field, type T.
Mesh 1
Create a user defined mesh with a very fine resolution close to the anode.
1
In the Model Builder window, under Component 1 (comp1) click Mesh 1.
2
In the Settings window for Mesh, locate the Sequence Type section.
3
From the list, choose User-controlled mesh.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Extra fine.
Size 1
1
In the Model Builder window, right-click Edge 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundary 1 only.
5
Locate the Element Size section. Click the Custom button.
6
Locate the Element Size Parameters section.
7
Select the Maximum element size checkbox. In the associated text field, type 1e-8.
8
Select the Maximum element growth rate checkbox. In the associated text field, type 1.1.
9
Click  Build All.
The finalized mesh should now look as follows:
Study 1
Step 1: Time Dependent
Modify the default solver to store all steps taken by the solver, and solve the problem.
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type 0 3600.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
Store the actual steps taken by the solver to avoid interpolation issues in the stored solution.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, locate the General section.
4
From the Times to store list, choose Steps taken by solver.
5
In the Study toolbar, click  Compute.
Results
Concentration, Na, 2D (tcd)
Some useful plots are generated by default. Polish them to illustrate the results more clearly.
Surface 1
1
In the Model Builder window, expand the Results > Concentration, Na, 2D (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Amber.
Arrow Surface 1
1
In the Model Builder window, click Arrow Surface 1.
2
In the Settings window for Arrow Surface, locate the Coloring and Style section.
3
From the Color list, choose White.
4
In the Concentration, Na, 2D (tcd) toolbar, click  Plot.
Concentration, H (tcd)
1
In the Model Builder window, under Results click Concentration, H (tcd).
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Time selection list, choose Interpolated.
4
In the Times (s) text field, type range(0,600,3600).
Line Graph 1
1
In the Model Builder window, expand the Concentration, H (tcd) node, then click Line Graph 1.
2
In the Settings window for Line Graph, click to expand the Legends section.
3
Select the Show legends checkbox.
4
In the Concentration, H (tcd) toolbar, click  Plot.
Concentration, H, 2D (tcd)
1
In the Model Builder window, under Results click Concentration, H, 2D (tcd).
2
In the Settings window for 2D Plot Group, click to expand the Title section.
3
From the Title type list, choose Label.
4
Locate the Color Legend section. Select the Show units checkbox.
Surface 1
1
In the Model Builder window, expand the Concentration, H, 2D (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Amethyst.
Arrow Surface 1
1
In the Model Builder window, click Arrow Surface 1.
2
In the Settings window for Arrow Surface, locate the Arrow Positioning section.
3
Find the R grid points subsection. In the Points text field, type 20.
4
Find the Z grid points subsection. In the Points text field, type 20.
5
Locate the Coloring and Style section. From the Arrow length list, choose Logarithmic.
6
In the Range quotient text field, type 10000.
7
In the Concentration, H, 2D (tcd) toolbar, click  Plot.
Concentration, Cl (tcd)
1
In the Model Builder window, under Results click Concentration, Cl (tcd).
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Time selection list, choose Interpolated.
4
In the Times (s) text field, type range(0,600,3600).
5
Locate the Legend section. From the Position list, choose Lower right.
Line Graph 1
1
In the Model Builder window, expand the Concentration, Cl (tcd) node, then click Line Graph 1.
2
In the Settings window for Line Graph, locate the Legends section.
3
Select the Show legends checkbox.
4
In the Concentration, Cl (tcd) toolbar, click  Plot.
Surface 1
1
In the Model Builder window, expand the Results > Concentration, Cl, 2D (tcd) node, then click Surface 1.
2
In the Settings window for Surface, locate the Coloring and Style section.
3
From the Color table list, choose Lagoon.
Arrow Surface 1
1
In the Model Builder window, click Arrow Surface 1.
2
In the Settings window for Arrow Surface, locate the Arrow Positioning section.
3
Find the R grid points subsection. In the Points text field, type 13.
4
Find the Z grid points subsection. In the Points text field, type 13.
5
Locate the Coloring and Style section. From the Color list, choose White.
6
In the Concentration, Cl, 2D (tcd) toolbar, click  Plot.
pH
1
In the Model Builder window, right-click Concentration, H (tcd) and choose Duplicate.
2
In the Settings window for 1D Plot Group, type pH in the Label text field.
3
Locate the Legend section. From the Position list, choose Lower right.
Line Graph 1
1
In the Model Builder window, expand the pH node, then click Line Graph 1.
2
In the Settings window for Line Graph, locate the y-Axis Data section.
3
In the Expression text field, type pH.
4
In the pH toolbar, click  Plot.
Electrode reaction current densities
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Electrode reaction current densities in the Label text field.
Point Graph 1
1
Right-click Electrode reaction current densities and choose Point Graph.
2
Select Boundary 1 only.
3
In the Settings window for Point Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck > Electrode kinetics > tcd.iloc_er1 - Local current density - A/m².
4
Click to expand the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dash-dot.
5
Click to expand the Legends section. Select the Show legends checkbox.
6
From the Legends list, choose Manual.
7
In the table, enter the following settings:
Legends
Oxygen evolution
Point Graph 2
1
Right-click Point Graph 1 and choose Duplicate.
2
In the Settings window for Point Graph, locate the y-Axis Data section.
3
In the Expression text field, type tcd.iloc_er2.
4
Click to expand the Title section. From the Title type list, choose None.
5
Locate the Legends section. In the table, enter the following settings:
Legends
Chlorine evolution
Point Graph 3
1
Right-click Point Graph 2 and choose Duplicate.
2
In the Settings window for Point Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Tertiary Current Distribution, Nernst–Planck > Electrode kinetics > tcd.itot - Total interface current density - A/m².
3
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Solid.
4
Locate the Legends section. In the table, enter the following settings:
Legends
Total electrode current density
5
In the Electrode reaction current densities toolbar, click  Plot.