Augmented Ant Colony Algorithm for Virtual Drug Discovery

Docking is a fundamental problem in computational biology and drug discovery that seeks to predict a ligand’s binding mode and aﬃnity to a target protein. However, the large search space size and the complexity of the underlying physical interactions make docking a challenging task. Here, we review a docking method, based on the ant colony optimization algorithm, that ranks a set of candidate ligands by solving a minimization problem for each ligand individually. In addition, we pro-pose an augmented version that takes into account all energy functions collectively, allowing only one minimization problem to be solved. The results show that our modiﬁcation outperforms in accuracy and eﬃciency.


Introduction
A major challenge in computational chemistry and drug discovery is the docking problem, which consists of identifying the optimal binding orientation and position between ligand molecules and a receptor protein, as well as ranking molecules according to their binding affinity [1][2][3].This problem is of great importance, as it can aid in revealing the docking mechanism at the atomistic level and in developing novel drug candidates [4][5][6][7].
From a mathematical point of view, the docking problem is a minimization problem, i.e. a type of optimization problem in which the objective is to find the minimum value of a given function.In the context of the docking problem, the target is to minimize an energy function that describes the strength of the binding between the ligand and the receptor, where the global minimum represents a particularly stable configuration that allows the activation of biochemical mechanisms of the receptor [8].Unfortunately, finding the global minimum is hampered by the presence of local minima and energy barriers that prevent the exploration of the energy function (trapping problem).A second obstacle is the high dimensionality of the energy functions, and thus the vast size of the space in which the solution resides.Indeed, energy functions in docking problems have at least six degrees of freedom: three rotational and three translational degrees of freedom that describe the ligand's orientation and position with respect to the receptor.Further, more sophisticated models incorporate the torsions of molecules' bonds that are not associated with a ring-chain, adding even more degrees of freedom.
In a nutshell, finding the global minimum is a hard computational task that requires advanced algorithms and techniques to optimize computational resources without losing accuracy.Ongoing research in this field has proposed a vast spectrum of methods and strategies to improve the accuracy and efficiency of docking algorithms.For instance, the programs GOLD [9] and AutoDock [10] use genetic algorithms, ICM [11] and QXP [12] use Monte Carlo minimization, PRO LEADS [13] implements simulated annealing, evolutionary programming, and tabu search, and FlexX [14] and DOCK [15] use fragments algorithms.
In this article, we review the stochastic optimization method PLANTS (Protein-Ligand ANT System) [16][17][18] implemented in the software platform VirtualFlow [19,20] for virtual screening and ligand preparation.Additionally, we propose an augmented version of PLANTS, which we will refer to as PLANTS+.PLANTS is an Ant Colony Optimization algorithm [21], i.e. an algorithm inspired by the behavior of ants that search for short paths between the anthill and food sources, guided by pheromone trails.This kind of algorithm was first designed for discrete optimization problems that can be represented by a graph made of nodes and edges.A typical example is the traveling salesman problem, where the shortest path connecting two network nodes is sought.The application to continuous function minimization problems is possible by discretizing the domain of the function in disjoint subsets that are analogous to the network nodes in the traveling salesman problem.The artificial ants, however, do not search for the shortest path between two subsets, but rather search the subset of the domain that contains the global minimum, and then deposit the pheromone in areas with low energy values.
Given a set of ligand candidates for a specific receptor, PLANTS is used to determine the global minimum of each binding energy function.The ligands are then ranked according to their energy values when they are in their optimal position and orientation, and the best ligand for the receptor is the one that takes the lowest energy value.PLANTS+ extends the search for the global minimum from the space of degrees of freedom to the space of binding energy functions by analyzing all energy functions simultaneously, rather than individually.Indeed, the artificial ants deposit the pheromone in the regions of the space of degrees of freedom where the energy value is lowest among all energy functions, and the algorithm returns not only the coordinates of the global minimum but also the energy function that takes the lowest value from those coordinates.Then, the main difference between PLANTS and PLANTS+ is that the former solves as many minimization problems as the number of ligands, whereas the latter solves one unique minimization problem, as illustrated by the workflows in fig. 1.
In this work, we review PLANTS and describe the new version PLANTS+ in section 2. Afterward, we present an illustrative example that highlights the differences between the two methods and present the results of a comparative analysis 3 that show under which conditions PLANTS+ is more accurate and more efficient than PLANTS.Section 4 concludes the article with a short summary of our findings as well as an outlook for future research.

