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.
- finite element method (FEM)
- Magnetic Hyperthermia
- modeling of the temperature field
- magnetite nanoparticles
The Magnetic Hyperthermia is one of the most promising therapies in the cancer treatment . 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 .
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 R1 and the healthy region has the shell thickness R2 − R1. At the center O of this geometric structure, a ferrofluid volume Vf was injected with the volumetric flow rate (ferrofluid infusion rate) Qv (μl/min) using a needle of a syringe with the radius ro. The ferrofluid flow within this geometry starts from the center O where is localized the injection site (IS).
In this analysis, MNP with different sizes were considered. A lognormal distribution was defined by the following distribution function :
R is the particle radius, ln[R0] 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 . The velocity of the magnetic particles within tissues was computed as a solution of the continuity equation in the spherical coordinates:
The radial velocity of the particles—the component of the velocity vector: is given by:
where the variable r defines the radial distance from the IS localized at the center of this geometry. The constant was computed considering the ferrofluid velocity at the tip of the needle as:
is the needle cross sectional area. The radial velocity of the particles depends on the volumetric flow rate Qv (μl/min) and the radial variable r:
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 :
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 P1(r):
and P2(r) in the healthy region:
was computed solving the Darcy’s equation (3) for the concentric tissues. On the external border of the geometry, the pressure is zero, P2(r = R2) = 0. The pressure developed in this geometric configuration depends significantly on the parameter Qv, ferrofluid viscosity μ and tissues characteristics (porosity εi and permeability Ki).
The expression describes the particle convection, the term - the particle diffusion and - 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 . The coefficients
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 vr of the particles, tissues characteristics (porosity εi and permeability Ki) 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)):
where the expressions mi and Ai are:
The collector efficiencyas 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) C2 = 0 on the external boundary of the geometry (r = R2);
(ii) Neumann boundary condition at the all inner interfaces:
(iii) at the injection site (IS), at the top of the needle (r = ro) the concentration has the particular expression C1 = Cmax.
The radial distribution of the MNP concentrations computed as solutions of Eq. (5) depends mainly on: (i) the ferrofluid infusion rate Qv; (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:
|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/m3)||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|
ρ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 :
2.2. The temperature field
The temperature field within malignant and healthy tissues is described by the solutions: Ti = Ti(r, t), (i = 1, 2) of the bioheat transfer equation (Pennes equation) in the living tissues :
with the following thermal characteristics (Table 2): ρi—the mass density, ci—specific heat capacity, ki—the thermal conductivity, ρb—mass density of the blood, cb—specific heat capacity of the blood, Tart—blood temperature, —metabolic heat production and Qi(r) (W/m3)—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/m3)||5180||1160||1060|
|Specific heat capacity (J/kg K)||670||3600||3600|
|Thermal conductivity (W/mK)||40||0.4692||0.512|
|Anisotropy constant K (kJ/m3)||9||–||–|
|Frequency f (kHz)||100–650||–||–|
|Magnetic field amplitude H (kA/m)||0–15||–||–|
As a result of the spatial distribution of the particles, the total volumetric heating rate Qi(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 :
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
where the volumetric heating rate released by one particle P[R] and susceptibility χ″ [R] depend strongly on the particle radius:
H0 is the intensity of the magnetic field, f—frequency of the magnetic field, μ0 = 4π ∙ 10−7 H/m the permeability, Ms is the saturation magnetization and kB = 1.38 ∙ 10−23 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:
τ0 is the average relaxation time, η carrier liquid viscosity, 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 :
with For i = 1, 2, the expressions are:
Eq. (8) was solved in Appendix 2 at the thermal equilibrium. The general solutions are obtained:
The expressions v1, v2, v3 and v4 are:
with the following notations:
The integration constants: c1, c2, c3, c4 were computed from the following boundary conditions.
(i) The temperature T1 is finite at the center (r→ 0) of the geometric structure (Figure 1). As a result, the constant c1is zero. Therefore, the temperatures: T1(r) and T2(r) are (Appendix 2):
(ii) Dirichlet boundary condition was considered on the external surface of the healthy tissue:
(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:
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: (H0 and f) verify the criterion of exposure safe and tolerable for the ablation of the whole tumor according with Hergt condition: H0 f < 5 ∙ 109Am−1 s−1 . 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 Qv 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 Qv influences significantly the particles velocities within tissues. Figure 2(b) shows the radial dependence of the pressure for the same values of the parameter Qv. 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.
The radial distribution of the particles and temperature field were analyzed for different values of the parameter Qv in the range 10–40 μl/min. Figure 3(a) shows the evolution with distance from IS of the volume fraction of the particles Φ1(r) within tumor for different values of Qv. The value of the concentration at IS was Cmax = 10 mg/cm3. Consequently, the maxim value of the volume fraction at IS was Φmax = 1.93 ∙ 10−3. Variation of the parameter Qv 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 Qv during the ferrofluid infusion process.
The evolution with the parameter Qv of the deposition rate coefficient of the particles 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 decreases with the increase of the value of the parameter Qv. Also the deposition of the particles decreases with the increase of the tissue porosity. The particles with high velocities (high values of Qv) as a result of the high pressure gradient have no capability to remain deposited on the solid matrix.
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 Φ1(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.
The values of the magnetic field parameters: (H0 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 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.
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.
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,
The porous media contains the spherical collector cells with diameters dc ranged from 0.05 to 0.50 mm. The coefficients depend on the particle diameter D, collector diameter dc, porosities of the malignant and healthy tissues εi, and the radial velocity vr. The collector efficiencydescribes 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:
The single collector contact efficiency 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 :
The following non dimensional coefficients are defined :
—the aspect ratio;
—Peclet number which defines the ratio of the convective transport to diffusive transport;
—van der Waals number—the ratio of van der Waals interaction energy to the particle’s energy;
—attraction number—combined influence of van der Waals attraction forces and fluid velocity on particle deposition rate due to interception;
—gravity number—the ratio of Stokes particle settling velocity to approach velocity of the fluid;
; γ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 (i = 1, 2) represents the fractional reduction in deposition rate of the particles due to the presence of the electrostatic repulsive energies . Bai and Tien (1999), Chang and Chan  computed the expression of αi. The analytical expression was compared successfully with the experimental data. They found the following correlation equation:
In the absence of the electrostatic double layer forces, αi become 1. The nondimensional coefficients: NLO, NE1NE2 and NDL have the following expressions:
is London number;
is the first electrokinetic parameter;
is the second electrokinetic parameter;
NDL = 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 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 Ci = Ci(r) as solution of Eq. (5)
At equilibrium, in the steady-state: and Eq. (5) becomes:
where the deposition rate coefficients of the particles are given by the relations (A1.1):
with the following constants:
Solution of Eq. (A1.9) has the following form:
Considering Eq. (A1.10) can be written as:
which is equivalent with
The solutions of Eq. (A1.13) are given by the following expressions:
and are modified Bessel functions I and K of the order . The expressions are the variables of these functions. The general solutions of Eq. (5) are computed using the expressions (A1.14) in the expression (A1.11):
(const1)i and (const2)i are the four integration constants which are determined from the following four boundary conditions:
C2 = 0 on the external boundary of the geometry (r = R2);
Neumann boundary condition at the all inner interfaces;
at the injection site (IS), at the top of the needle (r = ro) the concentration has the particular expression C1 = Cmax.
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:
Using substitution Ri = rTi, these equations become:
These expressions contain the following notations:
The solutions of Eq. (A2.1) are:
c1, c2, c3, c4 are the integration constants and βi (i = 1, 2) have the expressions: