Modeling of the Temperature Field in the Magnetic Hyperthermia

Abstract

The numerical and/or analytical modeling of the temperature field developed by the magnetic systems in the external alternating magnetic fields is essential in the Magnetic Hyperthermia. Optimization of the all parameters involved in the burning process of the malignant tissues can be realized more efficiently using a mathematical model. The analytical models can be used for the validation of any numerical complex models of the heating processes. This work focuses on the parameters which influences the therapeutic temperature field developed by the magnetic systems within the malignant tissues when the magnetic field is applied. An analytical model was developed to predict and control the bioheat transport within a malignant tissue. This model was compared with a numerical model which was developed in the same conditions of the thermal analysis. Infusion of a diluted suspension of magnetic nanoparticles (MNP) into liver tissue was modeled using the Darcy’s equation. The MNP concentration and the temperature field were computed for different parameters as: (i) ferrofluid infusion rates, (ii) particle zeta potential and (iii) magnetic field parameters. The convection-diffusion-deposition of the particles within tissues was considered in this analysis. This study indicates the essential role of these parameters to predict accurately the hyperthermic temperature field. The model presented in this paper predicts the optimum MNP dosage and the temperature at every point within the malignant tissue.

1. Introduction

The Magnetic Hyperthermia is one of the most promising therapies in the cancer treatment [1]. The malignant tissues are destroyed when their temperature reach a therapeutic hyperthermic temperature in the range 40–46°C [2, 3]. The main problem of this technique is to understand and to control as precisely as possible the temperature field developed by the magnetic systems injected within malignant tissues when the external alternating magnetic fields are applied.

Some experimental data realized in a tissue equivalent (agarose gel) evidences the particle diffusion within the tissue after their injection [4, 5]. The diffusion-convection and deposition of the particles have a strong influence on the radial particle distribution within the tissue volume and as a consequence on the temperature field in tissues [6].

In this paper, the temperature within a malignant tissue surrounded by a healthy tissue was studied considering the radial magnetic nanoparticles (MNP) concentration as an effect of the ferrofluid injection at the center of tumor. The MNP with different sizes having a lognormal particle size distribution were considered. The temperature developed by the magnetic systems in the external time-dependent magnetic field was analyzed for different values of the parameters. During the injection process of the particles within the tissues, their convection and deposition influences strongly the concentration of the particles. An analytical model was developed to predict the temperature field for different important parameters as: (i) ferrofluid infusion rates, (ii) particle zeta potential and (iii) optimal particle dosages. The results were compared with a numerical model in Comsol Multiphysics and Matlab. The values of temperatures computed using the analytical and numerical models in the same conditions were in very good agreement.

2. The analytical temperature model

Modeling the heat transport within the malignant tissues is an important element in the Magnetic Hyperthermia. The development of an analytical model for a complex system or process is usually an important breakthrough with a strong impact in the field.

In this case we are presenting an analytical model which will be a powerful analysis tool for hyperthermia treatment planning due to its ability to predict accurately the temperature field within the malignant tissues. The model will provide a tool to study the main parameters which influence significantly the temperature field and to give a tool to optimize them. For each patient, an individual therapy planning is required in correlation with the tumor location, geometry, shape and size. The elaboration of a patient model by segmentation of images from computerized tomography or magnetic resonance imaging scans is the first step—and the most important one—of hyperthermia treatment planning. This segmentation is used to compute the temperature field in the tumor located in a patient organ.

In our simulations a spherical concentric configuration composed of a malignant tissue surrounded by a healthy tissue was considered (Figure 1). The malignant tissue has a radius R₁ and the healthy region has the shell thickness R₂ − R₁. At the center O of this geometric structure, a ferrofluid volume V_f was injected with the volumetric flow rate (ferrofluid infusion rate) Q_v (μl/min) using a needle of a syringe with the radius r_o. The ferrofluid flow within this geometry starts from the center O where is localized the injection site (IS).

Figure 1.

The geometric configuration of a malignant and healthy tissue structure.

In this analysis, MNP with different sizes were considered. A lognormal distribution was defined by the following distribution function [7]:

$$g\left[\mathrm{R}\right]=\frac{1}{\mathit{\sigma R}\sqrt{2\pi }}exp\left[-\frac{{\left[ln\left(\raisebox{1ex}{$R$}\!\left/ \!\raisebox{-1ex}{$R_0$}\right.\right)\right]}^2}{2{\sigma }^2}\right]\text{and}{\int }_0^{\propto }g\left[\mathrm{R}\right]\mathit{dR}=1$$

R is the particle radius, ln[R₀] is the median and σ is the standard deviation of R.

2.1. The radial distribution of the MNP concentration

The ferrofluid (composed by small magnetite particles and water) was considered an incompressible diluted colloidal fluid with the small concentration (c ≤ 5% by volume). The presence of the small magnetite particles does not significantly affect the transport properties of the fluid [8]. The velocity of the magnetic particles within tissues was computed as a solution of the continuity equation in the spherical coordinates:

$$\nabla \bullet \overrightarrow{v}=0$$E1

The radial velocity of the particles—the component of the velocity vector: $\overrightarrow{v}\left(v_r,0,0\right)$ is given by:

$$v_r=\frac{B}{r^2}$$

where the variable r defines the radial distance from the IS localized at the center of this geometry. The constant $B=\frac{Q_v}{\pi }$ was computed considering the ferrofluid velocity at the tip of the needle as:

$$U=\frac{Q_v}{S_{\text{needle}}}=\frac{B}{r_0^2}(1.2)$$

$S_{\mathit{ne}\mathit{\text{edle}}}=\pi r_0^2$ is the needle cross sectional area. The radial velocity of the particles depends on the volumetric flow rate Q_v (μl/min) and the radial variable r:

