Inglês (pdf)
Artigo em XML
Referências do artigo
Tradução automática
Enviar este artigo por email
Citado por SciELO 
Citado por Google 
Similares em
SciELO 
Similares em Google 
Permalink
Introduction
Mathematical modeling techniques are a tool to provide solutions to complex systems in different fields of medicine and physiology (Bratus et al., 2017). In these fields, it is possible to obtain a mathematical model of any system using the principles of physics, chemistry, and biology (Selijteanu et al., 2015). However, these models represent the approximate physical problem (Akhmedova and Semenkin, 2013), a behavior that suggests the implementation of optimization algorithms that correlate mathematical models and experimental data in order to obtain a better understanding of the system (Rasdi et al., 2016).
A considerable number of researchers use optimization algorithms such as the Levenberg-Marquardt (LM) and Particle Swarm Optimization (PSO) methods to obtain information from systems that are difficult to measure (Chen et al., 2010; Pereyra et al., 2013). An optimization problem consists of minimizing or maximizing an objective function, with the purpose of finding the best available information to solve a problem. In the case of a complex problem with a unique solution, an optimization algorithm might find a local minimum that does not represent the feasible solution. The PSO method avoids these local minimums, allowing to find the global optimum of the objective function (Matajira-Rueda et al., 2018; Zhang, 2003).
Classical algorithms extract information sequentially, exploring the solution space in a unique direction (Chuang et al., 2012), whereas optimization algorithms can explore the solution space in different directions, increasing the probability of finding a feasible solution (Cornejo and Rebolledo, 2016).
Optimization algorithms do not need to detail the structure or behavior of the system: their function is to make random changes to the probable solutions, using an adjustment function in multiple variables in order to decide which solutions are optimal and computationally efficient. The implementation of optimization algorithms allows estimating the value of a variable in order to understand the behavior of a system's parameters under different conditions.
Radiofrequency (RF) hyperthermia therapy is a treatment for the partial or total elimination of cancer cells. It consists of increasing the tissue temperature between 40 and 45 °C via the induction of radiofrequency waves. This therapy is used as adjuvant therapy in traditional cancer treatments (Colombo et al., 2003; Curto, 2010; Horsman and Overgaard, 2007). Several researchers have conducted simulations of RF hyperthermia therapies in order to understand its behavior. These include the use of the finite element method (Gas, 2010; Gas and Miaskowski, 2015; Kurgan and Gas, 2009, 2010, 2011, 2015, 2016; Lv et al., 2005; Miaskowski et al., 2010; Miaskowski and Krawczyk, 2011; Miaskowski and Sawicki, 2013; Paruch and Turchan, 2018; Sawicki and Miaskowski, 2014; Yang et al., 2005), the boundary element method (Majchrzak, Drozdek, et al., 2008; Majchrzak, et al., 2008; Majchrzak and Paruch, 2009, 2010), and the finite-difference method (Deng and Liu, 2002). Other works have performed state variables and parameter estimation, which correlates numerical simulations and experimental data through algorithms based on Bayesian inference.
In the study by Bermeo et al. (2015)a, particle filters were used for the temperature field and heat source estimation. Lamien et al. (2017) and Bermeo et al. (2016a, 2016b) employed simultaneous parameters and state variables estimation while using the Liu and West filter, and Pacheco et al. (2020) the used Kalman filter to estimate temperature distribution. These analyses only consider a single frequency for the treatment. In the study by López et al. (2020), the estimation of electrical conductivity with three frequencies was performed via the LM method, with acceptable results in terms of accuracy.
Mathematical models in RF hyperthermia involve Maxwell's equations (Maxwell, 1865) and the bioheat transfer equation (Pennes, 1948). The modeling process poses a difficulty in the case of electrical conductivity in Maxwell's equations, as the values are different in intracellular and extracellular fluids, and, when the tissue is exposed to a uniform field, this behavior indicates that human tissue is inhomogeneous at a cellular scale (Ohmine et al., 2004; Peters et al., 2001). However, electrical conductivity can be assumed as constant in specific frequency ranges of the treatment because the differences are negligible, i.e., for each frequency value, there is a value of electrical conductivity. This assumption allows the evaluation of the model under certain conditions, considering that its environment is homogeneous (Gabriel, et al., 1996; Gabriel, Lau, et al., 1996; Haueisen et al., 1997; Kurup et al., 2012; Schepps and Foster, 1980).
Thus, in this work, electrical conductivity is estimated via the PSO algorithm, and the results are compared with the estimation process performed by means of the LM method. The inverse problem was performed via simulated temperature measurements in RF hyperthermia therapy with three different frequencies and powers. The results include the analysis of the sensitivity coefficients of the system and the statistical summary of the estimation process.
Particle Swarm Optimization method (PSO)
PSO is a stochastic method for the continuous optimization of parameters based on population behavior. The algorithm imitates the interactive behavior of the flocks of birds and banks of fish to reach an optimal solution (Kennedy and Eberhart, 1995). The method solves by improving a candidate solution, updating its velocity and position according to its individual experiences, as well as those of its neighbors (Chen et al., 2010). The PSO is a widely used tool for continuous optimization problems (Aazim et al., 2017; Alfi, 2011; Bermeo et al., 2015b; Muñoz et al., 2008; J. Tang et al., 2021).
In the PSO algorithm, the i-th particle is treated as a point within a space of N dimensions, and it is represented by a vector x . The best position found by the i-th particle is the position that produces the best value of the objective function, represented by Pbest i . ,j, while the best position found by the entire population is represented by Gbest.. The velocity for the i-th particle is represented by V . The velocity and position of the particles are calculated via Equations (1) and (2) (Kennedy and Eberhart, 1995):
where C 1 and C 2 are positive constants called the acceleration coefficients, and R 1 and R 2 are uniform random values for the interval [0;1] (Lashkari and Moattar, 2016). Equation (1) calculates the new velocity at which the particle moves in the solution area as a function of the current velocity and the position. w is the inertia coefficient that influences particle velocity, which affects the exploration capability of each individual towards the best location found by it (local best) and by the swarm in the search space (global best) (Kennedy and Eberhart, 1995). It allows particles to compensate their exploration ranges between particular and global in order to locate the optimal solution within a reasonable number of iterations. Algorithm 1 shows the canonical form of the conventional PSO algorithm, where S is the objective function, Np is the number of particles population, Pbest is the particle position, and Gbest is the best position of the particle.
According to Li et al. (2019), the implementation of an adaptive inertia coefficient (w) allows the particle population to adjust more precisely and quickly the exploration of the optimal solution, reducing the search times for particular and global values in the conventional implementation of the PSO algorithm. On the other hand, the variation in the acceleration coefficients is used to improve the speed of convergence towards the optimal global solution, thus allowing to reach this solution efficiently. Z. Tang and Zhang (2009) proposed a temporal variation of the acceleration coefficients in order to improve the exploration of the global optimum in the first stage of the optimization process and thus facilitate the convergence of the particles towards it. Equations (3)-(5) describe this proposal.
The modified algorithm is shown in Algorithm 2.
The Levenberg-Marquardt method (LM)
The LM method is used to iteratively solve nonlinear least-squares and linear problems for parameter estimation in highly complex inverse problems (Huang and Huang, 2007). This method combines the Gauss-Newton and gradient descent methods (Dattner and Gugushvili, 2018; Kanzow et al., 2005) controlled by a damping factor (¡j). If the damping factor tends to zero, the speed of convergence increases, and, if this factor tends to infinity, its convergence is slow and stepwise (Zhang, 2003). Equations (6) and (7) describe the parameter update and the objective function, respectively (Rouquette et al., 2007).
where P k+1 is the candidate parameter, J is the sensitivity matrix, Ω is the diagonal of the matrix JTJ S k+1 is the likelihood function, Y is the experimental temperature, and T pk is the numerical temperature. This method for estimating electrical conductivity is described in Algorithm 3.
Physical problem and mathematical formulation
The physical problem considered in this work involves a cylindrical domain with the physical properties of human muscle. The heating is generated using RF through a copper coil of three spirals with a radius of 8 mm and a cross-section of 70 mm, as presented in Figure 1. The heat in the domain is due to the electrical and magnetic loss of the therapy. The boundary conditions for the material are considered to be convective heat flux in the lower and upper surfaces and adiabatic boundary conditions around the domain. The solution of the direct (forward) problem was obtained using COMSOL Multiphysics® 5.3 and verified against the results presented by Hand et al. (1982) and Paruch and Turchan (2018). These processes were performed on a computer with 16 GB RAM and an i7-9750H processor.
The electromagnetic and thermal models of the domain were obtained via Maxwell's equations (Gratiy et al., 2017) and the bioheat transfer equation, respectively (Charny, 1992; Nakayama and Kuwahara, 2008; Pennes, 1948). Maxwell's equations are described in Equations (8)-(11).
where H is the magnetic field, J is the electric current density, B is the magnetic flux density, A is the magnetic vector potential, E is the electric field, σ is the electrical conductivity, and D is the electric displacement field.
The biological heat transfer is defined by Pennes (1948) as shown in Equation (12). Equation (13) presents the boundary conditions in the domain length L and the directions r, z.
where ρ is the density, cb is the specific heat capacity, k is the thermal conductivity, T0 is the initial temperature, ω b is the blood perfusion, Q m is the metabolic heat source, and Q is the heat source, as determined through Equations (14)-(16).
where Q rh is the electrical loss, Q ml is the magnetic loss, and Re is the real part of the losses (Lakhssassi et al., 2010).
The heat source inside of the copper coil to avoid overheating is calculated by means of Equation (17).
where M t is the mass flow, cpw is the heat capacity of water, Tin is the initial temperature, r is the inner radius, and A c is the cross section of the coil.
Results
In the estimation process of the electrical conductivity, the PSO and LM methods were evaluated using different frequencies and powers (Table 1 ). The simulated temperature measurements were considered in the inverse problem available at the center of the domain every 30 s. To avoid an inverse crime (Kaipio and Somersalo, 2004), the simulated temperature measurements were produced on a grid of 9 128 elements, while the inverse problem solution was performed on a grid of 3 091 elements. The initial conditions of the hyperthermia model were as follows: the temperature of the environment was T = 24°C, the initial temperature of the cooling pad was Tin = 20°C, the convection coefficient was h = 10 W/m2 K, the temperature in the domain was 24 °C (r = 0,04; -0,04 < z < 0,04), and the thermal properties of the biological tissue Q m , C b , and ω b , were zero because the domain was simulated as a solid material. The parameter values were assumed constant because the effects of temperature variation on the thermal and electrical properties are negligible (Rossmann and Haemmerich, 2014). This assumption is not possible at 10,0 MHz because the temperature increase is too high. Table 1 presents the electrical properties of the phantom with three different frequencies and powers, and Table 2 shows the physical properties of the domain.
Table 2 Physical properties of the domain
Source: Comsol Multiphysics (2012) and Hasgall et al. (2018)
Temperature field
Figures 2 to 4 show the simulated temperature measurements and numerical simulation at the center of the domain {0,0; 0,0} at 0,1, 1,0, and 10,0 MHz. These results indicate an increase in the temperature as the frequency increases. For 0,1 MHz, the temperature increase was 0,13 °C; 3,23 °C for 1,0 MHz; and 45 °C for 10,0 MHz. In terms of power, the temperature rises as the power increases.
To establish the appropriate frequency and power to perform RF hyperthermia, it is necessary to consider the thermal dose in the treatment. Experimental studies have revealed that patients are comfortable when the heating rate is 1,0 °C/min, which does not cause considerable pain or severe damage (de Oliveira, 2014). In this study, the mean rate at 0,1 MHz was 0,02 °C/min and 0,64 °C/min at 1,0 MHz. These rates are suitable for performing mild RF hyperthermia therapy.
For the case of 10,0 MHz, the rate was 1,8°C/min, which does not allow performing therapy. However, it is possible to modify the power to reach an adequate thermal dose, e.g., a frequency of 10,0 MHz at 100 W in approximately 9 minutes (López and Bermeo, 2021).
To establish the appropriate frequency and power to perform RF hyperthermia, it is necessary to consider the thermal dose in the treatment. Experimental studies have revealed that patients are comfortable when the heating rate is 1,0 °C/min, which does not cause considerable pain or severe damage (de Oliveira, 2014). In this study, the mean rate at 0,1 MHz was 0,02 °C/min and 0,64 °C/min at 1,0 MHz. These rates are suitable for performing mild RF hyperthermia therapy. For the case of 10,0 MHz, the rate was 1,8°C/min, which does not allow performing therapy. However, it is possible to modify the power to reach an adequate thermal dose, e.g., a frequency of 10,0 MHz at 100 W in approximately 9 minutes (López and Bermeo, 2021).
Analyzing the sensitivity coefficients
The sensitivity coefficients were analyzed in order to determine the parameters that influenced the temperature field in terms of dependence and magnitude (Özisik and Orlande, 2018). The coefficients were calculated using the finite-difference approximation method, as presented in Equation (18), which indicates the changes that occur in the temperature due to a low variation in the value of the parameters.
where J is the sensitivity matrix, i is the time instant, j is the number of the parameter, and ε is the perturbation of the parameter, represented at 10% of its nominal value.
Figure 5 shows the behavior of the sensitivity coefficients of k (thermal conductivity), ρ (density), σ (electrical conductivity), μ r (relative permeability), εr (relative permittivity), and cp (heat capacity). Note that those of μ r and £r have small magnitudes, indicating a low influence on the temperature field, and the remaining parameters exhibit linear dependence and large magnitudes. This behavior reveals that a single parameter can be estimated. As mentioned in this paper, electrical conductivity was estimated despite the limited information on this parameter in the literature.
Results of the estimation using the LM Method
The LM method was implemented with three different frequencies and powers (Table 1). The method was executed 50 times, and a statistical analysis was conducted in order to determine the electrical conductivity. The stopping criterion was likelihood. The study included the Shapiro-Wilk test (Shapiro and Wilk, 1965). Then, the mean, the standard deviation, and a confidence interval of 95% were calculated. Table 3 summarizes the estimation values. Note that the p-value in the Shapiro-Wilk test for the total cases is higher than 0,05, indicating that the results of 50 performances correspond to a normal distribution. In this sense, the estimation of electrical conductivity was represented by the mean, observing that the value was very close to the exact value, which indicates that the method accurately performed the estimation. Concerning the confidence interval, note that the exact and estimated values are within this interval.
Results of the estimation using the PSO Method
The PSO method was implemented using the same condition as the LM method. The algorithm was executed 50 times with a population of 50 particles for every computer simulation. The stopping criterion was the maximum number of iterations (50 by default). The study included the Shapiro-Wilk normality test, and the results are presented in Table 4. Note that the results of the 50 performances correspond to a normal distribution. Thus, the mean of a global best was very close to the exact value. The results of estimates are similar to those found by the LM method.
Discussions
Tables 3 and 4 show the statistical analysis for the methods implemented in RF hyperthermia therapy, demonstrating that the results obtained from the estimation are close to the exact value of electrical conductivity, with a difference of 0,1942%. Note that the maximum standard deviation of the total cases was 0,0130 S/m, indicating that the estimation values have good precision. This behavior was verified with the 95% confidence interval, which reveals that, within the range obtained, the exact value of the parameter lies within 95% certainty.
In the case of the LM method, 50 estimations were performed. Note that the mean values were very close to the estimated parameter. The estimated values were as follows: 0,3604, 0,5055, and 0,6175 S/m for 0,1, 1,0, and 10,0 MHz, respectively. These results indicate that the algorithm has a good convergence to the exact value. In terms of the computational time, the mean time of each estimation process was approximately 22 minutes, indicating that it is a fast-convergence algorithm. The PSO method was evaluated on a population of 50 particles and 50 iterations as a stopping criterion. Also, 50 estimation processes were performed, where the solution was obtained from the position of the best particle. The estimated values were 0,3614, 0,5033, and 0,6166 S/m for 0,1, 1,0, and 10,0 MHz, respectively. Note that this algorithm achieved an excellent convergence of the electrical conductivity value, with a difference of 0,0967%. In terms of the computational time, the mean time of each estimation process was approximately 327 minutes, which indicates that it is an acceptable convergence algorithm. This time could be lower if another stopping criterion were considered instead of the number of iterations, (e.g., the mean square error).
These results aim to predict the electrical conductivity in biological tissues using mathematical simulations and optimization algorithms on phantoms, where the final goal is to improve RF hyperthermia in terms of frequency and exposure time values, which may help medical doctors with the planning of individual treatment protocols. This also opens the possibility to design real-time control strategies.
Conclusions
This paper presents the application of two methods for estimating electrical conductivity in RF hyperthermia therapy for cancer treatment. The estimation problem was performed using the LM and PSO methods on a rectangular 2D axisymmetric system with temperature measurements available at one single location inside the domain. The estimation process was executed with three different frequencies and powers. The results were more accurate and smoother with the PSO method, and an excellent agreement with the exact values was obtained. In terms of computational time, the LM method performed better because its stopping criterion depends on the likelihood, whereas, for the PSO method, it is the number of iterations.