Methods
Consider a ligand-receptor system, where the receptor is kept fixed in the Euclidean space, while the ligand can rotate or translate with respect to the receptor determining 3 rotational degrees of freedom and 3 translational degrees of freedom.The system is characterized by a 6-dimensional binding energy surface where Ω denotes the space of all possible positions and orientations of the ligand with respect to the receptor, while the coordinates Ψ i and r i denote respectively the rotational and translational degrees of freedom.Note that we restricted the domain of the translational degrees of freedom to positive values belonging to R ≥0 .Consider now a set of N ligands ligands, each with an energy function where the index α = 1, 2, ..., N ligands denotes the αth ligand of the set, and the notation used for the coordinates Ψ α,i and r α,i remarks that each ligand has its own degrees of freedom.We are interested in finding the energy function that takes the lowest value in its global minimum.

PLANTS
In Plants, the global minimum of each binding energy function is determined using the Ant Colony Optimization algorithm.In order to solve problems involving continuous function minimization, the domain must be discretized.
Then, for a docking problem, the rotational degrees of freedom Ψ i and the translational degrees of freedom r i need to be discretized ∀i = 1, 2, 3. Various discretizations are possible and affect the efficiency and accuracy of the algorithm.However, when no a priori knowledge of the system is available, it is recommended that the degrees of freedom be discretized into equal intervals.For example, a rotational degree of freedom Ψ with Ψ(0) = −π and Ψ(N Ψ ) = π can be discretized into N Ψ equal disjoint subsets such that while a translational degree of freedom r with r(0) = 0 and r(N r ) = r max can be discretized into N r disjoint subsets such that where r max ∈ R ≥0 is the maximum distance between ligand and receptor appropriately chosen.This gives rise to a discretization of the space Ω into where the notation of the indices l i and k i indicates that for each degree of freedom, we can choose different indices l and k in order to create every possible combination.The discretization subsets of the space Ω of the degrees of freedom are analogs to the network nodes in the traveling salesman problem, and the artificial ants move from one subset to another one looking for the subset A * ∈ Ω that contains the coordinates Ψ * i , r * i , ∀i = 1, 2, 3, such that the energy function is minimized: Consequently, in this context, an artificial ant is a possible solution to the optimization problem, i.e. a specific combination of degrees of freedom that minimizes the energy function.The ants are randomly drawn in accordance with six probability vectors p Ψi and p ri , ∀i = 1, 2, 3, respectively generated for each rotational and translational degree of freedom.The probability vectors have a number of entries equal to the number of subsets in which the degrees of freedom have been discretized, and are obtained by normalizing six corresponding pheromone vectors τ Ψi and τ ri , ∀i = 1, 2, 3.

Pheromone update rule
Initially, all pheromone vector entries are equal to one, so that the probability vectors approximate uniform distributions and the N ants solutions Ψ j i , r j i , ∀i = 1, 2, 3 are generated with equal probability.Here the index j denotes one of N ants possible solutions.These solutions are then locally minimized: Local minimization can be achieved by means of several algorithms.PLANTS uses the Nelder-Mead algorithm [22], a simplex algorithm that does not need the calculation of derivatives, but other choices are available.Afterwards, the best artificial ant, i.e. the jth solution with coordinates Ψ * i and r * i with the lowest energy value, is identified and used to update the vectors τ Ψ and τ r .This ensures that the global minimum is approached step by step by generating new solutions in accordance with the probability vectors.
The rule for updating the pheromone vectors τ Ψ and τ r is the same as described in [20].Given the best solutions Ψ * and r * , the lth and kth entries of the pheromone vectors respectively associated with the rotational degree of freedom Ψ i and the translational degree of freedom r i , are updated as where the functions 1 Ψi (l i ) and 1 ri (k i ) are given by In eq. 8, ρ is an input parameter between 0 and 1 called evaporation rate.The function of the evaporation rate is to partially reset the pheromone vector to prevent the algorithm from being trapped in a local minimum.Finally, the corresponding probability vectors are obtained by normalization: The procedure runs for N steps iterative steps by generating new solutions which converge to the global minimum.The algorithm is finally applied to all of the ligands under consideration, which are ranked based on their global minimum energy value.Depending on the infrastructure, this can be implemented using a for loop or by parallelizing the process.The complete algorithm is described in alg. 1.

