Viral systems : a new bio-inspired optimisation approach

1 Viral Systems: a new bio-inspired optimisation approach Pablo Cortés*, José M. García, Jesús Muñuzuri and Luis Onieva Seville University Ingeniería Organización. Escuela Superior Ingenieros, Camino de los Descubrimientos s/n. Sevilla 41092. SPAIN Abstract.- The paper presents a new approach to deal with combinatorial problems. It makes use of a biological analogy inspired by the performance of viruses. The replication mechanism, as well as the hosts’ infection processes is used to generate a metaheuristic that allows the obtention of valuable results. The Viral System (VS) theoretical context is described and it is applied to a library of medium-to-large-sized cases of the Steiner problem for which the optimal solution is known. The method is compared with the metaheuristics that have provided the best results for the Steiner problem. The Viral System provides better solutions than genetic algorithms and certain tabu search approaches. For the most sophisticated tabu search approaches (the best metaheuristic approximations to the Steiner problem solution) VS provides solutions of similar quality. Keywords: virus, optimisation, metaheuristic, artificial intelligence, Steiner tree 1 Introduction In real life, most problems are combinatorial. For this type of problem, the available algorithms usually present weaknesses. The mathematically rigorous algorithms (as branching or cutting techniques) tend to be too slow while ad hoc heuristics often produce poor solutions. Since the last decades of the past century, metaheuristics and Artificial Intelligence (AI) have been trying to deal with these difficulties. This is the case of Simulated Annealing (see Kirkpatrick [1] and Rumelhart and McClelland [2] for the first references), Genetic Algorithms (firstly described by Holland [3], Holland et al. [4] and Goldberg [5]), Tabu Search (defined by Glover [6] and [7]) or more recently the use of Agent Technology (see Van Dyke [8], and Jennings et al. [9] for a survey including wide bibliography in the field). The use of agents to deal with optimisation problems is a new and incipient framework. The Network of Excellence “AgentLink” funded by the European Commission under its Fifth Framework Information Society Technologies programme has edited a report describing the current state- * Corresponding author. Fax: +34 95 448 73 29. E-mail: [email protected] 2 of-the-art of agent technologies and identifying trends and challenges that will need to be addressed over the next 10 years to progress in the field and realise the benefits (Luck et al. [10]). Recently, Artificial Immune Systems (AIS) have arisen. They are computer algorithms inspired by the principles and processes of the vertebrates immune system. The algorithms typically exploit the immune system's characteristics of learning and memory to solve a problem. They are coupled to artificial intelligence and closely related to genetic algorithms. Processes simulated in AlS include pattern recognition, hypermutation and clonal selection for B-cells, negative selection of T-cells, affinity maturation and immune network theory. AIS began in the eighties with Farmer et al. [11]. In this paper, we propose the use of a biological analogy based on viral infection. Not much research has been carried out using viral analogies. We have been able to find very few papers proposing the use of viruses , all of them as part of genetic algorithms. For instance, Kubota et al. [12] propose them as part of a specific operator in genetic algorithms, and Saito [13] has described the use of genetic algorithms which make use of a virus evolutionary theory (GAV), and an algorithm based on the conception of horizontal evolution caused by virus infections. GAV is carried out by attacking a chromosome by a number of viruses, and having the genes of the chromosome recombined by the attack. The infection is allowed when the evaluation value goes up, but it falls into local minima easily. In order to escape from these local minima, an infection which makes the evaluation value worse in a small rate under small probability is allowed as well. All these approaches do not fit with our definition of Viral System as a new metaheuristic. Here, we introduce a new optimisation approach that makes use of several ideas from multi-agent systems (MAS), as well as from other well-known AI approaches. In our proposal, we consider that the viruses are part of a general infection, where each virus tries to behave to its own benefit, but resulting in the benefit of the Viral System. A generic optimisation problem can be defined as (1), and its complexity is typically NP-complete (Garey and Johnson [14]). {} nixgxf i,,1 , 0)(:)(Min =∀≤ ( 1 ) Where x is the set of feasible solutions of the problem, f(x) is the objective function, and gi(x) ≤ 0 are the problem constraints. The objective function and the constraints are not necessarily linear, and the variables could be integer or continuous. 3 Based on this problem, the next section of the paper follows with the description of the biological analogy that we will use to solve the problem in section 2, the description of the approach in section 3, and the results for a collection of trial problems in section 4. Finally, we review the main aspects in the conclusion section. 2 Biological analogy description 2.1 Generic definition of virus Viruses are intracellular parasites shaped by nucleic acids, such as DNA or RNA, and proteins. The protein generates a capsule, called a capsid, where the nucleic acid is located. The capsid plus the nucleic acid shape the nucleus-capsid, which defines the virus. Figure 1 depicts different types of viruses. FIGURE 1 Each type of virus has different capabilities with respect to infecting cells, including the mechanisms of replication and spread/inoculation, and strategies to weaken the immune response of the host. Furthermore, each one of the different classes of viruses is able to be lodged in hosts that can be of animal, vegetal or bacterial type. Hence, they can represent unicellular or multi-cellular organisms. 2.2 Virus replication mechanism Although the replication mechanism of viruses depends on the type of virus, we have considered the replication mechanism for phages, which are viruses that infect bacteria. The phage life cycle (see figure 2) follows the following process: 1. The virus is adhered to the border of the bacterium. Subsequently the virus manages to penetrate the border being injected inside this one, (a) and (b) in figure 2. 2. The infected cell stops the production of its proteins, beginning to produce the phage proteins. So, it starts to replicate copies of the virus nucleus-capsids, (c) and (d) in figure 2. 3. After replicating a number of nucleus-capsids, the bacterium border is broken, and the new viruses are released, (e) in figure 2, which can infect other close cells, (f) in figure 2. FIGURE 2 The life cycle of the viruses can be developed in one step, through this lytic replication, or it can include more than one step. Some viruses are capable of lodging in cells giving rise to the lysogenic replication. In this case the process (see figure 3) is as follows: 1. The virus infects the host cell, being lodged in its genome, (a) and (b) in figure 3. 2. The virus remains hidden inside the cell during a while until it is activated by any cause, for example ultraviolet irradiation or X-rays, (c) in figure 3. 4 3. The replication of cells altered, with proteins from the virus, starts. FIGURE 3 The virus lodged in the genome of the host cell can cause alterations in the cell genome, (e) in figure 3. Occasionally certain genes participating in growth processes as the protooncogenes can give birth to tumour cells. Note that the virus enters the cell and then, the infected cell (bacteria for phagocyte case) is known as the organism’s host. 2.3 Antigenic properties in certain viruses and infection end Some viruses have the property of leading an antigenic response in the infected organism. In these situations an immune response is originated causing the creation of antibodies. The end of the virus infection can be reached in two ways: the organism beats the virus implying the host recovery, or the virus beats the defence capabilities of the organism and the host death takes place. 2.4 Natural optimization in Viral Systems and its computational equivalency In nature, viruses are active organisms that reach their major success when infecting individuals with a low level of health (note that the objective function f(x) represents the cell’s health). The viral infection acts in an efficient way by means of developing a continuous process in time searching for individuals, or bacteria in the case of phagocyte infections, that are “suitable for infection”, precisely those with a lower level of health (or low value of f(x)in computational terms). This way, viruses optimize their objective, and infecting less healthy individuals (poorer solutions) viruses are propagated through a major size of the population. In fact, in some occasions when the virus infection capabilities decrease, the virus using the cell genome (solution encoding) mutates to a different form with the intention of achieving higher success in the infection process, see figure 3. This describes a continuous process of the infection’s systematic progress lodging on unhealthy cells that can lead the organism up to two different ends. The organism dies (it corresponds to reaching the optimum or beating the gap) or the organism resists the infection according to the virus isolation case (it corresponds to reaching the maximum number of iterations). In fact, the concept is close to reality because “more successful epidemics” are those where a larger number of deaths took place. For example, this was the case of the Black Death killing between a third and two thirds of Europe's population in the mid-late 14th century, the Spanish Flu Pandemic causing between 50 million and 100 million people worldwide killed in 1918 and 1919, or more recently the African Ebola with mortality rates ranging from 50% to 90%. So, following this data the death of the organism would lead to the objective (optimum) of a “successful infection”. 5 3 Viral System description 3.1 Viral System components The use of natural analogies is quite common when dealing with optimization in NP-complete problems. In this line, multiagent systems (MAS) have contributed a large number of methods taken from the natural world. This is the case of ants for path planning or brood sorting, termites for nest building, wasps for task differentiation, birds and fish for flocking, and wolves for surrounding preys (see Van Dyke [8] for an in-depth explanation). Here, we present viruses infecting organisms. Our Viral System (VS) is an approach capable of dealing with any optimisation problem. To present the VS we follow the Van Dyke suggestion. Each VS is defined by three components: a set of viruses, an organism and an interaction between them: VS = <Virus, Organism, Interaction> The Virus component of the VS is a set consisting of single viruses: Virus = { Virus1, Virus2 , … , Virusn } Each virus is defined in four components: Virusi = <Statei, Inputi, Outputi, Processi> Where each component means:  Statei characterizes the virus. It defines the cell infected by the virus. It is typically the mathematical encoding of the solution in computational terms, which we also call genome. A concrete virus, Virusi, can produce the infection of a cell of the organism providing a host. Additionally, the evolution of the residence time of the virus inside the cell can be defined by the number of nucleus-capsids replicated for the lytic replication (NR) or the number of iterations for the lysogenic replication (IT). So, the three-tuple genome-NR-IT defines the Statei for the Virusi.  Inputi identifies the information that the virus can collect from the organism. This information is always collected in the proximity of the virus. Inputi represents the input’s interaction with the organism (organism’s information Æ virus). It corresponds to the neighbourhood of the cell in computational terms.  Outputi identifies the actions that the virus can take. Outputi represents the output’s interaction with the organism (virus Æ organism). It corresponds to the selection mechanism of the type of virus replication in computational terms. 6  Processi represents the autonomous behaviour of the virus, changing the Statei. It corresponds to the replication operator process in computational terms. The Organism component of the VS is defined by two components: Organism = < Stateo, Processo> Where each component means:  Stateo characterizes the organism state in each instant. It consists of the clinical picture and the lowest healthy cell (the best solution found of the optimization problem). The set of feasible solutions in a specific space ℜn is given by the problem constraints (2). {} nixgx i,,1 , 0)(:K =∀≤= ( 2 ) Where each one of the feasible solutions of the problem (2), x∈K, has been called a cell. The genome is the mathematical encoding of each cell or feasible solution. When a virus infects a cell, this cell enters the population of infected cells. The total amount of infected cells constitutes the infected part of K for each time instant, and it is named “clinical picture”. It contains the overall information of the infection needed by the algorithm in each instant, t. Thus, the clinical picture consists of every three-tuple genome-NR-IT defining the Statei of each Virusi. In the same way, the overall clinical picture plus the lowest value of f(x) defines the Organism Stateo. Figure 4 depicts the State concept for the organism and the viruses. FIGURE 4  Processo represents the autonomous behaviour of the organism that tries to protect itself from the infection threat, consisting of antigen liberation. Medically, an antigen is any substance that elicits an immune response. The antigens generate an immune response by means of antibodies trying to fight the virus infection. The computational mission of the antigens is to liberate space in the population of infected cells (clinical picture), trying to maintain free record memory in the clinical picture to incorporate new infected cells (new feasible solutions). Thus, due to the antigens’ activity, infected cells (in the clinical picture) can be recovered (removed) and cells in the organism that could be infected are not infected due to this antigenic substance. In computational terms, the antigens can be viewed as destructive agents, as described in Talukdar et al. [15]. 7 3.2 Viral System Interaction The Interaction component of the VS is conditioned by the Input and Output actions that lead to a Process of every virus and the corresponding Organism response. A Virusi process implies a resulting change in the organism, and the same applies for an Organism’s process. The interaction is the union of both actions. 3.2.1 Virus input sensor: Neighbourhood identification The input sensor of each virus, Inputi, collects information from the organism. The sensors map the genome of the cell and detect the set of cells close to the infected one. This set is named the neighbourhood of the feasible solution x, V(x). The neighbourhood depends on the shape of the constraints of the problem, gi(x). 3.2.2 Virus output ejector: replication type selection The ejector, Outputi, selects the type of evolution of the virus. We consider two options: the lytic replication if the life cycle of the virus is executed in only one step; and the lysogenic replication if the life cycle of the virus is executed in two steps. So, in order to select the type of replication, we consider lytic replication (with probability plt) or lysogenic replication (with probability plg), where plt + plg = 1. 3.2.3 Process: Lytic replication For the lytic replication case we have considered two different evolutions for different types of infection: a parsimonious infection, selecting the cells to be infected; and a devastating infection where a massive number of cells is infected. An example of the first case is the HIV virus, which through a step-by-step evolution destroys the immune system during a process that can take years. On the contrary, an example of the second case is the Ebola virus with a rapid and massive infection that very often produces the death of the patient in a few days. 3.2.3.1 Virus process in massive infection The lytic replication starts only after a specific number of nucleus-capsids have been replicated. So, each time instant (iteration t) a number of virus replications (NR) takes place. The number of replications per iteration is calculated as function of a binomial variable, Z, adding its value to the total NR. After a specific number of nucleus-capsids has been replicated inside the cell (LNR), the bacterium border is broken, liberating the lodged viruses. All these viruses are active and prepared to infect new cells. The value of LNR depends on the cell’s health conditions. So a healthy cell (with high value of f(x)) will have low probability of getting infected, and therefore the value of LNR will be higher. In the opposite it will have a lower value of LNR. The following equation (3) shows the calculation procedure for LNR in a cell x: 8 () () () () LNRfor valueinitial theis LNR analysed. being cell infected theis and ) of (in terms problem theofresult known best theproduces that cell theis where LNRLNR 0 0 -cell x xfx xf xfxf x    ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛− ⋅= ( 3 ) The number of nucleus-capsids replicated each iteration can be approximated by a Binomial distribution given by the maximum level of nucleus-capsids replicated, LNR, and the single probability of one replication, pr,: Z = Bin (LNR , pr). Once the distribution has been stated, we can calculate the probability of replicating exactly z nucleus-capsids, P(Z=z), as well as the average, E(Z), and variance, Var(Z), equations (4-6). () z r z rpp z z− − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ == LNR 1· LNR )Z(P ( 4 ) LNR·)Z(E r p= ( 5 ) () ·LNR1·)Z(Var rr pp −= ( 6 ) Once the bacterium border is broken liberating the viruses, each one of the viruses has a probability, pi, of infecting other new cells of the neighbourhood. Let the neighbourhood cardinality of the feasible solution, x, be |V(x)|. Then, we can define Y as a binomial random variable representing the cells infected by the virus in the neighbourhood. The parameters of the binomial distribution are Y = Bin (|V(x)|, pi). In this situation, the probability of infecting exactly y nucleus-capsids, P(Y=y), as well as the expected number of new infected cells, E(Y), and the variance, Var(Y), will be (7-9): () yxV i y ipp y xV y− − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ == |)(| 1· |)(| )Y(P ( 7 ) |)(·|)Y(E xVpi = ( 8 ) () |)(|·1·)Y(Var xVpp ii −= ( 9 ) Under these conditions, the growth of the total number of infected cells follows a geometric progression and there is no apparent limit to the growth of the population infected. So it must be bounded by another 9 phenomenon in order to maintain the VS bounded-in-computation (this is described in the Organism Process section as the antigenic response). 3.2.3.2 Organism process in massive infection The organism responds to the viral infection by releasing antigens. The antigens generate an immune response by means of antibodies. The Organism Process is defined by the antigenic capacity of the organism, which depends on the Stateo and leads to a new situation of the Organism State. Each one of the infected cells in the clinical picture has a probability of developing antibodies against the infection. The probability of generating antibodies for any cell of the infected population can be given by a Bernoulli probability distribution. Therefore, the probability of generating antibodies is pan: A(x) = Ber (pan), where x is an infected cell. Hence, the total population of infected cells generating antibodies is characterized by a Binomial distribution of parameters: the size of the population infected, n, and the probability of generating antibodies, p an: A(population) = Bin (n, pan). On the other hand, we know that when a virus reaches the LNR limit, the bacterium border is broken and new viruses try to infect new cells in their neighbourhoods. Under these conditions, we must compute the antigenic capacity for every cell in the neighbourhood of an active virus. It is estimated as a Bernoulli probability distribution given by the probability of generating antibodies, pan: A(x’) = Ber (pan) : x’∈V(x). Therefore, the total number of cells in the neighbourhood with antibodies will follow a Binomial probability distribution given by the total size of the neighbourhood for all the active viruses, |V(x)|, and the probability of generating antibodies, pan: A = Bin (|V(x)|, pan). Then, the probability of finding exactly a immune cells, P(A=a), as well as the expected value of immune cells, E(A), and the variance, Var(A), is given by (10-12). () () a an a an pp a x a− − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ == V 1· V )A(P ( 10 ) () xpan V·)A(E = ( 11 ) ()() xpp anan V·1·)A(Var −= ( 12 ) For this massive infection case, we have to guarantee two opposed tendencies. Firstly, we have to allow enough capability for the infection expansion (infecting as many cells in the neighbourhood as possible). Secondly, we have to keep the infection evolution bounded in order to have enough gaps in the clinical picture to guarantee 16 If CP(c).NR > CP(c).LNR Then Interaction_Virus(c) End If End Procedure Procedure Interaction_Virus(c) Get_ Neighbourhood (c , Vc) {Get the probability of generating antibodies} A = Get_Rnd_Binomial_Probability(|Vc| , pan) a = Get_Value_Binomial(A, |Vc|, pan) {Get the probability of infection in the neighbourhood} Y = Get_Rnd_Binomial_Probability(|Vc| , pi) y = Get_Value_Binomial(Y , |Vc| , pi) Get_Best_Neighbourhood_Solution(Vc , y , a , cNEW) c= cNEW Replicat_Type(CP(c)) = Get_Rnd_Replication_Type (plt) End Procedure Procedure Get_Best_Neighbourhood_Solution(Vc , y , a) j = Get_Rnd_Starting_Evaluation(|Vc|) i = 0 Do i = i + 1 Solution = Get_Solution(Vc , j) if health(cNEW ) < health(Solution) then cNEW = Solution end if j = j+1 Loop Until i = y-a End Procedure Procedure Lysogenic(c) CP(c).IT = CP(c).IT + 1 If CP(c).IT > CP(c).LIT Then Mutation(c, cNEW) c = cNEW Replicat_Type(CP(c)) = Get_Rnd_Replication _Type(PLytic) End If End Procedure 4 An application example for Viral Systems: the Steiner problem We have used the Steiner problem, a well-known NP-Complete problem, to test the efficiency of our Viral System described above. The Steiner problem can be stated as follows:  Given a non-directed graph G = (N,A) with |N| nodes and |A| arcs with costs cij ∀(i,j)∈A; and a subset T⊆N with |T| nodes called terminals or targets, with the rest of the nodes in N called Steiner nodes, 17  find a network GT ⊆ G joining all the terminal nodes in T at minimum cost. This network can include some of the Steiner nodes but does not have to include all the Steiner nodes. The Steiner problem has been widely dealt with in the scientific literature. There exists a wide range of heuristics giving good approximations (e.g. Takahashi and Matsuyama [16], and Beasley [17]), as well as other well-proven direct methods like branch and bound, branch and cut, cutting plane algorithms, etc. (e.g. Wong [18], Chopra et al. [19], and Grötschel et al. [20]). Also, metaheuristics have been applied to the Steiner problem successfully. Gendreau et al. [21] for Tabu Search, as well as the Genetic Algorithms from Esbensen [22], and Voss and Gutenschwager [23] are three of the most notable approaches. 4.1 VS characterization for the Steiner problem The Organism state is depicted by the clinical picture representing the infected part of the Steiner problem hull, K, containing all the feasible solutions. A common formulation (21) for K is the coverage one, e.g. Koch et al. [24]. ()() ( ) ( ) () integer , , 10 \ , , , 1δ: x Ajix TWNTWNWWXK ij ∈∀≤≤ ∅ ≠ ∅ ≠ ⊂∀≥ ∩∩ ( 21 ) Where δ(X) denotes the cut induced by X⊆N, that is, the set of arcs with one node in W and one in its complement, and () AFxFX Fji ij ⊆∀= ∑ ∈ , )( , . It is easy to see that there is a one-to-one correspondence between Steiner trees in G = (N,A) and {0,1} vectors satisfying K. We represent the genome of the cells by a bit string of size equal to |N| in which each bit position i corresponds to the node i in the graph. A 1 means that the node i is connected, while the bit is set to 0 otherwise. As all the terminals must be in the Steiner tree, it is sufficient to use a bit string of size |N-T| including only the Steiner nodes belonging to the Steiner tree. Under these conditions, a Steiner tree can be constructed by means of a minimum spanning tree (MST) that contains all the terminal nodes (set T), the subset of Steiner nodes in the bit string fixed to 1 and, perhaps, some artificial arcs if the set is disconnected. In case of necessity of introduction of artificial arcs due to the disconnection of the tree, there will be diverse possibilities. We make use of the graph construction mechanisms described in Gendreau et al. [21]. 18 Once we have stated the cell genome we can define the Virus state. The three-tuple formed by the genome of each cell infected plus the number of replicated nucleus-capsids (in the case of lytic replication) or the number of generations (in the case of lysogenic replication) defines the virus state. The entire infected cell population, which is the clinical picture, and the best solution complete the Organism and therefore the Virus state. The Output ejectors of the Virus component of the VS are clearly defined by the type of replication. On the contrary, the Input sensors must be carefully stated. In fact, a key decision in every VS is to state an adequate cell neighbourhood for the virus in the lytic replication process and a genome alteration process for the lysogenic replication. In the case of the Steiner problem, for the lysogenic replication and given a feasible solution x∈K, the genome alteration is made by flipping a bit in the string. For the lytic replication and given a feasible solution x∈K, the neighbourhood of such a solution consists of the set of bit strings that can be obtained by the removal or the addition of a single Steiner node from/to the current cell encoding. In order to be efficient, the new MSTs must be found by manipulating a rooted tree data structure carefully. Gendreau et al. [21] describes the procedure in-depth. Finally, the Virus and Organism components are completed by the specification of the Process. The Organism Process consists of the antigenic response and it is mainly determined by the determination of the parameter pan , as was stated in section 3.2.3.4. The Virus Process consists of the determination of the type of replication that is conditioned by the parameters plt and plg. Additionally, the Virus Process depends on the parameters of replication, pr, infection, pi, and the limits LNR0 and LIT0. Due to the special encoding for the Steiner problem solutions the neighbourhood size is constant and equal to the number of Steiner nodes. It must be noted that the neighbourhood is set by changing the value of each bit from 0 to 1 and vice versa. The Interaction takes place after the selection of the Virus Process. It is basically determined by the random evolution of the viral infection and the antigenic capacity of response. 4.2 Specific pseudocode procedures for the Steiner problem The general pseudocode functions and procedures were described in section 3.4. We have included here the more specific procedures for the Steiner problem. They are mainly the neighbourhood characterization and the Steiner-oriented mutation for the lysogenic replication. Such procedures will have to be specifically oriented to the problem being analysed. The following pseudocode describes them. 19 Procedure Get_ Neighbourhood (c, Vc) Vc = ∅ {Neighbourhood} For i = 1 to NumSteiners If i∈Solution(c) then NeighbourNEW = Extract_from_Solution(c, i) Else NeighbourNEW = Introduce_in_the_Solution(c, i) End if Vc = Vc + NeighbourNEW Next End Procedure Procedure Get_Mutation_Solution(c, cNEW) SNode = Get_Rnd_SteinerNode(NumSteiners) If SNode∈Solution(c) then cNEW = Extract_from_Solution(c, SNode) Else cNEW = Introduce_in_the_Solution(c, SNode) End if End Procedure 4.3 Computational results In order to evaluate the quality of the VS solutions, we have compared the VS results with the optimal solutions of a test problems’ collection of the OR-Library (series C, D and E, each one of them including 20 problems). Steiner problems series C consists of trials with 500 nodes, a number of arcs varying from 625 to 12,500, and terminals from 5 to 250, while series D consists of problems with 1,000 nodes, arcs varying from 1,250 to 25,000, and terminals from 5 to 500; and finally series E includes trials of 2,500 nodes, arcs varying from 3,125 to 62,500, and terminals from 5 to 1,250. The collection is accessible in the Web page http://people.brunel.ac.uk/~mastjjb/jeb/info.html. For the statistical calculations, due to the large size of the problems, we were forced to replace the binomial probabilistic distribution by its corresponding normal approximation. Firstly, Steiner OR-Library problems were pre-processed according to the graph reduction rules (Winter [25]) that reduce the graph size. The graph characteristics are defined by the number of nodes, terminals and arcs. One of the parameters that determine the difficulty of obtaining the optimum is the percentage of terminals with respect to the total number of nodes, as Cortés et al. [26] have previously stated, although other factors like the total number of nodes and arcs also affect the computing time. The argument is based on the fact that problems with a low density of terminals or with a high density of terminals have a small feasible region; however problems with a medium density of terminals have a large feasible region. 20 Using this rationale, we divided the problems into three groups: group 1 for problems with less than 15% of terminal nodes; group 2 for problems with a percentage of terminal nodes between 15% and 30%; and group 3 for problems with more or equal than 30% of terminal nodes. According to the special characteristics of the Steiner problem we realised that the selective infection provided a more adequate approach. We used an experimental design based on a fractional factorial design at two levels in order to analyse the impact of the VS parameters on solution quality. The fractional factorial design is a factorial experiment in which only an adequately chosen fraction of the treatment combinations required for the complete factorial experiment is selected to be run. In fact, for a number of factors k in a full factorial design, 2k runs must be computed, which can quickly become very large. The solution to this problem is to use only a fraction of the runs specified by the full factorial design. The advantages of choosing fractional factorial designs for 2-level experiments are the desirable properties of being both balanced and orthogonal. See, for example, Montgomery [27] for an in-depth explanation on fractional factorial design. The main two-level parameters of the VS according to the selective infection were encoded for the fractional factorial design as:  ITER: number of iterations. Fractional factorial design levels: 50,000 (+) / 10,000 (-)  POB: clinical picture. Fractional factorial design levels: 100 (+) / 50 (-)  PLITI: probability for the lytic cycle. Fractional factorial design levels: 0.7 (+) / 0.3 (-)  LNR: Limit of iterations for the lytic cycle (number of replications inside the cell to be broken). Fractional factorial design levels: 20 (+) / 10 (-)  LIT: Limit of iterations for the lysogenic cycle. Fractional factorial design levels: 20 (+) / 10 (-)  PZ: Probability of replicating z nucleus-capsids. Fractional factorial design levels: 0.7 (+) / 0.3 (-) The experimental design was done with the Stein-C problems from the OR-Library. The fractional factorial design at two levels generated a total number of 26 parameter settings that were executed 4 times in order to reduce the random component. This gave a total amount of 256 trial parameter settings for each Stein-C problem. Since there are 20 test problems in the Stein-C set, we solved 5,120 instances. We grouped the parameter settings in the Stein-C set into three different categories according to the percentage of terminals (a critical parameter when solving the Steiner problem). Thus, Stein-C group 1 includes parameter settings with percentage of terminals under 15%, Stein-C group 2 parameter settings with percentage of terminals between 15% and 30%, and Stein-C group 3 for parameter settings with percentage of terminals higher than 30%. Hence, the group 1 is formed by steinc1, c2, c6, c7, c11, c12, c16, and c17 cases (giving 2,048 instances), while the 21 group 2 is formed by steinc8, c9, c13, c14, c18, and c19 cases (giving 1,536 instances), and the group 3 is formed by steinc3, c4, c5, c10, c15, and c20 cases (giving 1,536 instances). We used Limdep econometric software in order to adjust the mode of the experimental design. In accordance with the selected codification, the variance analysis coincides with the minimum square error regression analysis. The results summarised for the three groups are shown below in Table 2, TABLE 2 We have used italics to highlight the significant parameters for each specific group. Group 1 shows that the number of iterations (ITER), the clinical picture size or population (POB), the lytic replication probability (PLITI), and the limit of iterations for the lysogenic cycle (LIT) are the most relevant parameters in the experimental design. Also, the interactions with a t-ratio higher than 1 or lower than -1 must be taken into account. For example, the independent analysis of parameter ITER will recommend a low value of the parameter, whereas in the case of parameter POB a high value is recommended; however, a positive combination of both parameters ITER and POB (which is also called the combined effect or the cocktail) is recommended in order to obtain a better fitting. This situation corresponds to a contribution ‘+ +’ or ‘- -’ for ITER_POB. We have to take into account all these possible combinations in order to select the best overall combination of parameters (the final selection is shown in Table 3). Finally, we highlight that we obtained a very good adjusted R-squared equal to 92.1% for the data of this group. Group 2 shows that the number of iterations (ITER), the clinical picture sizes (POB), the lytic replication probability (PLITI), and the limit of iterations for the lysogenic cycle (LIT) are the most relevant parameters in the experimental design. Also, the interactions with a t-ratio greater than 1 or lower than -1 must be taken into account. We obtained a very good adjusted R-squared equal to 93.1%. Finally, for group 3, the number of iterations (ITER), the clinical picture size (POB), the lytic replication probability (PLITI), the limit of iterations for the lytic cycle, and the limit of iterations for the lysogenic cycle (LIT) are the most relevant parameters in the experimental design. Also, the interactions with a t-ratio greater than 1 or lower than -1 were taken into account. We obtained a very good adjusted R-squared equal to 98.4%. Table 3 describes the best parameter selection for each group of Steiner problems. The selection of parameters was done by calculating the average error for each possible combination of sequences. TABLE 3 It is important to note that the VS efficiency is non-dependent on the probability of generating a great or low number of nucleus-capsids (parameter PZ) with respect to the Steiner problem case, so its performance showed 22 non-dependency from this parameter. The rest of parameters can affect in different ways the objective of reducing the error of the method. After executing the whole fractional factorial design at two levels we obtained the optimum for all the Stein-C set problems. Then we used the parameters from Table 3 to solve the more complex Stein-D and Stein-E collections. First, we executed four times the VS with the first set of parameters, and then we executed four times the VS with the second set of parameters. Table 4, Table 5 and Table 6 show the results (in error percentage with respect to the optimum) for the Stein-C, Stein-D and Stein-E problems and their comparison with the Tabu Search approach from Gendreau et al. [21] (their P-Tabu and F-Tabu methods), which is the best approach for the Steiner problem in terms of solution quality. The tables include the best Genetic Algorithm approaches that have solved the Steiner problem, which are due to Esbensen [22] (GA-E) and Voss and Gutenschwager [23] (GA-V). Finally, we also show the results obtained for the Minimum Path Heuristic (MPH) by Takahashi et al. [16] and their results reported by Gendreau et al [21]. This collection of methods provides an adequate basis of comparison in order to test the effectiveness of our VS. TABLE 4, TABLE 5, TABLE 6 According to the results of Table 4, Table 5 and Table 6 (a total of 60 problems), VS was the best approach in 48 times and beat clearly the GA-E (22 times), GA-V (9 times), P-Tabu (40 times) and MPH (25 times) approaches. Only F-Tabu showed better performance, being the best approach 51 times. However, VS provided a better solution for the C9, C14, D3, E3, E12, E15 and E19 problems. VS provided very valuable results taking into account that F-Tabu was processed after selecting the 100 best different trees found by the MPH algorithm, after executing the P-Tabu approach as an initial search and reprocessing it into the final Full Tabu Steiner (F- Tabu). So the quality of the F-Tabu results is very high but it is also very much conditioned by the very good seed that was provided (MPH). Furthermore, it is remarkable that MPH obtains better results than special oriented genetic algorithms such as GA-E and GA-V; it gives an idea of the seed quality for F-Tabu. Undoubtedly it is due to the heuristic being executed 100 times in order to obtain the best possible seed to be supplied to the P-Tabu and F-Tabu methods. On the contrary, we apply VS directly to the reduced graph without pre-processing it with any special previous heuristic as MPH (or the others that F-Tabu uses). Nevertheless, our initial clinical picture (the seed of our method) is wholly random-generated. We did not use a good seed provided by a good heuristic because we were interested in observing the quality of the VS evolution to the final solution, more than on outperforming previous heuristics. However, we realized that without searching for a 23 good seed we were obtaining results equivalent (in quality terms) to the best Steiner approach: the F-Tabu algorithm. We also have to state that when the problem was fully executed for the overall set of parameters with all the possible combinations, that is, what we did with the Stein-C, VS always obtained the optimum of the problem. Nevertheless, when we reduced the possible combinations of parameters (focusing them on Table 3) the results quality was maintained on a very high standard, but several (very few) times we were outperformed by F-Tabu, due to all our previous explanations (mainly the preparation of the seed by means of other previously executed heuristics). With respect to the time consumption, we have omitted the CPU times in Tables 7, 8 and 9 as they cannot be compared with other algorithms whose tests were run in a different computer. However, we can estimate the order of time consumption by the algorithm’s complexity, given in (22). ( ) 2 ~NgraphsNumSteinerITEROtime ⋅⋅ ( 22 ) where ITER is the maximum number of iterations, NumSteiners the number of Steiner nodes in the graph and Ngraph the total number of nodes in the graph. The evaluation of the neighbourhood of a cell is the process that requires the largest amount of computational time in the algorithm, and this process is the same used in the F-Tabu algorithm to reach the next solution in the search procedure. Therefore, the worst case would be to run lytic replications in all iterations of the algorithm. For that case, the degree of complexity would be equal to the number of neighbourhood solutions evaluated multiplied by the complexity of the evaluation process. The number of evaluations is obtained as the number of iterations multiplied by the number of Steiner nodes. The evaluation process is the calculation of a minimum spanning tree, which is O(Ngraph2). Furthermore, many combinations of VS parameters were analysed in our experimentation. Hence, we just indicate in Table 7 the average of the computational times (in CPU-seconds) on the best (ITER=10.000, PLITI=0.3 and LNR=20) and worst cases (ITER=50.000, PLITI=0.7 and LNR=10). TABLE 7 5 Conclusions In this paper we have presented a new approach to optimize combinatorial problems called Viral System due to its inspiration in the infection behaviour of viruses. The method was tested with the Steiner problem, which is NP-Complete. VS was applied to a library of medium-to-large-sized problems for which the optimal solution is 24 known and the results were compared with the best Tabu Search approach (and perhaps the best metaheuristic for the Steiner problem) and the best Genetic Algorithm approaches that we found in the scientific literature. VS clearly improves the results from the Genetic Algorithms (a bio-inspired evolutionary methodology close to our proposal) and for several cases VS obtains better results than the Tabu Search, even taking into account that the Tabu Search starts from a very good seed provided by the MPH heuristic. Our future research is focused on generating an optimisation library associated to different classes of viruses (Figure 1 showed several cases). Here we have only addressed the phagocyte case, which corresponds to the less complex virus infection. We want to test different type of infections with different set of problems trying to identify the most adequate one to each type of problem. 6 Acknowledgements The authors acknowledge two anonymous referees for their helpful comments on earlier versions of this manuscript that have significantly contributed to enrich the quality of the paper in its final form. 7 References [1] Kirkpatrick, S. Optimization by simulated annealing: Quantitative studies. Journal of Statistical Physics 1982; 34: 975-986. [2] Rumelhart, D.E., McClelland, J.L. Parallel Distributed Processing: Explorations in the Microstructure of Cognition, Volume I: Foundations, MIT Press, 1986. [3] Holland, J.H. Adaptation in Natural and Artificial Systems, University of Michigan Press, 1975. [4] Holland, J.H., Holyoak, K.J., Nisbett, R.E., Thagard, P.R. Induction: Processes of Inference, Learning and Discovery, MIT Press, 1986. [5] Goldberg, D.E. Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, 1989. [6] Glover, F. Tabu Search – Part I. ORSA Journal on Computing 1989; 1 (3): 190-206. [7] Glover, F. Tabu Search – Part II. ORSA Journal on Computing 1990; 2 (1): 4-32. [8] Van Dyke Parunak, H. “Go to the ant”: Engineering principles from natural multi-agent systems. Annals of Operations Research 1997; 75: 69-101. [9] Jennings, N.R., Sycara, K., Wooldridge,M. A roadmap of agent research and development. Autonomous Agents and Multi-Agent Systems 1998; 1: 7-38. 25 [10] Luck, M., McBurney, P., Preist, C. Agent Technology: Enabling Next Generation Computing, AgentLink, EC V FP Information and Service Technologies Programme, 2003. [11] Farmer, J.D., Packard, N. and Perelson, A. The immune system, adaptation and machine learning. Physica D 1986; Vol. 22: 187--204. [12] Kubota, N., Fukuda, T., Shimojima, K. Virus-evolutionary genetic algorithm for a self-organizing manufacturing system. Computers and Industrial Engineering 1996; 30 (4): 1015-1026. [13] Saito, S. A Genetic Algorithm by Use of Virus Evolutionary Theory for Combinatorial Problems. Optimization and Optimal Control, Vol.1, World Scientific Publishing Co. 2003; 10: 251-268. [14] Garey, M.R. and Johnson, D.S. Computers and Intractability. W.H. Freeman, New York, 1979. [15] Talukdar, S., Gove, A., De Souza, P. Asynchronous teams: cooperation schemes for autonomous agents. Journal of Heuristics 1998; 4: 295-321. [16] Takahashi, H., Matsuyama, A. An approximate solution for the Steiner problem in graphs. Mathematica Japonica 1980; 24: 573-577. [17] Beasley, J.E. A Heuristic for Euclidean and Rectilinear Steiner Problems. European Journal of Operational Research 1992; 58: 284-292. [18] Wong, R.T. A Dual Ascent Approach for the Steiner Tree Problems on a Directed Graph. Mathematical Programming 1984; 28: 271-287. [19] Chopra, S. Gorres, E.R., Rao, M.R. Solving the Steiner Tree Problem on a Graph Using Branch and Cut. ORSA Journal on Computing 1992; 4 (3): 320-335. [20] Grötschel, M., Martin, A., Weismantel, R. Packing Steiner Trees: A Cutting Plane Algorithm and Computational Results. Technical report SC92-09, Konrad-Zuse-Zentrum für Informationstechnik, Berlín, 1992. [21] Gendreau, M., Larochelle, J.-F., Sansò, B. A tabu search heuristic for the Steiner tree problem. Networks 1999; 34 (2):162-172. [22] Esbensen, H. Computing near-optimal solutions to the Steiner problem in a graph using genetic algorithm. Networks 1995; 26, 173-185. [23] Voss, S. and Gutenschwager, K. A chunking based genetic algorithm for the Steiner tree problem in graphs. In Pardalos, P.M., Du, D., (Eds.) Network Design: Connectivity and Facilities Location, DIMAC series in Discrete Mathematics and Theoretical Computer Science 40, AMS, Providence, 1999. p. 335-355. 32 Table 3. Adequate selection of parameters for each group of problems in the Steiner collection Group 1 Group 2 Group 3 Parameters 1st set 2nd set 1st set 2nd set 1st set 2nd set ITER 50,000 (+) 10,000 (-) 50,000 (+) 50,000 (+) 10,000 (-) 50,000 (+) POB 100 (+) 50 (-) 100 (+) 50 (-) 50 (-) 50 (-) PLITI 0.7 (+) 0.7 (+) 0.7 (+) 0.7 (+) 0.7 (+) 0.7 (+) LNR 15 (~) 15 (~) 15 (~) 15 (~) 10 (-) 10 (-) LIT 10 (-) 10 (-) 20 (+) 10 (-) 10 (-) 10 (-) Pz 0.5 (~) 0.5 (~) 0.5 (~) 0.5 (~) 0.5 (~) 0.5 (~) Table 4. Comparison on the solution quality for the reduced C graphs of the Steiner collection Problem Optimum GA-E GA-V MPH P-Tabu F-Tabu VS C1 106 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% C2 220 1,67% 0,83% 0,00% 0,00% 0,00% 0,00% C3 1565 0,13% 0,13% 0,00% 0,00% 0,00% 0,00% C4 1935 0,11% 0,04% 0,09% 0,00% 0,00% 0,00% C5 3250 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% C6 67 0,73% 1,09% 0,00% 0,00% 0,00% 0,00% C7 103 1,76% 2,75% 0,00% 0,00% 0,00% 0,00% C8 1072 0,63% 0,51% 0,00% 0,00% 0,00% 0,00% C9 1448 1,05% 1,30% 0,99% 0,14% 0,14% 0,00% C10 2110 0,26% 0,27% 0,09% 0,00% 0,00% 0,00% C11 29 1,88% 1,88% 0,00% 0,00% 0,00% 0,00% C12 42 1,30% 0,43% 0,00% 0,00% 0,00% 0,00% C13 500 1,01% 1,32% 0,78% 0,00% 0,00% 0,00% C14 667 0,87% 0,68% 1,24% 0,31% 0,31% 0,00% C15 1116 0,25% 0,22% 0,18% 0,00% 0,00% 0,00% C16 13 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% C17 23 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% C18 223 0,71% 0,71% 5,31% 0,88% 0,00% 0,00% C19 310 0,41% 0,82% 4,79% 0,68% 0,00% 0,00% C20 537 0,00% 0,00% 0,37% 0,00% 0,00% 0,00% Best approach 5 5 11 16 18 20 Table 5. Comparison on the solution quality for the reduced D graphs of the Steiner collection Problem Optimum GA-E GA-V MPH P-Tabu F-Tabu VS D1 106 0,57% 0,88% 0,00% 0,00% 0,00% 0,00% D2 220 0,00% 0,73% 0,00% 0,00% 0,00% 0,00% D3 1565 0,92% 1,25% 0,77% 0,26% 0,06% 0,00% D4 1935 0,52% 0,63% 0,16% 0,00% 0,00% 0,00% D5 3250 0,12% 0,19% 0,06% 0,00% 0,00% 0,00% D6 67 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% D7 103 1,94% 3,62% 0,00% 0,00% 0,00% 0,00% D8 1072 1,55% 2,28% 1,59% 0,47% 0,37% 0,47% D9 1448 0,50% 1,15% 0,83% 0,41% 0,21% 0,69% D10 2110 0,13% 0,44% 0,38% 0,00% 0,00% 0,00% D11 29 2,07% 1,84% 0,00% 0,00% 0,00% 0,00% D12 42 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% D13 500 0,56% 1,48% 2,00% 0,00% 0,00% 0,00% D14 667 0,30% 0,75% 0,75% 0,15% 0,15% 0,15% D15 1116 0,16% 0,39% 0,36% 0,00% 0,00% 0,00% D16 13 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% D17 23 0,00% 0,00% 0,00% 0,00% 0,00% 0,00% D18 223 1,26% 1,43% 6,28% 1,35% 0,90% 0,90% 33 D19 310 1,03% 1,20% 5,81% 0,65% 0,32% 0,65% D20 537 0,15% 0,16% 0,19% 0,00% 0,00% 0,37% Best approach 5 4 8 15 19 16 Table 6. Comparison on the solution quality for the reduced E graphs of the Steiner collection Problem Optimum GA-E MPH P-Tabu F-Tabu VS E1 111 0,00% 0,00% 0,00% 0,00% 0,00% E2 214 0,93% 0,00% 0,00% 0,00% 0,00% E3 4013 0,00% 1,07% 0,42% 0,32% 0,24% E4 5101 0,02% 0,18% 0,00% 0,00% 0,00% E5 8128 0,00% 0,02% 0,00% 0,00% 0,00% E6 73 0,00% 0,00% 0,00% 0,00% 0,00% E7 145 0,00% 2,07% 2,07% 0,00% 0,00% E8 2640 0,23% 1,63% 0,49% 0,42% 1,14% E9 3604 0,19% 1,17% 0,42% 0,14% 0,47% E10 5600 0,00% 0,21% 0,04% 0,04% 0,14% E11 34 0,00% 0,00% 0,00% 0,00% 0,00% E12 67 1,49% 1,49% 1,49% 1,49% 0,00% E13 1280 0,70% 1,88% 0,78% 0,63% 1,33% E14 1732 0,23% 1,04% 0,29% 0,23% 0,64% E15 2784 0,00% 0,22% 0,11% 0,11% 0,00% E16 15 0,00% 0,00% 0,00% 0,00% 0,00% E17 25 0,00% 0,00% 0,00% 0,00% 0,00% E18 564 3,37% 7,62% 2,66% 1,60% 2,66% E19 758 1,26% 4,35% 1,19% 1,19% 1,18% E20 1342 0,00% 0,67% 0,00% 0,00% 0,15% Best approach 12 6 9 14 12 Table 7. Computational times (in seconds) Set of Problems Group Best Case Worst Case Stein-C 1 1 2 Stein-C 2 71 205 Stein-C 3 54 112 Stein-D 1 4 13 Stein-D 2 1,471 4,329 Stein-D 3 166 750 Stein-E 1 6 33 Stein-E 2 4,692 19,892 Stein-E 3 480 3,292