$$v_r=\frac{Q_v}{\pi r^2}$$E2

Local velocity of the particles (2) depends on the pressure gradient developed within geometry as result of the ferrofluid injection process. The ferrofluid flow through tissues was modeled using the Darcy’s equation [8]:

$$\nabla P_i=-\frac{{\epsilon }_i\mu }{K_i}\overrightarrow{v}\left(i=1,2\right)$$E3

In this analysis, the index i = 1 defines the tumor (malignant tissue) and the index i = 2 defines the healthy tissue. The expression of the pressure in the malignant tissue is P₁(r):

$P_1\left(r\right)=\frac{{\epsilon }_1\mu }{K_1}\left[\frac{1}{r}-\frac{1}{R_1}\right]\frac{Q_v}{\pi }+\frac{{\epsilon }_2\mu }{K_2}\left[\frac{1}{R_1}-\frac{1}{R_2}\right]\frac{Q_v}{\pi },r_o\le r\le R_1.$E4.1

and P₂(r) in the healthy region:

$$P_2\left(r\right)=\frac{{\epsilon }_2\mu }{K_2}\left[\frac{1}{r}-\frac{1}{R_2}\right]\frac{Q_v}{\pi }{R}_1\le r\le R_2$$E4.2

was computed solving the Darcy’s equation (3) for the concentric tissues. On the external border of the geometry, the pressure is zero, P₂(r = R₂) = 0. The pressure developed in this geometric configuration depends significantly on the parameter Q_v, ferrofluid viscosity μ and tissues characteristics (porosity ε_i and permeability K_i).

The mass concentrations of the particles, C_i = C_i(r) (i = 1, 2) (expressed in mg/cm³) are the solutions of the modified convection-diffusion equation [8, 9]:

$$\frac{\partial C_i}{\partial t}+\nabla ·\left({\overrightarrow{u}}_i{C}_i\right)=\nabla ·\left(D_i^{\ast }\nabla C_i\right)-k_f^i{C}_i$$E5

where $k_f^i$ represent the values of the deposition rate coefficients of the particles within the malignant and healthy tissues which were computed using the relations (A1.3)–(A1.5) from Appendix 1.

The expression $\text{∇}\bullet \left({\overrightarrow{u}}_i{C}_i\right)$ describes the particle convection, the term $\text{∇}\bullet \left(D_i^{\ast }\text{∇}{C}_i\right)$ - the particle diffusion and $k_f^i{C}_i$ - the mean volumetric deposition rate of the particle on the solid phase.

The deposition processes play an important role in the spatial distribution of the particles. The temperature field within geometry is strongly dependent on the volumetric deposition rate of the particles on the solid phase $k_f^i{C}_i$. The coefficients $k_f^i$

$$k_f^i=\frac{3\left(1-{\epsilon }_i\right)}{2{\epsilon }_i{d}_c}{\eta }_s^i{v}_r,\left(\mathrm{i}=1,2\right)$$

were computed considering the superposition of the effects developed by the hydrodynamic forces, van der Waals interactions, gravity effect and the repulsive electrostatic double layer forces. The particle deposition on the cellular structure of tissues depends on the particle diameter D, the radial velocity v_r of the particles, tissues characteristics (porosity ε_i and permeability K_i) and particle zeta potential ζ_p(Appendix 1) [10, 11, 12, 13, 14, 15, 16]. At equilibrium, Eq. (5) have the following solutions (Appendix 1 - relations (A1.15)):

$$C_i\left(r\right)=\left(\frac{{\mathrm{e}}^{-\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}}}{\sqrt{\mathrm{r}}}\right)\left\{{\left(\mathbf{\text{const}}1\right)}_{\mathrm{i}}\text{Bessel}\mathrm{I}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}},\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]+{\left(\mathbf{\text{const}}2\right)}_{\mathrm{i}}\text{Bessel}\mathrm{K}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}},\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]\right\}$$E6

where the expressions m_i and A_i are:

$$m_i=-\frac{3\left(1-{\epsilon }_i\right)}{2{\epsilon }_i{d}_c}{\eta }_s^i\frac{Q_v}{\pi D_i^{\ast }}\text{and}{A}_i=\frac{Q_v}{\pi D_i^{\ast }}.$$

The collector efficiency${\eta }_s^i$ as result of the electrostatic repulsive forces is given by the relation (A1.2) from Appendix 1.

Boundaries conditions: The constants, (const1)_i and (const2)_i, were computed using the following boundary conditions:

(i) C₂ = 0 on the external boundary of the geometry (r = R₂);

(ii) Neumann boundary condition at the all inner interfaces:

$${\mathrm{C}}_1\left(\mathrm{r}={\mathrm{R}}_1\right)={\mathrm{C}}_2\left(\mathrm{r}={\mathrm{R}}_1\right)$$

$$D_1^{\ast }{\left.\frac{\text{∂}{\mathrm{C}}_1}{\text{∂}\mathrm{r}}\right|}_{\mathrm{r}={\mathrm{R}}_1}=D_2^{\ast }{\left.\frac{\text{∂}{\mathrm{C}}_2}{\text{∂}\mathrm{r}}\right|}_{\mathrm{r}={\mathrm{R}}_1}$$

(iii) at the injection site (IS), at the top of the needle (r = r_o) the concentration has the particular expression C₁ = C_max.

The radial distribution of the MNP concentrations computed as solutions of Eq. (5) depends mainly on: (i) the ferrofluid infusion rate Q_v; (ii) tissue and particle characteristics as porosity, permeability and particle size (Table 1). The convection, diffusion and deposition of the particles influence strongly the spatial distribution of the particles after their injection within tissues. At equilibrium, the temperature field depends significantly on the spatial distribution of the particles. In this analysis, the volume fractions of the particles were considered:

