WO2003105080A2

WO2003105080A2 - Metapopulation genetic algorithm for combinatorial optimisation problems

Info

Publication number: WO2003105080A2
Application number: PCT/BE2003/000105
Authority: WO
Inventors: Michel C. Milinkovitch; Allan R. Lemmon; Cédric SZPIRER; Philippe Gabant
Original assignee: Universite Libre De Bruxelles
Priority date: 2002-06-10
Filing date: 2003-06-10
Publication date: 2003-12-18
Also published as: EP1520254A2; WO2003105080A3; AU2003240321A1

Abstract

The present invention relates to a method for obtaining an efficient solution to a combinatorial optimisation problem, comprising the following steps: a) selecting different groups (A, B, C) of population able to interact with each other; b) applying specific selection parameters to each group of population; c) obtaining for each group of population, a consensus selection of optimal solution portions and portions of solution to be optimised, said solutions satisfying the said parameters; d) comparing by a consensus pruning the solutions obtained from the different groups of population, the optimal solution portions present in the majority of the different groups being not modified; e) repeating steps c) and d) for the other solution portions to be optimised until said portions are optimised; f) selecting from said preceding steps one or more complete optimal solutions for said combinatorial optimisation problem. The present invention relates also to the algorithm for performing the method according to the invention.

Description

METAPOPULATION GENETIC ALGORITHM FOR COMBINATORIAL OPTIMISATION PROBLEMS

Field of the invention

[0001] The present invention is related to a new method for selecting efficient solutions of combinatorial optimisation problems in different industrial processes such as transport and telecommunication networks, designing optimal production processes, designing optimal management, logistic and production decisions, DNA-array technology, neuronal systems, etc.

Background of the invention [0002] Many of the optimisation problems are "NP-complete" , which means that no future computer will ever be able to solve them exactly in practical computing time because the number of possible solutions increases explosively with the number of elements of the problem. Hence, many heuristic approaches (i.e., methods that are fast but find approximate solutions) have been developed, with varying degrees of success .

[0003] Optimality-criterion based philogeny inference is a perfect example of "unsolvable" problem. Indeed, the number of solutions increases explosively with the number of taxa (taxa are groups of individuals (i.e., populations), and can be species, genera, families, etc., depending on the question addressed by the philogeny investigation) . The total number of possible unrooted, bifurcating tree topologies (i.e., solutions) among T terminal taxa is:

[0004] This corresponds to nearly 32 billion different trees for 14 taxa and 3X10⁸⁴ trees (i.e., more than the number of atoms in the known universe) for 55 taxa. Philogeny inference is a combinatorial optimisation problem of "NP-complete" type { e . g. , (1)) because (i) no known algorithm can solve it in polynomial time (i.e., with a running time growing no faster than some fixed power of the number of elements, here, the number of taxa) , and (ii ) demonstrating the existence of such an algorithm would imply that all NP problems have a polynomial time solution. As most mathematicians expect that no such algorithm exists, one is forced to admit that no future civilisation will ever build a computer capable of solving the problem while guaranteeing that the optimal solution has been found.

[0005] Computation time has rapidly become a central and widespread practical problem in philogeny inference for two main reasons. First, major advances in molecular biology and biotechnology have caused a dramatic increase in the number of DNA sequences available, stimulating researchers to increase their taxon sampling when performing tree reconstruction. Incidentally, several authors ( e . g. , (2)) have suggested that the accuracy of philogeny inference increases with increased taxon sampling. Second, many simulation studies in the last 10 years ( e . g. , (3)) have identified the maximum likelihood

(ML) criterion (4) as one of the best for philogeny inference. Reasons for this include (i ) statistical consistency; ML estimators tend to converge to true parameters values when the number of characters is increased, (ii) robustness; violations of the ML model assumptions have only a moderate impact on the tree inference accuracy, (Hi) the ability to compare different trees within a statistical framework, and (iv) the ability of ML to make full use of the original character matrix. Unfortunately, the analytical power of ML has a cost : computation time. Hence, application of this model-based criterion makes the use of exact philogeny inference methods impractical for more than a trivial number of taxa

(about 12) while the number of living species is estimated between 10 and 100 millions) .