PLANTS+
The PLANTS algorithm searches for the optimal position and orientation of each candidate ligand by analyzing each energy function individually.The best candidate is finally determined by evaluating the energy value of each ligand in its global minimum.This corresponds to solving N ligands optimization problems, i.e. the number of ligands considered, where each optimization problem consists of minimizing an energy function E(Ψ i , r i ) with 6 degrees of freedom.Alternatively, we consider a single optimization problem in which the function to be minimized is constructed as the sum of all energy functions defined in eq.2: This function, which we will refer to as the global energy function, has 6 × N ligands degrees of freedom, so the new problem consists of minimizing a 6 × N ligands -dimensional function.
end for 21: end for 22: return ligand with the lowest energy in its global minimum.
To solve this problem, we propose an augmented version of the original PLANTS algorithm, which we will refer to as PLANTS+.In practice, PLANTS is modified as follows.The outermost for-loop, over the candidate ligands, is removed, and the simplex algorithm is applied to the global energy function defined in eq.11.Then, the pheromone and probability vectors are updated according to the rule defined in eq. 8.The crucial difference from the original algorithm is that the solutions (artificial ants) are minimized with respect to the global energy function.This entails that the pheromone and probability vectors contain information on all ligands taken together, and not individually as in the original PLANTS.The algorithm is fully illustrated in alg.2, where we denote in line 15 the best solution as Ψ * L and r * L in alg. 2 since it is calculated with respect to all ligands.
The limitation of this approach is the high number of dimensions of the energy function defined in eq.11, which grows with the number of ligands examined.It is in fact well known that the simplex algorithm loses accuracy with the dimensionality of the function to be minimized [23].On the other hand, because the ligands' energy functions (eq. 1) are independent of each other, there is no need to apply the simplex algorithm to the global energy function (eq.11).Instead, it is more convenient to apply the local minimization algorithm to the individual energy functions either within a for loop, as in alg.2, or by paralleling the calculation.// The following double for loop can be parallelized 9: for i = 1 to N ligands do 10: for a = 1 to N ants do p Ψ ← normalize(τ Ψ ); 18: p r ← normalize(τ r ) 19: end for 20: return ligand with the lowest energy in its global minimum.

Illustrative example
To present the main differences between PLANTS and PLANTS+, we consider the problem of minimizing five two-dimensional functions representing the binding energy functions of five ligands interacting with a receptor.The energy functions are written as where µ 1 , µ 2 , µ 3 , µ 4 , µ 5 , µ 6 ; ϵ, r 0 , ν 1 , ν 2 , ν 3 , ν 4 are parameters randomly generated from a uniform distribution (see supplementary information for details).
The function E(Ψ, r) can be decoupled into two components E Ψ (Ψ) and E r (r) which depend on two degrees of freedom corresponding to a rotational and a translational degree of freedom.The pairs of five functions E Ψ (Ψ) and E r (r) are illustrated respectively on the left and the right column of fig. 2 and are characterized by irregular sequences of valleys and ridges to make the search for the global minimum more difficult.In particular, the functions E Ψ (Ψ) of the rotational degree of freedom are periodic and are realized by superposing sinusoidal functions with low and high frequencies.Similarly, the translational degree of freedom functions E r (r) are realized by means of the Lennard-Jones potential, which well describes the interaction between non-bonding molecules, perturbed by sinusoidal functions to generate local minima.The functions greatly differ from each other due to the randomness of the parameters and the global minima, marked by the green dots, are located in different positions.We also mark that the height of the barriers of the functions E Ψ (Ψ) are considerably higher than the barriers of the function E r (r): approximately 3 kJ mol −1 versus 10 kJ mol −1 .PLANTS and PLANTS+ were applied with the same input parameters: the evaporation rate was ρ = 0.1, the number of solutions, i.e. the artificial ants generated at each iteration was N ants = 5, and the total number of iteration was N steps = 80.The range [−π, π] rad of the rotational degree of freedom and the range [0, 10] nm of the translational degree of freedom were both discretized into 200 equal subsets.In order to validate the results, we used the brute force algorithm scipy.optimize.brute,included in the SciPy library [24], to find the exact global minima of all the functions.