$${\mathrm{\Phi }}_{\mathrm{i}}\left(r\right)=\frac{C_i\left(r\right)}{{\rho }_{\mathit{MNP}}}$$E7

Hamaker constant A (J)	3 ∙ 10^−21 – 4 ∙ 10^−20
Particle radius R (nm)	5 – 30
Collector diameter dc (mm)	0.05 – 0.50
Tissue porosity: εi(i = 1, 2); ε1—malignant tissue porosity; ε2—healthy tissue porosity	ε1 = 0.1–0.8; ε2 = 0.2
Permeability Ki(i = 1, 2); K1—malignant tissue permeability; K2—healthy tissue permeability	K1 = 10^−14; K2 = 5 ∙ 10^−13;
Water mass density ρW (kg/m^3)	1000
Boltzmann coefficient kb (J/K)	1.38 ∙ 10^−23
Absolute fluid viscosity μ (kg/(s m))	0.001
ζp—particle zeta potential (mV)	−10 to −50
ζc—collector zeta potential (mV)	−20

Table 1.

Parameters values used in simulations [9, 10].

ρ_MNP is the mass density of the magnetic particles. At the injection site - at the center of the geometry - the maximum value of the volume fraction was ${\text{Φ}}_{\text{max}}=\frac{{\text{C}}_{\text{max}}}{{\rho }_{MNP}}$:

2.2. The temperature field

The temperature field within malignant and healthy tissues is described by the solutions: T_i = T_i(r, t), (i = 1, 2) of the bioheat transfer equation (Pennes equation) in the living tissues [6]:

$${\rho }_i{c}_i\frac{\partial T_i}{\partial t}=\nabla \left[k_i\nabla T_i\right]+{\rho }_b{w}_b{c}_b\left[T_{\mathit{art}}-T_i\right]+Q_{\mathrm{met}}^i+Q_i\left(r\right),$$E8

with the following thermal characteristics (Table 2): ρ_i—the mass density, c_i—specific heat capacity, k_i—the thermal conductivity, ρ_b—mass density of the blood, c_b—specific heat capacity of the blood, T_art—blood temperature, $Q_{\mathrm{met}}^i$—metabolic heat production and Q_i(r) (W/m³)—power density (volumetric heating rate) dissipated by the magnetic particles within geometric configuration when the magnetic field is applied.

Thermal and magnetic characteristics	Magnetite	Tumor tissue	Healthy tissue
Mass density (kg/m^3)	5180	1160	1060
Specific heat capacity (J/kg K)	670	3600	3600
Thermal conductivity (W/mK)	40	0.4692	0.512
Magnetization (kA/m)	446	–	–
Anisotropy constant K (kJ/m^3)	9	–	–
Frequency f (kHz)	100–650	–	–
Magnetic field amplitude H (kA/m)	0–15	–	–

Table 2.

The thermal and magnetic characteristics [9, 10].

As a result of the spatial distribution of the particles, the total volumetric heating rate Q_i(r) depends on the radial dependent volume fractions of the particles Φ_i(r). When MNP with different sizes are injected within tissues, the volumetric heating rate of the ferrofluid is [7]:

$${\mathrm{Q}}_{\mathrm{i}}\left(\mathrm{r}\right)={\mathrm{\Phi }}_{\mathrm{i}}\left(r\right)\overline{P}$$E9

The size distribution of particles influences strongly the heating rate. For a magnetic system which contains the particles with different sizes, the volumetric heating rate is given by

$$\overline{P}={\int }_0^{\propto }P\left[R\right]g\left[R\right]\mathit{dR},$$

where the volumetric heating rate released by one particle P[R] and susceptibility χ^″ [R] depend strongly on the particle radius:

$$P\left[R\right]={\mu }_{0}f\pi {\chi }^{″}\left[R\right]H_0^2\text{and}{\chi }^{″}\left[R\right]=\frac{{\mu }_0{M}_s^{2}V\left[R\right]}{3k_{B}T}\frac{2\pi f\tau \left[R\right]}{1+{\left(2\mathit{\pi f}\tau \left[R\right]\right)}^2}$$E10

H₀ is the intensity of the magnetic field, f—frequency of the magnetic field, μ₀ = 4π ∙ 10⁻⁷ H/m the permeability, M_s is the saturation magnetization and k_B = 1.38 ∙ 10⁻²³ J/K is Boltzmann constant. The effective relaxation time contains the Brown relaxation time τ_B[R] and the Néel relaxation time τ_N[R] as function of the particle radius:

$$\tau \left[R\right]=\frac{{\tau }_N\left[R\right]{\tau }_B\left[R\right]}{{\tau }_N\left[R\right]+{\tau }_B\left[R\right]}$$

$${\tau }_B\left[R\right]=\frac{3\eta V_H\left[R\right]}{k_{B}T};{\tau }_N\left[R\right]={\tau }_0\frac{\sqrt{\pi }}{2}\frac{exp\left[\frac{KV\left[R\right]}{k_{B}T}\right]}{\sqrt{\frac{KV\left[R\right]}{k_{B}T}}}.$$

τ₀ is the average relaxation time, η carrier liquid viscosity, $V_H\left[\mathrm{R}\right]=V\left[R\right]{\left(1+\frac{\delta }{R}\right)}^3$ is the hidrodinamic volume of the particles, δ is the surfactant layer thickness and K is the anisotropy constant of the magnetic particles. Using the spherical symmetry, Eq. (8) can be written as [17]:

$$\frac{k_i}{r^2}\frac{\partial }{\partial r}\left[r^2\frac{\partial T_i}{\partial r}\right]+M_i\left[T\right]={\rho }_i{c}_i\frac{\partial T_i}{\partial t}$$E11