[0006] As with many NP-complete problems, the only practical approach to large philogeny inference is the use of heuristics, i.e., methods that find one or several solution (s) faster than exact searches while sacrificing a guarantee that the optimal solution will be found. In response to this trend, much of the current research in computational philogenetics concentrates on the development of more efficient heuristic approaches.

[0007] Many existing heuristic methods are handicapped with the need to optimise model parameters

(such as branch lengths) on each solution (here, each tree) examined, a procedure that significantly slows computation time. One such heuristic, which is implemented in the most widely used inference packages (PAUP* (5) , PHYLIP (6) , etc.), is based on "hill climbing": an initial solution

(most commonly obtained using the "stepwise addition" algorithm (7)) is subjected to a perturbation (here, a topological rearrangement known as "branch-swapping" (7)), followed by an optimisation of model parameters. If the new solution has a better score than the starting solution, it is kept and used as the new starting solution - otherwise the proposed solution is rejected. The procedure is repeated in a systematic fashion, until no swapping on the current best solution can improve the score . Due to the need to optimise parameters at every step (intra-step optimisation) , simple perturbation methods remain computationally impractical for more than 25-50 taxa, depending on the complexity of the ML model . [0008] Two classes of solutions have been developed in response to the problems associated with optimising parameters on large trees. The first of these classes includes solutions that partition the large problem into many small sub-problems whose solutions are then combined into a consensus global solution. For example, the "quartet puzzling" method (8) works by reconstructing the best ML tree of each possible quartet of taxa, then combining

(i.e., puzzling) the quartet trees to construct an overall

T-taxon tree. The advantage of quartet puzzling is that it avoids numerous optimisations on a large number of elements, optimising instead numerous small sets of elements, which is computationally much simpler. Despite the initial appeal, the method has proven to be only moderately more accurate (8) than simple clustering method that requires significantly less computation time.

[0009] The second class is comprised of stochastic heuristics that avoid optimisation of numerous solutions entirely. Instead, they incorporate methods that allow model parameters to be optimised as the search proceeds - taking an inter-step optimisation strategy. Stochastic Simulated Annealing (SSA (9) ) , for example, is based on the simple perturbation algorithm described above, but incorporates a method that perturbs model parameters at each iteration instead of requiring optimisation of each potential solution. SSA avoids local optima by accepting changes that decrease the likelihood of the solution with a probability inversely proportional to the reduction in likelihood. Another increasingly popular approach is the

Bayesian method (10) based on the Markov-Chain Monte-Carlo

(MCMC) algorithm. MCMC-based methods also benefit from avoiding intra-step optimisation, although they have a slightly different aim: sampling the distribution of the space of solutions instead of only finding optimal solutions .

[0010] Finally Matsuda (11) , Lewis (12) , and Katoh et al . (13) have recently applied the genetic algorithm

(GA) , a type of evolutionary computation method (14) , to philogeny inference. GAs implement a set of operators that mimic processes of biological evolution such as mutation, recombination, selection, and reproduction. After an initial step of generating a population, the individuals

(specific solutions) within that population are (i ) subjected to mutation and/or recombination, and (ii ) allowed to reproduce with a probability that is a function of their relative fitness value. In the case of a philogenetic inference problem, individuals are typically composed of topologies and model parameters ( e . g. , branch lengths, the transition/transversion ratio, rate heterogeneity parameters, etc.) that need to be optimised. A mutation is a stochastic alteration of one component of the individual ( e . g. , a topological rearrangement, a change in one branch length, or a random modification of a model parameter) , and the fitness of an individual is a function of the score for that tree. Since selection retains the changes that improve the value of the optimality function, the mean score of the population of trees improves over time, i.e., across "generations". Aims of the invention