PLANTS
We applied PLANTS, in its original form, seeking the global minimum of all five functions.For each energy function, we generated N ants = 5 artificial ants, i.e. five solutions Ψ and r uniformly distributed respectively in the range [−π, π] rad and [0, 10] nm.The initial solutions are marked with red dots in fig.3-a, where we also plot the solutions generated at each generation (black lines) and the coordinates of the best solutions Ψ * and r * (blue lines), i.e. the solutions which take the lowest energy value after being minimized.The green dots denote the global minima.Correspondingly, fig.3-b shows the time evolution of the energy E(Ψ * , r * ) of the best solutions.
Analyzing the results, we observe that the algorithm converges slower with rotational than translational degree-of-freedom functions.Indeed, the coordinates Ψ * stabilize after approximately 30 iterations, while the coordinates r * need approximately 10-20 iterations.The different behavior is due to how the probability vectors p Ψ (Ψ) and p r (r), associated with the pheromone vectors, update during the search for minima.In the beginning, the probability vectors describe a uniform distribution, and ideally, they converge over time to a delta function, with the peak located in a neighborhood of the global minimum.Thus, solutions generated subsequently have a higher probability of belonging to a subset that contains the global minimum.However, if the function under consideration has several local minima that assume almost the same energy value as the global minimum, then the pheromone will be distributed over a wider spectrum and new solutions will be generated accordingly.This is what happens with the energy functions E Ψ (Ψ) under consideration.As an example, we report in fig. 4 the time evolution of the probability vectors p Ψ (Ψ; t) (left side) and p r (r; t) (right side) related to the second energy function, where the index t denotes a specific iteration step.We observe that for the rotational degree of freedom, the convergence of the probability vector takes more time despite the average height of the barriers of E Ψ (Ψ) being greater than E r (r).What makes the determination of the global minimum of E Ψ (Ψ) difficult is the presence of four local minima at Ψ ≈ −2.2, −0.5, 0.7, and 2.5 rad that take approximately the same energy value.Only after 40 iterations does the probability associated with the exact global minimum at Ψ ≈ −0.5 rad dominate the rest.In contrast, the global minimum of the function E r (r) is less ambiguous, and the probability p r (r) converges with less than 10 iterations, identifying the global minimum in the range r ∈ [2,4] nm.Taking a closer look, we identify three close peaks in this interval due to the serrated pattern that characterizes the second function E r (r), but, after a further 30 iterations, only the peak associated with the exact global minimum survives.
In this example, after about 40 iterations, the global minimum for each function is identified and the functions can be ranked according to their energy value shown in fig.3-b.Thus, we conclude that the second ligand, represented by a red line in the figure, is the most suitable for the receptor in consideration.

PLANTS+
We applied PLANTS+ to the same five functions E(Ψ, r) examined in the previous section.However, unlike PLANTS, we did not use PLANTS+ for finding the global minimum of each energy function individually, but for the global minimum of the global energy function defined in eq.11.For this, we defined only one probability vector p Ψ (Ψ) for the rotational degree of freedom, and one probability vector p r (r) for the translational degree of freedom, which were updated at each iteration accordingly to the best solution (Ψ * L , r * L ) with respect to the orientation, the position, and the energy function.The probability vectors were used to generate 5 solutions (Ψ, r) over the range [−π, π] rad and [0, 10] nm, that converge to the global minimum.
In fig.5-a, we show the time evolution of the coordinate Ψ * L of the best solution on the left side, and of the coordinate r * L on the right side.During each iteration, the lines change color to emphasize which energy function assumes the lowest energy value.In this example, initially, the lines are colored blue since the first ligand takes the least value, then the color becomes red between steps 2 and 7, then it turns blue again.After the 25th iteration step, the color remains red, as the second energy function takes the lowest value.
Fig. 5-b shows the time evolution of the energy.Precisely, this graph shows the energy value taken by the best solution (Ψ * L , r * L ), at a given iteration step, using the energy function that returns the lowest value.In accordance with the graph depicting the evolution of the solution, the best energy function is used to color the line.At first, the algorithm identifies the first ligand as optimal for the receptor, but at the end of the execution, the algorithm returns the second ligand, in accordance with results obtained by PLANTS.
We also observe that between steps 0 and 25, the competition between energy functions involves ligands 2, 4, and 5. Similarly, we observed the same competition with PLANTS in fig. 3.This is due to the fact that the global minima of each of these ligands assume approximately the same value, as shown in fig. 2.