with $M_i\left[T\right]={\mathrm{\omega }}_{\mathrm{b}}^{\mathrm{i}}{\rho }_b{c}_b\left[T_e^i\left(\mathrm{r}\right)-T_i\right].$ For i = 1, 2, the expressions $T_e^i\left(\mathrm{r}\right)$ are:

$$T_e^1\left(\mathrm{r}\right)={\mathrm{T}}_{\mathrm{B}}^1+\frac{Q_1\left(r\right)}{{\mathrm{\omega }}_{\mathrm{b}}^1{c}_b{\rho }_b}\text{and}{T}_e^2\left(\mathrm{r}\right)={\mathrm{T}}_{\mathrm{B}}^2+\frac{Q_2\left(r\right)}{{{\omega }^2}_b{c}_b{\rho }_b}$$

where

$${\mathrm{T}}_{\mathrm{B}}^1={\mathrm{T}}_{\mathrm{b}}+\frac{Q_{\mathrm{met}}^1}{{\mathrm{\omega }}_{\mathrm{b}}^1{c}_b{\rho }_b}\text{and}{\mathrm{T}}_{\mathrm{B}}^2={\mathrm{T}}_{\mathrm{b}}+\frac{Q_{\mathrm{met}}^2}{{\mathrm{\omega }}_{\mathrm{b}}^2{c}_b{\rho }_b}.$$

Eq. (8) was solved in Appendix 2 at the thermal equilibrium. The general solutions are obtained:

$$T_1\left(\mathrm{r}\right)={\mathrm{c}}_1\frac{cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{c}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_1\frac{cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}(\text{A}2.3)$$

and

$$T_2\left(\mathrm{r}\right)={\mathrm{c}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{c}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}$$

The expressions v₁, v₂, v₃ and v₄ are:

$${\mathrm{v}}_1=\frac{1}{{\mathrm{\beta }}_1}\int \mathrm{r}\left[{\mathrm{a}}_1+{\mathrm{b}}_1{\mathrm{\Phi }}_1\left(r\right)\right]sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)\mathrm{dr}\text{and}{\mathrm{v}}_2=-\frac{1}{{\mathrm{\beta }}_1}\int \mathrm{r}\left[{\mathrm{a}}_1+{\mathrm{b}}_1{\mathrm{\Phi }}_1\left(r\right)\right]cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)\text{dr}$$

$${\mathrm{v}}_3=\frac{1}{{\mathrm{\beta }}_2}\int \mathrm{r}\left[{\mathrm{a}}_2+{\mathrm{b}}_2{\mathrm{\Phi }}_2\left(r\right)\right]sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)\mathrm{dr}\text{and}{\mathrm{v}}_4=-\frac{1}{{\mathrm{\beta }}_2}\int \mathrm{r}\left[{\mathrm{a}}_2+{\mathrm{b}}_2{\mathrm{\Phi }}_2\left(r\right)\right]cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)\mathrm{dr}$$

with the following notations:

$$a_1={\beta }_1^2{\mathrm{T}}_{\mathrm{B}}^1\text{and}{a}_2={\beta }_2^2{\mathrm{T}}_{\mathrm{B}}^2$$

$${\text{b}}_1=\frac{{\beta }_1^2\overline{P\mathrm{ }}}{{\text{ω}}_{\text{b}}^1\mathrm{ }{c}_b\mathrm{ }{\rho }_b}{\text{and b}}_2=\frac{{\beta }_2^2\overline{P\mathrm{ }}}{{\text{ω}}_{\text{b}}^2\mathrm{ }{c}_{b\mathrm{ }}{\rho }_b}$$

$${\beta }_1^2=\frac{{\mathrm{\omega }}_{\mathrm{b}}^1{c}_b{\rho }_b}{k_1};{\beta }_2^2=\frac{{\mathrm{\omega }}_{\mathrm{b}}^2{c}_b{\rho }_b}{k_2}$$

The integration constants: c₁, c₂, c₃, c₄ were computed from the following boundary conditions.

Boundary conditions:

(i) The temperature T₁ is finite at the center (r→ 0) of the geometric structure (Figure 1). As a result, the constant c₁is zero. Therefore, the temperatures: T₁(r) and T₂(r) are (Appendix 2):