[0011] The present invention aims at providing a simple and efficient method as well as an algorithm used in said method able to find and obtain excellent solutions to combinatorial optimisation problems in a suitable fraction of time (approximately 100 to 1000 times faster of that required by existing methods of the state of the art.

Summary of the invention [0012] The present invention relates to a method for obtaining an efficient solution to a combinatorial optimisation problem, comprising the following steps: a) selecting different groups (A, B, C) of population able to interact with each other; b) applying specific selection parameters to each group of population (necessary for performing a consensus selection of different solutions) ; c) obtaining for each group of population, a consensus selection of optimal solution portions and portions of solution to be optimised, said (tested) solutions satisfying the requirements of said parameters; d) comparing by a consensus pruning (communication tools among the populations) the solutions obtained from the different groups of population, the optimal solution portions present in the majority of the different groups being not modified; e) repeating steps c) and d) for the other solution portions to be optimised until said portions are optimised; f) selecting from said preceding steps one or more complete optimal solutions (when the portions of solution can not be further optimised) for said combinatorial optimisation problem. [0013] Advantageously, the method can be applied on populations submitted to random modifications of their structure (for instance different parts of a machine, different groups of a chemical molecule, different parameters of a process, different possible interactions between elements, etc.) before consensus selection (step b)) .

[0014] According to a preferred embodiment, the populations are different genetic sequences (structure of the population) of a group of individuals (such as species, genera, families, etc.), wherein the random modifications of structure are mutations in the genetic sequence and wherein the solutions are philogenic bifurcating trees representing the interactions between said populations. [0015] Preferably, one of the parameters of the consensus selection may require that the bifurcating trees exist in most of or all groups of population. [0016] The present invention relates also to the algorithm for performing the method according to the invention.

[0017] The present invention will be described in detail in the following preferred embodiment of the present invention hereafter called "metapopulation genetic algorithm" which finds suitable solution (s) of a ^"genetic philogenic problem.

Short description of the drawings

[0018] Fig. 1 represents the principle of Consensus Pruning (CP) applied in the method according to the invention.

[0019] Fig. 2 represents the times required for single rounds of perturbation algorithms (NNI , SPR and TBR branch swapping) . [0020] Fig. 3 represents:

— Part a: the score vs. time for genetic algorithm (one population) and metapopulation genetic algorithm

(strict CP with 4 populations of 4 individuals each) runs (80 taxa) ;

- Part b: the score of best tree (320 taxa) for each of four populations run in parallel (no Consensus Pruning) ;

— Part c: the number of consensus portions among solutions from four independent 80-taxa populations running under the GA vs . number of consensus among four 80-taxa populations interacting through Consensus Pruning during a metapopulation genetic algorithm run. [0021] Fig. 4 represents: — Part a: run times for a one-population genetic algorithm search as well as for 2; 4, 6, 8 and 16 population metapopulation genetic algorithm searches;

- Part b: run time (160 taxa) vs. percent error for 2, 4, 6, 8 and 16 population metapopulation genetic algorithm searches under probability Consensus Pruning.

Detailed description of a preferred embodiment of the invention

[0022] The present invention relates to a new genetic algorithm named the "metapopulation genetic algorithm" (metaGA) that vastly improves the speed and efficiency with which solutions are found (such that nucleotide sequence data sets incorporating hundreds or thousands of taxa can be analysed in practical computing times) and yields a probability index for each branch. The metaGA procedure has been incorporated into a computer program for philogeny inference, METAPIGA (Philogeny Inference using the metapopulation Genetic Algorithm) . Analyses on simulated and real data sets demonstrated (15) the efficiency of the metaGA procedure .

Material and methods The GA cycle of evaluation-selection-mutation Evaluation

[0023] At the beginning of each generation, each tree is evaluated (without optimisation of branch lengths or other model parameters) and its log likelihood (hereafter called score) is recorded. Minimal updating (MU) is used, i.e., a specific application of the pruning algorithm of Felsenstein (4) , when calculating the score of each tree . The pruning algorithm takes advantage of the fact that likelihood calculation at a particular node depends only on the values at the two nodes connected to it in the direction of the terminal nodes. MU extends this by taking advantage of the fact that perturbation of the tree topology and branch lengths during the GA search (see below) forces recomputation of the likelihood only at nodes along the path from the changed part(s) of the tree to the centre.

Selection [0024] Once all individuals (trees) have been scored, selection is imposed on the population. Four types of selection are implemented: rank, tournament, replacement, and improve. In rank selection, individuals are assigned a probability of leaving an offspring (i.e., a copy of themselves) as a function of their position in a list in which they are ranked by their score. In METAPIGA, a rank selection identical to that described in Lewis (12) is implemented, i.e., in a population of n individuals ranked by their InL, the probability for the ith individual of leaving an offspring to the next generation is equal to (n - i + l) n (n + l)