Comparative Analysis
With the next numerical experiment, we investigate how PLANTS and PLANTS+ depend on the input parameters, and we seek under which settings PLANTS+ is more accurate and more efficient than PLANTS.For this purpose, we randomly generated a set of 100 energy functions E(Ψ, r) using eq.12 and we searched for the energy function which takes the least energy value at its global minimum.The choice of the energy functions was made in order to make difficult the optimization problem, not only generating random sequences of local minima and barriers but distributing the global minima over a large range, as illustrated in fig.6.We studied 64 different sets of parameters defined by the combination of the following input parameters: A total of 20 runs of each algorithm were performed on each set of parameters in order to produce meaningful statistical results.For each set of parameters we measured and compared the execution time and the success rate.The execution time is the real-time in seconds, averaged over all 20 repetitions, required to complete the instructions described in algs. 1 and 2. To evaluate the success rate, we first applied the brute force algorithm scipy.optimize.bruteto each energy function to find the exact global minima.Then we defined a prediction, realized by PLANTS or PLANTS+, as "successful" if the relative error was smaller than 5%, and the success rate of a specific set of parameters as the percentage of successful predictions over all 20 repetitions.

PLANTS
Fig. 7-a reports the results for the algorithm PLANTS.The leftmost panel shows 64 points representing the execution time (horizontal axis) versus the success rate (vertical axis) of each set of parameters.The points have been colored according to the number of generated solutions N ants , thus we identify four clusters: using N ants = 1, PLANTS never converges, with N ants = 5 the success rate varies between 40% and 80%, with N ants = 10 the success rate is between 80% and 90%, and with N ants = 20 is between 85% and 95%.Correspondingly the execution time increases from approximately 10 seconds to 150 seconds.This result is explained by the fact that more solutions require more calculations, especially during the execution of the simplex algorithm.However, a large number of solutions also facilitates the exploration of the domain and the search for the global minimum.
Likewise, we could color the points to highlight the dependence on the parameters ρ and N steps , but for these two parameters, the figures would be unreadable.Thus, we decided to consider only the set of 16 algorithm executions realized with N ants = 5 solutions, and we plot in the middle and rightmost panel respectively the dependence on the evaporation rate ρ and on the iteration steps N steps .We observe that ρ = 0.4 (green) and ρ = 0.3 (red) have no impact on execution time, while the success rate is approximately 60% and 70%.This suggests that small evaporation rate values favor accuracy.Intuitively, the reason for this is that if the evaporation rate is too high, the pheromone and probability vectors are reset too quickly, and the information acquired during the first few iterations is lost.By further decreasing the evaporation rate to ρ = 0.2, the success rate is improved to 80% only if N steps = 30, 40, 50.With ρ = 0.1, the dependence on N steps is further emphasized, and the success rate reaches 85% only if N steps is equal to 50.Thus, low evaporation rate values require longer algorithm runs to get out of any local minima in which the system could be trapped.Analyzing the rightmost panel, we deduce that the number of iteration steps mainly impacts the execution time and the difference in success rate is due to the parameter ρ.
In essence, we conclude the following: ❼ Increasing the number of solutions improves accuracy, but the required execution time exponentially increases.❼ Small evaporation rates improve the accuracy, but only if the number of iteration steps is large enough.❼ Conversely, the number of iteration steps affects the execution time, but it improves the accuracy only in combination with a small evaporation rate.

PLANTS+
We out the same analysis with PLANTS+, whose results are reported in fig.7-b.The color of the points in the graphs highlights how success rate and execution time depend on the parameter N ants (left panel), ρ (middle panel), and N steps (right panel).As for PLANTS, we fixed the number of solutions N ants = 5 to better illustrate the dependence on ρ and N steps .
With regard to the number of generated solutions N ants , we observe the same behavior as PLANTS: accuracy increases with N ants , but at the expense of an exponentially increased execution time.Despite the difficulty of reading the second panel, we can still discern a cluster associated with ρ = 0.3 (red) above ρ = 0.4 (green).Thus, reducing the evaporation rate increases accuracy, but execution time depends on the other parameters.Lastly, the number of iterations N steps affects execution time primarily.Then, we conclude that PLANTS and PLANTS+ are dependent on the input parameters in the same manner.
In fig.7-c, we compare the success rate and execution time of the two algorithms, in order to determine the algorithm that is more efficient and accurate.The graph reports the difference (success rate of PLANTS+) − (success rate of PLANTS) , on the vertical axis, and (execution time of PLANTS+) − (execution time of PLANTS) , on the horizontal axis.Hence, a positive value on the vertical axis and a negative value on the horizontal axis indicate that PLANTS+ is respectively more accurate and more efficient than PLANTS using the same parameters N ants , ρ and N steps .In all cases, PLANTS+ executions were more efficient than PLANTS: the gain varies from 5 seconds with N ants = 5 (black dots) to 20 seconds N ants = 15 (green dots).Accuracy worsened in 16 cases, those below the dashed line, but improved in all the other 48 cases.In particular, we observe that in 5 cases with N ants = 1, in 12 cases with N ants = 5, in 15 cases with N ants = 10, and in all cases with N ants = 15, the success rate improves.Then, PLANTS+ benefits from increasing the number of solutions N ants both in accuracy and efficiency.The rightmost panel reveals another important result.With the same settings, PLANTS+ is more efficient when the number of iteration steps is increased, and the gain in efficiency is proportional to N steps .In other words, with PLANTS+ the N parameter has less impact on execution time than with PLANTS.