$$\begin{array}{c}{T}_1\left(\mathrm{r}\right)={\mathrm{c}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_1\frac{cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}\\ T_2\left(\mathrm{r}\right)={\mathrm{c}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{c}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}\end{array}$$E12

(ii) Dirichlet boundary condition was considered on the external surface of the healthy tissue:

$$T_2\left[\mathrm{r}={\mathrm{R}}_2\right]={37}^{°}\mathrm{C}$$E13

(iii) The Newman boundary conditions are considered at malignant—healthy tissue interface. The heat flux coming from the malignant tissue is completely received by the healthy region. The continuity condition of the heat fluxes is imposed at tumor—healthy region interface:

$$k_1{\left[\frac{\partial T_1}{\partial r}\right]}_{R_1}=k_2{\left[\frac{\partial T_2}{\partial r}\right]}_{R_1}$$E14

$$T_1\left[r=R_1\right]=T_2\left[r=R_1\right]$$E15

3. Results and discussions

In this analysis, the magnetite system with sizes in the range (5–30) nm was considered. The temperature field within a liver tissue was computed for different values of the magnetic field parameters. The malignant and healthy tissues are two concentric domains having the diameter of 20 mm and 100 mm, respectively (Figure 1). The temperature values depend strongly on the particle size and magnetic field parameters (frequency and amplitude). The values of the magnetic field parameters: (H₀ and f) verify the criterion of exposure safe and tolerable for the ablation of the whole tumor according with Hergt condition: H₀ f < 5 ∙ 10⁹Am⁻¹ s⁻¹ [7]. Eqs. (1), (3), (5) and (8) were solved in a numerical model using the finite element method (FEM). Their numerical solutions were compared by the previous analytical solutions in the same mathematical conditions.

Figure 2(a) shows the radial dependence of the particle velocity for values of the ferrofluid infusion rate Q_v in the range of 5–30 μl/min. The velocity of the particles on radial direction, decreases with the distance from IS according with the relation (2). The parameter Q_v influences significantly the particles velocities within tissues. Figure 2(b) shows the radial dependence of the pressure for the same values of the parameter Q_v. In agreement with relations (4.1) and (4.2), the pressure decreases with the distance from IS. The pressure developed within geometry as result of ferrofluid infusion depends strongly on the Qv. Higher values of Qv determines higher values of pressure and faster movements of the particles within tissues. As result of the convection process, the particles move on larger distances or remain in the small vicinity of the IS (at small radial distances). The parameter Qv and implicitly the pressure generated within tissue by the ferrofluid infusion influences strongly the convection—diffusion—deposition processes of the particles and implicitly the spatial distribution of the particles.

Figure 2.

(a) The radial velocity of the magnetite particles within tissues; (b) the pressure on radial direction within malignant tissue.

The radial distribution of the particles and temperature field were analyzed for different values of the parameter Q_v in the range 10–40 μl/min. Figure 3(a) shows the evolution with distance from IS of the volume fraction of the particles Φ₁(r) within tumor for different values of Q_v. The value of the concentration at IS was C_max = 10 mg/cm³. Consequently, the maxim value of the volume fraction at IS was Φ_max = 1.93 ∙ 10⁻³. Variation of the parameter Q_v determines different spatial distributions of the particles within the malignant tissue which influences strongly the temperature field (Figure 3(b)). Practically, the temperature values within the malignant tissue can be controlled in the therapeutic range (42–46)°C by using an optimum value of Q_v during the ferrofluid infusion process.

Figure 3.

The influence of the parameter Qv on the radial dependent volume fraction of the particles (a) and the temperature field (b).

The evolution with the parameter Q_v of the deposition rate coefficient of the particles $k_f^1$ was studied for different porosities of the malignant tissue, in order to understand the influence of this parameter on the deposition process of the particles on the solid porous matrix (Figure 4). The coefficient $k_f^1$ decreases with the increase of the value of the parameter Q_v. Also the deposition of the particles decreases with the increase of the tissue porosity. The particles with high velocities (high values of Q_v) as a result of the high pressure gradient have no capability to remain deposited on the solid matrix.

Figure 4.

The evolution with the parameter Qv of deposition rate coefficient of the particles, $k_f^1$ within the malignant tissue.

The repulsive electrostatic double layer forces influences the particle deposition process and implicitly the spatial distribution of the particles. Temperature field depends strongly on the particle zeta potential ζ_p. Figure 5(a) shows the evolution with radial distance from IS of the volume fraction of the particles Φ₁(r) for different values of the particle zeta potential ζ_p =  − 10 to − 40 mV. As a result of the strong repulsive electrostatic double layer forces, a number of the particles are deposited in the solid structure of tissue. This effect influences significantly the spatial distribution of the temperature (on radial direction) as Figure 5(b) shows. As a consequence, the temperature gradients become smaller in the case of smaller repulsive electrostatic double layer forces. The repulsive electrostatic interactions (as a result of the repulsive electrostatic double layer (EDL) forces) influence strongly the mass concentration of the particles and the spatial temperature field. The particle zeta potential ζ_p can be controlled in the ferrofluid (as liquid medium) due to the ionic conditions measured by pH and ionic strength.

Figure 5.

The influence of the parameter—zeta potential ξ_p on the radial dependent volume fraction of the particles (a) and the temperature field (b).

The values of the magnetic field parameters: (H₀ and f) are essential in the optimization of the Magnetic Hyperthermia therapy. In the following, the increase of the temperature on the radial direction with the frequency of the magnetic field was followed.

Figure 6(a) and (b) show a 3D and 2D view of the radial temperature field within the malignant tissue for different values of the frequency of the magnetic field. It was considered a small value of both Qv and the tissue porosity in order to analyze a temperature field with strong non-uniformity (and implicitly high thermal gradients).

Figure 6.

The evolution with the frequency of the magnetic field of the temperature field on radial direction Q_v = 10 μl/min; ɛ₁ = 0.2 and ξ_p = −30 mV. (a) 3D view and (b) 2D view.

Figure 7(a) shows the values of the main parameters Qv and f which determines the same temperature on the radial direction. Figure 7(b) shows the isothermal surfaces for different values of values of the main parameters Qv and f.

Figure 7.

The isothermal surfaces for different parameters. (a) Isothermal surfaces for different values of the parameter Qv and frequency of the magnetic field f and (b) isothermal surfaces for different values of the tissue porosity and frequency of the magnetic field f.

The analytical temperature was compared with the numerical results given by the finite element method (FEM) in Comsol Multiphysics using the same parameters in the same mathematical conditions. Good agreements were found between the predictions of the analytical model and numerical results.

The simulations allow (i) the optimization of the main parameters which influences strongly the heating of the tumor in the therapeutic temperature range and (ii) provide better temperature control through treatment preplanning.

4. Conclusion

The model developed in this paper analyzes the essential role of the ferrofluid infusion rate in the radial MNP distribution after their injection within a malignant tissue. Analytical correlations between the following parameters: (i) the particle velocity, (ii) the pressure developed in geometry during the ferrofluid infusion and (iii) the particle concentration were done in order to understand and predicts the temperature field within tissues when an external magnetic field is applied. The temperature field is concentrated within the malignant tissue. The temperature on the tumor border (approximately 38–39°C) not affects the healthy tissue.

The thermal energy deposited within the malignant tissue provides from the MNP distributed as result of convection-diffusion-deposition of the particles after their injection inside tissue. The ferrofluid infusion rate influences significantly the radial distribution of the particles and consequently the temperature field.

The temperature field within the malignant tissues can be controlled by the control of the ferrofluid infusion rate Qv during the infusion process. The particles having higher velocity moves on larger distances on radial direction from the injection site within tumor. As a result, the particles which not remain in the vicinity of the injection site are distributed in the tumor volume. This important effect determines a temperature field with small temperature gradients. The model developed in this paper can be used as a planning tool to compute the temperature field for different parameters.

2.1. Convection-Diffusion-Deposition of the particles

1) The computation of the deposition rate coefficients of the particles, $k_f^i$

