Abstract
In the paper, a multi-objective optimization of an electrode pair for electrochemotherapy is performed. The aim is to design the electrode pair geometry in order to obtain an effective treatment, which means a homogeneous and suitable electric field distribution. To solve this problem, two recently developed algorithms, M-NSGA and 𝜇-BiMO, are applied.
Introduction
Electrochemotherapy (ECT) is a cancer therapy currently used in clinical practice to treat skin tumors with a local approach. The principle of this therapy is based on the reversible cell membrane permeabilization that can be induced by pulses of electric field [1–3]. The permeabilization of cell membrane can improve the cell uptake of chemotherapeutic drugs delivery to cancer. This therapy is currently applied in clinical practice to melanoma, other skin tumors and chest wall recurrence of breast cancer [4,5]. In ECT, to treat large surface tumors, flexible electrodes, composed of an array of needles, could be studied for the application of the treatment to human tissues [6,7]. The positioning of the needles, their dimensions and their relative inclinations strongly affect the distribution of the electric field [8,9] as well as the tissue inhomogeneity and the applied voltage.
Since the real tumor tissue could show inhomogeneous composition like described in [8] where fibrous tissue could be interposed in fat tissue. It is well known that different tissues show different electrical conductivity [10,11], then the distribution of the electric field is different with respect to the homogenous tissue [8,12,13]. At experimental point of view in [8,12], the authors propose an experimental model that uses potato coupled to a gel for which different conductivity could be suitably designed. Moreover, electroporated potato become dark where the electric field strength overcome electroporation threshold [9]. This way the distribution of the electric field in a bi-material model in [8,12] is compared to the one in homogenous model evaluating the difference in terms of extension of the dark area.
The electric field distribution in an inhomogeneous tumor tissue could be studied using the Finite Element simulation [14–18]. It is well known that in electroporation the electric field intensity has to overcome a prescribed threshold in order to permeabilize cell membrane [1,19], nevertheless the inhomogeneity of tissue can modify the electric field strength in some treatment area [8]. This way Finite Element analysis is coupled to optimization algorithm in order to improve the therapy set-up. The aim of optimization is the modification of electric field intensity in order to guarantee the electroporation of the tissue in the critical area in the treatment volume. The proposed bi-material model, with potato and gel to simulate the tissue inhomogeneity, derives from experimental model in [8,12].
In the paper, in order to obtain a homogenous and effective ECT treatment, an electrode pair is optimized, by evaluating the electric field distribution with finite element models, considering non-uniform tissues.
The forward problem
The forward problem, used to simulate the inhomogeneity effect on the electric field distribution, is a two-materials 3D geometry where three volumes with different conductivity are considered. In particular, the computational model is a parallelepiped (5 × 3.5 × 1.5 cm) made of potato (homogeneous case) or potato with gel (non-homogeneous case); the gel is highlighted in Fig. 1a. Potato material has been chosen since it is used in experiments because it shows a dark colour after few hours from electroporation procedure [20]. Two cylindrical electrodes (Fig. 1a), with distance D (can vary between 7.3 and 25 mm), diameter d (can vary between 0.2 and 1.8 mm) are supplied by a voltage U (that can vary between 500 and 3000 V) [14–17]. The potato is characterized by a non-linear electrical conductivity σ (E) [21]:
(a) 3D electrode pair geometry for numerical computation, (b) conductivity versus electric field curve.
The electric field E is computed solving a static conduction problem by means of Finite Element Analysis. In particular, the Laplace equation in electric potential V, was used using a commercial software [17,23]:
The electric field distribution found by finite-element models in the homogeneous and non-homogeneous cases are reported in Fig. 2 [22].
Electric field distribution, (a) homogenous case, (b) and (c) non-homogeneous cases with gel A and B, respectively.
The design variables considered for the optimization are the distance between electrodes D, the electrode diameter d and the applied voltage U. Considering a sampling surface S in the xy plane at the midline of the electrodes with length L (200 × 200 sampling points, see Fig. 3), and the distribution of the electric field E along that surface, the objective functions are defined as follows:
Surface considered for electric field sampling with D distance between needles and – L length of the needle.
f 1 - number of points on S in which E > 900 V/cm, to be minimised;
f 2 - number of points on S in which E < 300 V/cm, to be minimised.
For solving this problem, two multiobjective algorithms recently developed are applied, Migration Non-dominated Sorting Genetic Algorithm M-NSGA [24,25] and micro-Biogeography-inspired Multi-Objective Optimization 𝜇-BiMO [26].
The Migration-NSGA (M-NSGA) optimization method is based on the classical NSGA-II multi-objective algorithm [27]. This multi–objective optimization method starts from an initial population of N i individual, where usually N i = 20. The M-NSGA version modify NSGA_II in order to allow the migration of a new population of N m individuals (with N m < N i ) inside the current population of Ni individuals. The idea is based on the concept that new individuals could modify the genetic heritage of the current population. This way an improved population could be found. The event of the migration occurs periodically during the genetic process. Nevertheless, the time of the migration event is modified by the algorithm based on the entity of front movement with respect to the previous step. In practice, if the position of the front is changed with respect to the previous step, the migration event is less frequent with respect the case in which the front doesn’t move a lot. This functionality is managed by the algorithm as described in [24,28,29].
The nature-inspired algorithm i.e. biogeography-inspired optimization algorithm (BiMO) was recently extended to the solution of multi-objective inverse problems [30,31] and to micro ecosystems [26] by the authors.
Each solution considered in a biogeography-inspired optimization algorithm is treated as a habitat (design vector) composed of suitability index variables (SIV, design variables), and each habitat exhibits a quality given by the habitat suitability index (HSI, objective function) [32]. In contrast to other algorithms e.g. the genetic ones, in BiMO a selection between the previous and the actual population is not done, but the ecosystem is progressively modified by means of two stochastic operators, i.e. migration and mutation. The migration operator improves the HSI of poor habitats by sharing features from good habitats; in turn, the mutation operator modifies some randomly selected SIV of a few habitats in view of a better exploration of the ecosystem (design space).
In particular, at each iteration, habitats are sorted from the best ones to the worst ones according to the relevant value of the generalized fitness value. For each SIV of each habitat a random event r j , rules the selection of immigration or emigration operator (Fig. 4).
Schemes of immigration and emigration events.
If immigration occurs, the current SIV is kept for the next ecosystem, while if emigration occurs the current SIV in the considered habitat (say the i-th habitat) goes extinct and the SIV of another habitat (say the k-th habitat) takes the same SIV position in the i-th habitat of the next ecosystem. The k-th habitat is selected depending on the emigration rate and the SIV.
In the paper, the multiobjective version of the algorithm, already presented by the authors in [30,31], is applied. It is based on the definition of the generalized fitness, which takes into account simultaneously two or more objective functions by exploiting the concept of non-dominated ranking of solutions in the objective space and crowding distance.
The BiMO algorithm showed good results in different field of applications [33,34]. However, it happens that, when the number of SIV and the number of islands is small, many duplicates of islands occur. In the original version of the algorithm, the duplicates are handled by randomly varying duplicated islands. This is like implementing a strong mutation on the islands and hence the exploration is improved while the exploitation is worsened.
To overcome this drawback, a new version of BiMO called 𝜇-BiMO is implemented, where the role of small rocks in the migration of individuals is taken into account. As in reality the small rocks help immigrants to colonize islands that otherwise would not be reached, with the concomitant loss of the individuals who would never reach the ground, in this method the rocks have the function not to waste habitats that otherwise would never characterize an ecosystem.
In particular, during the migration procedure it could happen that good habitats are replaced. To recover this, the discarded habitats are stored in a vector (rock vector) that tracks the habitats.
In BiMO, when the number of islands is very small, during the processes of immigration and emigration, the generation of duplicates is a frequent event. Instead of generating new habitats randomly, they are taken from the best habitats belonging to the rock vector. The proposed algorithm is described by the following pseudo-code:
(0) initialize randomly each SIV (design variable) of the habitats and rank habitats according to the generalized fitness;
(1) elitism: save the best two habitats;
(2) given the immigration/emigration curves, decide for each SIV whether emigrate or immigrate;
(3) if emigration occurs: store the lost habitat in the rock vector;
(4) select habitats for mutation: the lost habitats (before mutation) are stored in the rock vector;
(5) evaluate the objective functions;
(6) if duplicates occur, substitute them with the best habitats taken from the rock vector;
(7) rank habitats according to the generalized fitness;
(8) if the stopping criterion is not reached, go to (1), else stop.
The 𝜇-BiMO was applied with 7 habitats and 100 iterations, i.e. 700 objective function calls were done.
In Figs 5–7 the results obtained with both methods and for homogeneous and non-homogeneous cases are shown. In particular, the improved Pareto fronts found for the three different cases in Fig. 2 and obtained applying the M-NSGA and 𝜇-BiMO algorithms are compared. The stars represent the individuals on Pareto front obtained applying M-NSGA and the circles the ones obtained applying 𝜇-BiMO algorithm. In each figure small panel on the right are a zoom of a part of the Pareto front depicted on the panel at left. The starting points are taken randomly and the stopping criterion is the maximum number of objective function calls, set to 700.
Optimization results, homogeneous case like in Fig. 2(a). Pareto front obtained with M-NSGA and 𝜇-BiMO algorithm (stars represent the individuals on pareto front obtained applying M-NSGA and the circles the ones obtained applying 𝜇-BiMO algorithm). Small panels are a zoom of the Pareto front evidenced in the darker rectangle. EP1 and EP2 are points representing the two solutions reported in Table 1.
Optimization results, non-homogeneous case like in Fig. 2(b) with the gray rectangle made with the electrical characteristic of the Gel A. Pareto front obtained with M-NSGA and 𝜇-BiMO algorithm (stars represent the individuals on pareto front obtained applying M-NSGA and the circles the ones obtained applying 𝜇-BiMO algorithm). Small panels are a zoom of the Pareto front evidenced in the darker rectangle. EP1 and EP2 are points representing the two solutions reported in Table 1.
Optimization results, non-homogeneous case like in Fig. 2(c) with the gray rectangle made with the electrical characteristic of the Gel B. Pareto front obtained with M-NSGA and 𝜇-BiMO algorithm (stars represent the individuals on pareto front obtained applying M-NSGA and the circles the ones obtained applying 𝜇-BiMO algorithm). Small panels are a zoom of the Pareto front evidenced in the darker rectangle. EP1 and EP2 are points representing the two solutions reported in Table 1.
The values of design variables and objective functions of two end-points (EP) belong to Pareto front of each case (homogeneous, Gel A and Gel B) are shown in Table 1. The end points are defined as end points of the Pareto front obtained by merging the 𝜇-BiMO and M-NSGA results.
The 𝜇-BiMO algorithm is able to well approximate the end-points of the Pareto front in all the three cases; however, some solutions found by the algorithm are dominated by the solutions found by M-NSGA.
On the other hand, the M-NSGA algorithm is able to approximate the whole Pareto front but in some cases the end-points are not accurately found, in comparison with 𝜇-BiMO.
Figure 8 shows the maps of electric field intensity where the equifield lines are represented for the three points labelled EP2 in Figs 5–7. The panel (a) represents the homogenous case where the insertion is considered made on potato tissue, whereas in panel (b) and (c) the inhomogeneous cases are represented, where the insertion is made of gel A and gel B, respectively. In each map, the square represents the sampling area for electric field and the arrows show the line at 300 V/cm. In all the cases, the electric field distributes in a different way depending on the conductivity of the insertion. In the three cases, the optimized needle distances are different as well as the applied voltages. Considering the ratio between the voltage and the needle distance, reported in Table 2, it appears that if the minimum of f 1 is considered (EP1 points) the required voltage with respect the optimal distance (ratio U/D) is lower with respect to the ones considered minimizing the function f 2 (EP2 points). In particular, focusing on EP2 points, the optimized geometry in the case of Gel A shows a U/D ratio that is about twice with respect to the homogeneous case. Moreover, considering a less conductive insertion (gel B with conductivity 0.24 S/m) the applied voltage is a little bit higher with respect to the homogeneous case.
Values of the ratio U/D in [V/cm] for the cases in Table 1
From the optimization results and focusing on EP2 points, it appears that inhomogeneous cases require different applied voltages in order to obtain an electric field intensity that overcome the prescribed threshold of electroporation with respect to homogeneous case. Moreover, the values of the applied voltage depends on the conductivity of the material of the insertion with respect to the one of the bulk.
In the paper a multi-objective optimization for the design of the electroporation set-up in terms of applied voltage and needles position and diameter in order to improve the distribution of electric field is proposed. The target was a distribution of the electric field that overcame the electric field threshold in a tissue with inhomogeneity in electrical properties in the treatment area limiting the area interested by irreversible electroporation. The optimization results show that both different electrode geometry and voltage amplitude are suitable for different arrangement of tissue inhomogeneity in order to improve the electric field intensity in the analysed area. In particular, if the inhomogeneity has an higher conductivity the ratio between the applied voltage and the needle distance is larger with respect to the for a lower conductivity inhomogeneity.
The results obtained with M-NSGA and µ-BiMO algorithms are shown. Both the two algorithms found solution belong to the same Pareto front. In particular, 𝜇-BiMO found solutions on the extremes of the front, whereas M-NSGA is able to found intermediate solutions. The proposed algorithms show results considering both homogeneous and the non-homogenous cases.