[0025] In tournament selection, two individuals are drawn randomly from the population of n individuals, and one offspring is produced from the individual with the higher score . Both trees are then placed back into the mating population, and the whole process is repeated until n offspring have been generated.

[0026] In replacement selection, two individuals are drawn randomly from the population of n individuals and two copies of the better individual are returned to the mating pool (parents are discarded) . The process is repeated sn times, where s is the strength of selection, then the offspring population is generated as an exact copy of the postselection parent population. [0027] All three of the above-mentioned selection regimes tolerate the maintenance of poor trees in the evolving populations, an effect that may, in some cases, allow escape from local optima but generally lowers the efficiency with which the GA searches for the optimal tree (s) . The decreased efficiency is especially noticeable toward the end of the search when beneficial mutations are extremely rare and should never be ignored.

[0028] The improve selection method avoids this problem by allowing only those individuals that have scores better than that of the best tree from the previous generation to produce an offspring. Each individual that fails this test is discarded and replaced by a copy of the current best individual . The latter selection scheme greatly reduces the intra-population variability after each selection step. Local optima are avoided, however, through the metapopulation procedures described hereafter. Mutation

[0029] Each individual, with the exception of the current best individual from each population, is then subjected to a single mutation. Besides SPR and NNI , which are performed as described (3) , two other types of topological mutations were implemented: taxa swap and subtree swap (STS) . The former involves the exchange of two randomly selected terminal branches while the latter corresponds to a swap between two small subtrees. As in the GA of Lewis (12) , branch length mutations (BLM) were also allowed: the length of a random branch in the tree is multiplied by a factor drawn from a gamma distribution with a gamma shape parameter a specified by the user and a mean equal to 1. Exploratory analyses indicated that branch length optimisation yields external branch lengths that are very similar to those obtained through topology-constrained neighbour joining (NJ) , both on a NJ topology and a ML topology. Hence, an original operator has been implemented in the metaGA: BLMint, i.e., BLM is allowed to act on internal branches only (whose initial length is set to the average of all NJ-estimated internal branch lengths) while external branch lengths are fixed throughout the search to the NJ-estimated values. This method typically allows for more rapid exploration of search space. [0030] Finally, an enhanced version of the recombination operator implemented in GAML (12) was incorporated. In Lewis' recombination strategy, individuals

(trees) i and are recombined as follows: the subtree defined by a random branch in i is removed and grafted onto a random branch in individual j after pruning from all taxa represented in the subtree. It was found that the use of this recombination procedure in the GA was ineffective.

[0031] This procedure is indeed unlikely to yield a recombined tree with a better score than any of the two parents because it essentially performs SPR on a tree that has already been modified through pruning of some of the taxa.

[0032] Then, the inventors have defined a new recombination operator (RCM) that involves the identification of a consensus region between two trees from different populations (see the main text) , followed by the swapping of the subtrees defined by that consensus region. By definition, this procedure does not require any pruning of taxa because the consensus branch defines, in both trees, the same partitions. Hence, RCM is more comparable to real DNA recombination (with homologous regions being exchanged) .

Dynamic Operators

[0033] In addition to the possibility of assigning to the operators (NNI, SPR, taxa swap, STS, recombination operator, BLM/BLMint) fixed probabilities, the probabilities were allowed to be assigned dynamically. [0034] Under this option, the probability assigned to each operator is relative to the average contribution that operator has made to improving the population, within the last G generations . To prevent an operator from being switched off while it could become efficient later^" in the search, a lower probability bound (specified by the user) is placed on each operator.

Starting Trees

[0035] To allow for variation across populations, initial topologies are either random (each of the

I individuals in each of the P populations is a different random tree) or constructed through procedures called