Conclusion and outlook
In this article, we reviewed the optimization method PLANTS and proposed an augmented version, referred to as PLANTS+ in the text, for finding the optimal orientation and position of a ligand at the receptor's binding site.The main difference between the two methods is that PLANTS examines the binding energy functions of ligands individually and determines the best ligand for a given receptor by comparing the energy that each ligand takes at its optimal pose, while PLANTS+ solves a single optimization problem and searches for the optimal position and orientation by comparing the energy functions of the ligands at each iteration step.This modification permits better use of computational resources because the probability vector used to generate new possible solutions incorporates more information than PLANTS.In fact, the probability vector is not updated by taking into account only the degrees of freedom, but by evaluating at each iteration step which energy function takes the lowest value.
We studied how the input parameters affect the accuracy and the efficiency of both algorithms and our results show that the number of solutions generated, i.e. the input parameter N ants , is decisive for the success of the global minimum search, but that this parameter also significantly affects the execution time.However, with the same settings, PLANTS+ is always more efficient than PLANTS.Thus, the loss in efficiency that would occur with PLANTS using a high value for the N ants parameter is recovered by PLANTS+ which returns the same result in less time.
This work demonstrates that PLANTS+ is more accurate and efficient than PLANTS, and future work will focus on the implementation with Virtualflow [19], in order to contribute to the discovery of novel drug candidates and to the optimization of existing drugs.We emphasize that in this work we have considered the ligand-protein system as static.However, further development of the algorithm will take into account the inclusion of kinetic observables, e.g.transition rates, as a criterion for classifying ligands.This will provide a more accurate screening of ligands [25] and will make it possible to study the dependence of docking on external variables [26], e.g. the acidity of the cell membrane, a variable that is considered crucial in the development of undesirable effects [27,28].

Declarations
Ethics approval and consent to participate.This article does not contain any studies with human participants or animals performed by any of the authors.
Consent for publication.All authors gave their consent.

Fig. 1
Fig. 1 Overview of the workflow of the algorithm PLANTS and PLANTS+.

Fig. 3
Fig. 3 PLANTS.(a) Time evolution of the coordinates (black lines) for each potential: Ψcoordinate (left), r-coordinate (right); one ligand per row.The blue line marks the evolution of the best solution at each iteration.The red dots mark the initial solutions, the green dots mark the global minima.(b) Evolution of the energy of the best solution for each ligand.

Fig. 4
Fig. 4 PLANTS.Time evolution of the probability vectors p Ψ (Ψ; t) (left) and pr(r; t) (right) of the second energy function of the illustrative example.Dark and light blue are associated respectively with a low and high probability.

Fig. 5
Fig. 5 PLANTS+.(a) Evolution of the Ψ-coordinate (left) and r-coordinate (right) of the best solution with respect to all energy functions.The color denotes the energy function which takes the lowest value during a specific iteration step.(b) Time evolution of the energy of the best solution.

Fig. 7
Fig. 7 (a) Success rate (vertical axis) and execution time (horizontal axis) of PLANTS; (b) Success rate (vertical axis) and execution time (horizontal axis) of PLANTS+; (c) Difference in execution time and success rate between PLANTS+ and PLANTS.The points falling within the rectangle with dashed edges indicate that PLANTS+ is more accurate and more efficient than PLANTS using the same input parameters.The legend box of figure (a) applies also to figures (b) and (c).