The mean deposition rate coefficients of the particles are expressed by the following relation [10, 11, 12, 13, 14, 15, 16]:

$$k_f^i=\frac{3\left(1-{\epsilon }_i\right)}{2{\epsilon }_i{d}_c}{\eta }_s^i{v}_r,\left(\mathrm{i}=1,2\right)$$EA1.1

The porous media contains the spherical collector cells with diameters d_c ranged from 0.05 to 0.50 mm. The coefficients $k_f^i$ depend on the particle diameter D, collector diameter d_c, porosities of the malignant and healthy tissues ε_i, and the radial velocity v_r. The collector efficiency${\eta }_s^i$ describes the ratio of the particles captured by the solid surface to those brought into a unit structural cell of the porous medium [8, 13]. This coefficient is given by the expression:

$${\eta }_s^i={\alpha }^i{\eta }_0^i.$$EA1.2

The single collector contact efficiency ${\eta }_0^i$ is the fraction of the particles brought to the collector surface by the Brownian diffusion, interception and/or gravitational sedimentation. This coefficient was computed by N. Tufenkji and M. Elimelech considering the superposition of the effects developed by the hydrodynamic forces, van der Waals interactions and gravity effect [13]:

$${\eta }_0^i=\underset{\begin{array}{c}\mathit{\text{transport}}\mathit{by}\\ \mathit{\text{diffusion}}\end{array}}{\underbrace{{\eta }_D^i}}+\underset{\begin{array}{c}\mathit{\text{transport}}\mathit{by}\\ \mathit{\text{interception}}\end{array}}{\underbrace{{\eta }_I^i}}+\underset{\begin{array}{c}\mathit{\text{transport}}\mathit{by}\\ \mathit{\text{gravitation}}\end{array}}{\underbrace{{\eta }_G}}$$

with:

$${\eta }_D^i=2.4{\left(A_s^i\right)}^{1/3}{N}_R^{-0.081}{\left(N_{\mathit{Pe}}\right)}^{-0.715}{N}_{\mathit{vdW}}^{0.052};{\eta }_I^i=0.55A_s^i{N}_R^{1.675}{\left(N_A\right)}^{0.125};$$

$${\eta }_G=0.22N_R^{-0.24}{\left(N_G\right)}^{1.11}{N}_{\mathit{vdW}}^{0.053}.$$

The following non dimensional coefficients are defined [13]:

1. $N_R=\frac{D}{d_c}$—the aspect ratio;

2. $N_{\mathit{Pe}}=\frac{U{d}_c}{D_w}$—Peclet number which defines the ratio of the convective transport to diffusive transport;

3. $N_{\mathit{vdw}}=\frac{A}{\mathit{kT}}$—van der Waals number—the ratio of van der Waals interaction energy to the particle’s energy;

4. $N_A=\frac{A}{12\mathit{\pi \mu }{R}^{2}U}$—attraction number—combined influence of van der Waals attraction forces and fluid velocity on particle deposition rate due to interception;

5. $N_G=\frac{2}{9}\frac{R^2\left({\rho }_{\mathit{MNP}}-{\rho }_{\mathit{\text{ferrof}}}\right)g}{\mathit{\mu }\mathit{U}}$—gravity number—the ratio of Stokes particle settling velocity to approach velocity of the fluid;

6. $A_s^i=\frac{2\left(1-{\gamma }_i^5\right)}{2-3{\gamma }_i+3{\gamma }_i^5-2{\gamma }_i^6}$; γ_i = (1 − ε_i)^1/3—porosity dependent parameter of Happel’s model.

Deposition of the particles on the pore wall is influenced by the electrostatic repulsive forces. The attachment (collision) efficiency coefficient (filter coefficient) αⁱ (i = 1, 2) represents the fractional reduction in deposition rate of the particles due to the presence of the electrostatic repulsive energies [15]. Bai and Tien (1999), Chang and Chan [16] computed the expression of αⁱ. The analytical expression was compared successfully with the experimental data. They found the following correlation equation:

$${\alpha }^i=exp\left[\frac{1}{2}\left(ln\left({\alpha }_{C-C}^i\right)+ln\left({\alpha }_{B-T}^i\right)\right)\right],\left(\mathrm{i}=1,2\right).$$EA1.3

The expression of ${\alpha }_{C-C}^i$ is given by Chang and Chan [15, 16]:

$$\begin{array}{c}{\alpha }_{C-C}^i=0.024{\left(N_{\mathit{DL}}\right)}^{0.969}{\left(N_{E1}\right)}^{-0.423}{\left(N_{E2}\right)}^{2.880}{\left(N_{\mathit{LO}}\right)}^{1.5}+3.176{\left(A_s^i\right)}^{1/3}{\left(N_R\right)}^{-0.081}{\left(N_{\mathit{Pe}}^i\right)}^{-0.715}{\left(N_{\mathit{LO}}\right)}^{2.687}\mathrm{}\\ +0.222A_s^i{\left(N_R\right)}^{3.041}{\left(N_{\mathit{Pe}}\right)}^{-0.514}{\left(N_{\mathit{LO}}\right)}^{0.125}+{\left(N_R\right)}^{-0.24}{\left(N_G\right)}^{1.11}\left(N_{\mathit{LO}}\right)\end{array}$$EA1.4