"jack-knifed" NJ (JNJ) and noisy NJ (NNJ) . JNJ consists of generating all starting trees with the NJ algorithm, but on a different subset of the data set for each population. [0036] Two variations were implemented: nonoverlapping and overlapping JNJ. The former consists of randomly assigning each character to one of the P populations, such that the original data set is eventually divided into P nonoverlapping sets of characters . On the other hand, the overlapping JNJ consists of randomly assigning a proportion p of the characters to each of the populations. Because the process is independent (and performed from the original data set) for each population, the P sets of characters are typically overlapping and each character may be assigned to 0, 1, ..., or N of the populations. The second procedure, NNJ, uses the full data set but with a modified method for joining nodes: at each step, suboptimal nodes (i.e., nodes that are never joined in classical NJ because they exhibit a non-minimal value in the current distance matrix) are joined with a probability proportional to the inverse of their pairwise distance. [0037] The JNJ and NNJ procedures allow the search to start with trees whose scores are much better than those of random trees while still keeping enough variation among the initial populations for the metaGA to be effective in finding the ML topology. Initial internal branch ^" lengths are set either to those specified by the NJ algorithm or an arbitrary value specified by the user.

Software

[0038] METAPIGA is written in the Java programming language and supports Windows, Unix/Linux, and Macintosh operating systems. Java was chosen to allow for cross- platform compatibility and easy implementation of a user- friendly interface. [0039] Preliminary tests suggested that coding the software in Java instead of C++ would not have an appreciable (i.e., >15%) effect on the speed of the software . Resul ts and Discussion Population Size

[0040] For small (<60 taxa) data sets, varying the population size from two to 32 individuals in single population GA runs decreases accuracy but does not significantly change the mean running time. On the other hand, for large data sets (>60 taxa) , increasing the population size increases run time but has no significant effect on accuracy. In short, small populations yield both the quickest and most accurate results, regardless of the size of the data set.

Relative Efficiencies of the Starting Tree Options [0041] Although a random tree will generally exhibit a score vastly worse than that of a StepAdd (7) or a NJ (16) tree, the former constitutes an unbiased starting point for a search, such that different random trees are less likely (on average) to become systematically trapped in the same local optimum. One could therefore argue that the most conservative strategy for obtaining a globally optimal solution would be to repeat the heuristic search beginning with different starting trees (e.g., different random trees) and test whether most of these searches yield the same topology. Obviously, this approach has a drawback: computing time. Preliminary analyses indicate that the JNJ/NNJ starting tree strategies (described above) are both vastly faster than starting with random trees and efficient at avoiding local optima. Although more extensive analyses need to be performed, these methods seem to be effective compromises between speed and efficiency because they allow both for high starting scores and some random variation.

The metaGA Essence of the procedure

[0042] The inventors have designed a family of heuristic search strategies named the "metapopulation genetic algorithm" (metaGA) . This new approach relies on the coexistence of two or more populations interacting in a "metapopulation" setting. Both the number of populations, and the number of individuals per population can be specified by the user.

[0043] As the metaGA involves several parallel searches, a large amount of inter-population variation can be maintained, even when each population is subjected to strong selection pressures. However, the spirit of the metaGA is that the populations are not fully independent as they are forced to cooperate in the search for the optimal solution. Within each population, specific solutions (here, trees) are subjected to evaluation, selection, and mutation events as they would be in a "classical" single-population GA, but all mutational operators are guided through inter- population comparisons. These comparisons allow the metaGA to identify sub-solutions (i.e., portions of solutions) that are optimal (and should not be modified) and portions that still need to be optimised.

[0044] Communication among the populations is defined and controlled through a principle named "Consensus Pruning" (CP, (15)), i.e., consensus information gathered from comparison of the best solution across populations is used to identify the portions of each solution that can be subjected to mutations (Fig. 1) .

[0045] The concept of CP allows the elaboration of many specific inter-population communication procedures. For example, under "random" consensus pruning, each of the P populations is picked in turn and randomly paired with one of the remaining P-l populations. The consensus between the two solutions define the partitions that cannot be affected by topological mutations. For example, in the case of philogeny inference, if a consensus branch defines a partition between groups I and II, no swap between a taxon in group I and a taxon in group II is allowed, while topological changes within group I and within group II are allowed (Fig. 1) .

Principle of Consensus Pruning (CP)

[0046] Before a solution (here, a tree) is subjected to mutation, it is compared to the best solution (s) from other populations; the consensus branches (indicated in bold red) define the portions of the solution that can