The expression of ${\alpha }_{B-T}^i$ is given by Bai and Tien [15, 16]:

$${\alpha }_{B-T}^i=2.527\bullet {10}^{-3}{\left(N_{\mathit{LO}}\right)}^{0.7031}{\left(N_{E1}\right)}^{-0.3121}{\left(N_{E2}\right)}^{3.5111}{\left(N_{\mathit{DL}}\right)}^{1.352}$$EA1.5

In the absence of the electrostatic double layer forces, αⁱ become 1. The nondimensional coefficients: N_LO, N_E1N_E2 and N_DL have the following expressions:

1. $N_{\mathit{LO}}=\frac{4A}{36\pi \mu R^{2}U}$ is London number;

2. $N_{E1}=\frac{{\epsilon }_r{\epsilon }_0\left({\zeta }_p^2+{\zeta }_c^2\right)}{6\mathrm{\pi }\mathrm{\mu }\mathrm{R}U}$ is the first electrokinetic parameter;

3. $N_{E2}=\frac{2{\zeta }_p{\zeta }_c}{{\zeta }_p^2+{\zeta }_c^2}$ is the second electrokinetic parameter;

4. N_DL = 2 κ R is the double layer force parameter;

κ is Debye length for the colloidal suspension; ζ_p is particle zeta potential; ζ_c is collector zeta potential (Table 1) and $U=\frac{Q_v}{S_{\mathit{\text{needle}}}}$ is the ferrofluid velocity at the top of the needle.

The repulsive electrostatic double layer (EDL) forces appear in the liquid medium due to the ionic conditions measured by pH and ionic strength.

2. The computation of the MNP concentrations C_i = C_i(r) as solution of Eq. (5)

At equilibrium, in the steady-state: $\frac{\text{∂}{C}_i}{\text{∂}t}=0$ and Eq. (5) becomes:

$$\nabla ·\left(\overrightarrow{v}{C}_i\right)=\nabla ·\left(D_i^{\ast }\nabla C_i\right)-k_f^i{C}_i,$$EA1.6

where the deposition rate coefficients of the particles $k_f^i$ are given by the relations (A1.1):

$$\frac{1}{{\mathrm{r}}^2}\frac{\partial }{\partial \mathrm{r}}\left(v_r{r}^2{C}_i\right)-\frac{1}{{\mathrm{r}}^2}\frac{\partial }{\partial \mathrm{r}}\bullet \left(D_i^{\ast }{r}^2\frac{\partial C_i}{\partial r}\right)=-k_f^i{C}_i$$EA1.7

$$\frac{\partial }{\partial \mathrm{r}}\left(v_r{r}^2{C}_i-D_i^{\ast }{r}^2\frac{\partial C_i}{\partial r}\right)=M_i{C}_i$$EA1.8

$$\frac{\partial }{\partial \mathrm{r}}\left\{C_i{r}^2\left(\frac{A_i}{r^2}-\frac{1}{C_i}\frac{\partial C_i}{\partial r}\right)\right\}=m_i{C}_i$$EA1.9

with the following constants:

$$M_i=-\frac{3\left(1-{\epsilon }_i\right)}{2{\epsilon }_i{d}_c}{\eta }_s^{i}B;m_i=\frac{M_i}{D_i^{\ast }};A_i=\frac{B}{D_i^{\ast }};B=\frac{Q_v}{\pi }$$EA1.10

Solution of Eq. (A1.9) has the following form:

$$C_i\left(r\right)=C_0^i\left(r\right)exp\left[-\frac{A_i}{r}\right]$$EA1.11

Considering $C_0^i=C_0^i\left(r\right),$Eq. (A1.10) can be written as:

$$-\frac{\partial }{\partial \mathrm{r}}\left\{\left(r^{2}exp\left[-\frac{A_i}{r}\right]\right)\frac{\partial C_0^i}{\partial \mathrm{r}}\right\}=m_i{C}_0^{i}exp\left[-\frac{A_i}{r}\right]$$EA1.12

or:

$$r^2\left(\frac{{\partial }^2{C}_0^i}{\partial {\mathrm{r}}^2}\right)+\left(2r+A_i\right)\left(\frac{\partial C_0^i}{\partial \mathrm{r}}\right)+m_i{C}_0^i=0$$

which is equivalent with

$$\frac{\partial }{\partial \mathrm{r}}\left\{r^2\frac{\partial C_0^i}{\partial \mathrm{r}}\right\}+A_i\left(\frac{\partial C_0^i}{\partial \mathrm{r}}\right)+m_i{C}_0^i=0$$EA1.13

The solutions of Eq. (A1.13) are given by the following expressions:

$${\mathrm{C}}_0^{\mathrm{i}}\left(\mathrm{r}\right)=\frac{\left(\mathbf{\text{const}}\mathbf{1}\right){\mathrm{e}}^{\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}}\text{Bessel}\mathrm{I}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}} ,\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]}{\sqrt{\mathrm{r}}}+\frac{\left(\mathbf{\text{const}}\mathbf{2}\right){\mathrm{e}}^{\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}}\text{Bessel}\mathrm{K}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}} ,\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]}{\sqrt{\mathrm{r}}}$$EA1.14

$\text{Bessel}I\left[\sqrt{1-4m_i},\frac{A_i}{2r}\right]$ and $\text{Bessel}K\left[\sqrt{1-4m_i},\frac{A_i}{2r}\right]$ are modified Bessel functions I and K of the order $\sqrt{1-4m_i}$. The expressions $\frac{A_i}{2r}$ are the variables of these functions. The general solutions of Eq. (5) are computed using the expressions (A1.14) in the expression (A1.11):