(green arrows) and cannot (red arrows) be affected by mutations; i . e . , any operation exchanging a taxon across a consensus branch is prohibited. TXS, taxa swap; SPR, subtree pruning-regrafting.

Stopping rule

[0047] As the overall number of consensus partitions increases with search time, the metaGA procedure provides a convenient stopping rule: the search stops when (i) the best solutions of all populations are identical, or (ii ) when all mutational changes allowed by the latest consensus information have been attempted on the best solutions of all populations and these changes did not improve their scores. Given the high percentage of partitions that are fixed towards the end of the search, mutating to completion is swift even on very large sets of elements (here, large trees) . Efficiency of the metaGA

Simple GA and other classical heuristics versus the metaGA

[0048] Fig. 3 indicates that the time required by a metaGA search (i.e., Consensus Pruning with strict group- consensus and 4 populations) is much less than that required by stepwise addition (StepAdd) or by the classical

GA. The metaGA is up to 800 times faster than stepwise addition.

[0049] It is important to note that (i) the striking improvement in speed of the metaGA and GA over more classical search strategies is not at the expense of lower accuracy (see below) , and (ii) the ratios of StepAdd to metaGA run times increase with increasing number of taxa.

[0050] As the StepAdd algorithm yields a tree which is typically used as the starting point of a hill-climbing search, it is also shown on Fig. 2 the times required for single rounds of perturbation algorithms (NNI , SPR, and TBR branch swapping (7)), i.e., the time required to swap to completion the single best ML tree. [0051] Part a of Fig. 3 shows the relative run times vs . number of taxa for the GA (one population) , metaGA (with strict group-consensus among 4 populations) and stepwise addition algorithm (StepAdd, (7) ) as well as for a single round of NNI, SPR, and TBR swapping (dotted ^""lines) ; all under a complex ML model (HKY + rate heterogeneity) . Standard errors (among ten runs) are indicated for the GA and metaGA only. A typical heuristic run would comprise several rounds of SPR or TBR swapping on a starting tree obtained by stepwise addition. Even if a reasonable starting tree could be obtained instantaneously (which is far to be the case) , the minimum time required to TBR- or SPR-swap that tree would vastly exceed the time required by a full metaGA search. [0052] The efficiency of the metaGA can be attributed to two factors. First, it allows the stopping rule to be reached very quickly: near the end of the search, consensus-pruning-constrained mutations affect a greatly reduced number of subsolutions, hence swapping to completion is much faster than it would be on an unconstrained solution. Second, the consensus information shared among parallel populations allows them to increase their scores and the number of consensus partitions faster than if they were each searching in isolation (Parts b and c of Fig. 3) . Hence, despite the fact that a 4-population metaGA search requires evaluation of four times more solutions each generation than in a single GA search, the former completes the search much faster than the latter (part a of Fig. 3) .

[0053] Part a of Fig. 3 shows the score vs . time for GA (1 population) and metaGA (strict CP with 4 populations of 4 individuals each) runs (80 taxa) . The asterisks indicate when the stopping rule has been reached. [0054] _ Part b of Fig. 3 shows the score of best tree

(320 taxa) for each of four populations run in parallel (no

Consensus Pruning) . At the time indicated by the red arrow, one population (red) was allowed to use the consensus information from the three other populations (black) . This demonstrates that score increases faster when Consensus Pruning (i.e., the essence of the metaGA) is enabled. [0055] Part c of Fig. 3 shows the number of consensus portions among solutions from four independent 80-taxa populations running under the GA vs . number of consensus among four 80-taxa populations interacting through Consensus Pruning during a metaGA run. Accuracy of the metaGA and effect of the number of populations

[0056] Fig. 4 indicates that computing time increases with the number of populations involved in "Consensus Pruning" . As many as 10 populations are required to slow down the metaGA to computing times similar to those of single-population GA runs.

[0057] Part a of Fig. 4 shows the run times for a one-population GA search (dotted circle; ±SE indicated) as well as for 2-, 4-, 6-, 8-, and 16-population metaGA searches (10 runs, SE too small for being visible) . The vertical arrow indicates the number of populations for which a metaGA run takes the same time than a one- population GA search. Run time increases polynomially with the number of populations.