$$C_i\left(r\right)=\left(\frac{{\mathrm{e}}^{-\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}}}{\sqrt{\mathrm{r}}}\right)\left({\left(\mathbf{\text{const}}\mathbf{1}\right)}_{\mathrm{i}}\text{Bessel}\mathrm{I}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}} ,\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]\mathrm{}+{\left(\mathbf{\text{const}}\mathbf{2}\right)}_{\mathrm{i}}\text{Bessel}\mathrm{K}\left[\sqrt{1-4{\mathrm{m}}_{\mathrm{i}}} ,\frac{{\mathrm{A}}_{\mathrm{i}}}{2\mathrm{r}}\right]\right)$$EA1.15

(const1)_i and (const2)_i are the four integration constants which are determined from the following four boundary conditions:

1. C₂ = 0 on the external boundary of the geometry (r = R₂);

2. Neumann boundary condition at the all inner interfaces;

$$C_1\left(\mathrm{r}={\mathrm{R}}_1\right)=C_{\mathit{2}}\left(\mathrm{r}={\mathrm{R}}_1\right)$$

$$D_1^{\ast }{\left.\frac{\text{∂}{\mathrm{C}}_1}{\text{∂}\mathrm{r}}\right|}_{\mathrm{r}={\mathrm{R}}_1}=D_2^{\ast }{\left.\frac{\text{∂}{\mathrm{C}}_2}{\text{∂}\mathrm{r}}\right|}_{\mathrm{r}={\mathrm{R}}_1}$$

1. at the injection site (IS), at the top of the needle (r = r_o) the concentration has the particular expression C₁ = C_max.

The constants (const1)_i and (const2)_i were computed in the Wolfram Mathematica 10 software.

2.2. The temperature model

At the thermal equilibrium, Eq. (7) are:

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left[r^2\frac{\partial T_i}{\partial r}\right]+{\beta }_i^2\left(T_e^i\left(\mathrm{r}\right)-T_i\right)=0\text{or}\frac{{\partial }^2{T}_i}{\partial r^2}+\frac{2}{r}\frac{\partial T_i}{\partial r}-{\beta }_i^2{T}_i+{\beta }_i^2{T}_e^i\left(\mathrm{r}\right)=0$$EA2.1

Using substitution R_i = rT_i, these equations become:

$$\frac{d^2{R}_i}{d{r}^2}-{\beta }_i^2{R}_i+{\beta }_i^2\mathrm{r}{T}_e^{\mathrm{i}}\left(\mathrm{r}\right)=0$$EA2.2

with solutions

$$R_1\left(\mathrm{r}\right)={\mathrm{c}}_{1}cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)+{\mathrm{c}}_{2}sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)+{\mathrm{v}}_{1}cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)+{\mathrm{v}}_{2}sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)$$

$$R_2\left(\mathrm{r}\right)={\mathrm{c}}_{3}cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)+{\mathrm{c}}_{4}sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)+{\mathrm{v}}_{3}cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)+{\mathrm{v}}_{4}sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)$$

where

$${\mathrm{v}}_1=\frac{1}{{\mathrm{\beta }}_1}\int \mathrm{r}\left[{\mathrm{a}}_1+{\mathrm{b}}_1{\mathrm{\Phi }}_1\left(r\right)\right]sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)\mathrm{dr}\text{and}{\mathrm{v}}_2=-\frac{1}{{\mathrm{\beta }}_1}\int \mathrm{r}\left[{\mathrm{a}}_1+{\mathrm{b}}_1{\mathrm{\Phi }}_1\left(r\right)\right]cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)\mathrm{dr}$$

$${\mathrm{v}}_3=\frac{1}{{\mathrm{\beta }}_2}\int \mathrm{r}\left[{\mathrm{a}}_2+{\mathrm{b}}_2{\mathrm{\Phi }}_2\left(r\right)\right]sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)\mathrm{dr}\text{and}{\mathrm{v}}_4=-\frac{1}{{\mathrm{\beta }}_2}\int \mathrm{r}\left[{\mathrm{a}}_2+{\mathrm{b}}_2{\mathrm{\Phi }}_2\left(r\right)\right]cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)\mathrm{dr}$$

These expressions contain the following notations:

$$a_1={\beta }_1^2{\mathrm{T}}_{\mathrm{B}}^1\text{and}{a}_2={\beta }_2^2{\mathrm{T}}_{\mathrm{B}}^2$$

$${\text{b}}_1=\frac{{\beta }_1^2\overline{P\mathrm{ }}}{{\text{ω}}_{\text{b}}^1\mathrm{ }{c}_b\mathrm{ }{\rho }_b}{\text{and b}}_2=\frac{{\beta }_2^2\overline{P\mathrm{ }}}{{\text{ω}}_{\text{b}}^2\mathrm{ }{c}_{b\mathrm{ }}{\rho }_b}$$

The solutions of Eq. (A2.1) are:

$$T_1\left(\mathrm{r}\right)={\mathrm{c}}_1\frac{cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{c}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_1\frac{cosh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}+{\mathrm{v}}_2\frac{sinh\left({\mathrm{\beta }}_1\mathrm{r}\right)}{r}$$EA2.3

and

$$T_2\left(\mathrm{r}\right)={\mathrm{c}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{c}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_3\frac{cosh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}+{\mathrm{v}}_4\frac{sinh\left({\mathrm{\beta }}_2\mathrm{r}\right)}{r}$$EA2.4

c₁, c₂, c₃, c₄ are the integration constants and β_i (i = 1, 2) have the expressions:

$${\beta }_1^2=\frac{{\mathrm{\omega }}_{\mathrm{b}}^1{c}_b{\rho }_b}{k_1};{\beta }_2^2=\frac{{\mathrm{\omega }}_{\mathrm{b}}^2{c}_b{\rho }_b}{k_2}$$