[0058] Part b of Fig. 4 indicates the run time (160 taxa) vs . percent error for 2-, 4-, 6-, 8-, and 16- population metaGA searches under probability Consensus Pruning. Coordinates of the one-population GA run are indicated by the blue cross. The difference of speed between the GA and the metaGA are much larger with more complex ML models.

[0059] However, as part b of Fig. 4 indicates, two- population metaGA searches are not optimal as they yield trees with 2 to 4% error (in comparison with the best score) . In this case, uninterrupted interaction between the same two populations does not allow consensus pruning to differentiate between informative and random consensus, such that the populations rapidly converge towards each other, terminating the search prematurely. This negative effect decreases dramatically with 3 populations and virtually disappears for runs utilising >3 populations: the likelihood that all populations evolve the same suboptimal portions of solution decreases very rapidly with the number of populations. In all analyses performed so far, -population metaPIGA runs seem to give the best compromise between computing speed and accuracy of topological inference. The optimal combination of maximum speed and accuracy is "probability consensus" with 4 populations.

REFERENCES

I. Cook, S. (1971) in Proceedings of the 3rd Annual Symposium on the Theory of Computing, ACM, New York) , pp. 151-158. 2. Hillis, D. M. (1996) Nature 383, 130-131.

3. Huelsenbeck, J. P. (1995) Syst . Biol . 44, 17-48.

4. Felsenstein, J. (1981) J. Mol . Evol . 17, 368-376.

5. Swofford, D. L. (2002), (Sinauer Associates, Sunderland, Massachusetts) . 6. Felsenstein, J. (1993), (Distributed by the author. Department of Genetics, University of Washington., Seattle) .

7. Swofford, D. L., Olsen, G. J. , Waddell, P. J. & Hillis, D. M. (1996) in Molecular systematics, eds . Hillis, D. M. , Moritz, C. & Mable, B. K. (Sinauer Associates, Sunderland, Massachussetts) .

8. Strimmer, K. & von Haeseler, A. (1996) Molecular Biology and Evolu tion 13, 964-969.

9. Salter, L. A. & Pearl, D. K. (2001) Syst . Biol . 50, 7- 17.

10. Huelsenbeck, J. P., Ronquist, F. R. , Nielsen, R. & Bollback, J. P. (2001) Science 294, 2310-2314.

II. Matsuda, H. (1996) in Pacific symposium on biocomputing

' 96, eds. Hunter, L. & Klein, T. E. (World Scientific, London), pp. 512-523.

12. Lewis, P. O. (1998) Molecular Biology and Evolution 15, 277-283.

13. Katoh, K. , Kuma, K.-I. & Miyata, T. (2001) J.Mol . Evol . 53, 477-484. 14. Foster, J. A. (2001) Nature Reviews in Genetics 2 , 428- 436.

15. Lemmon, A. R. & Milinkovitch, M. C. (2002) submi tted.

16. Saitou, N. & Nei, M. (1987) Mol . Biol . Evol . 4, 406- 425.

Claims

1. Method for obtaining an efficient solution to a combinatorial optimisation problem, comprising the following steps : a) setting p populations (1, 2, 3, ...) of n individuals (e.g., random individuals), each individual being any specific solution (good or bad) to the problem at hand ; computing the quality of each of the n*p solutions,

b) applying specific selection parameters to individuals within each population, for example: keeping only the single best individual within each of the p populations; c) identifying (by a process called "Consensus Pruning" ) the common features among the p best individuals (one individual from each of the p populations) ; d) applying specific mutation parameters to individuals within each population but constraining (strictly or not) the common features not to be mutated; e) repeating steps b) to e) until all the p best individuals are identical or all mutation have been tried on all unconstrained parameters; f) keeping the p solutions (identical or not) as excellent solutions to the optimization problem at hand.

2. The method according to claim 1, wherein the populations are submitted to random modifications of their structure and are allowed to interact by "Consensus Pruning" .

3. The method according to claim 2, wherein the individuals within each population are different different solutions (to an industrial combinatorial optimisation problem) such as a specific transport or telecommunication networks, a specific production process, a specific management or logistic or production decision, a specific engine or any machine or molecule to be designed/optimized, etc.

4. An algorithm and software for performing the method according to any one of the preceding claims .