Imene Lahyani

Ismael Bouassida Rodriguez

Mohamed Jmaiel
email: mohamed.jmaiel@enis.rnu.tn

Christophe Chassot
email: chassot@laas.fr

Towards Self Healing Publish/Subscribe System on MANET

Keywords: MANET, publish/subscribe, crisis management applications, Quality of service, reconfiguration actions, Graph grammar

This paper presents a novel approach for Quality of Service aware publish/subscribe systems on Mobile Adhoc Network. Crisis management application seems to be a relevant case study to our approach. On the one hand, our proposed solution monitor the system and analyzes his state towards predicting Quality of Service degradations. Real time prediction is based on the Auto Regressive Integrated Moving Average Formula. On the other hand, our solution defines new reconfiguration actions to be applied once Quality of Service degradation are predicted. In order to evaluate the proposed reparation actions, we define a protocol describing different transformation rules for adaptation using graphs and graph grammar. The Graph Matching and Transformation Engine (GMTE) is used to apply rules to the publish/subscribe system.

I. INTRODUCTION

Mobile Ad Hoc Network (MANET) is a network consisting of a set of mobile hosts communicating with each other without the assistance of base stations [START_REF] Corson | Mobile Ad hoc Networking (MANET): Routing Protocol Performance Issues and Evaluation Considerations[END_REF]. MANET includes mobile nodes that can dynamically self organize into arbitrary ad-hoc network topologies, allowing people and devices to seamlessly work in areas without a preexisting communication infrastructure such as, disaster recovery environments. An efficient way to communicate nodes in such environments is publish/subscribe (pub/sub) paradigm [START_REF] Carzaniga | Design and evaluation of a wide-area event notification service[END_REF] due to the decoupling feature characterizing communicated entities.

A typical case study to our study is crisis management systems (CMS). CMS is composed of a set of mobile devices that communicate via a pub/sub architecture. CMS involve assessing potential threats and finding the best ways to avoid these threats. In such urgent situations, crisis management requires dealing with serious situations before their occurrence. Thus, CMS require high Quality of Service (QoS) and need to consider several criteria's that could affect the quality of nodes and links in the network.

Facing this problem, several studies [START_REF] Cugola | REDS: a reconfigurable dispatching system[END_REF], [START_REF] Denko | A Mobility-Aware and Cross-layer Based Middleware for Mobile Ad Hoc Networks[END_REF] react at the middleware level by introducing some reconfigurations actions in order to repair the system. Although the mentioned systems ensure a better system state by substituting faulty links, they do not give any guarantee to have a fault free system. The main point of failure of these systems remains failure occurrence since its undesirable effects on the system especially on CMS.

This paper provides quality of service (QoS) guarantees in an efficient and scalable manner for CMS where (pub/sub) systems are deployed on MANET. The paper elaborates on an architecture that can support QoS in case of dynamic topology. It presents a decisional model aiming to monitor the system and analyzes his state at runtime by means of a predictive approach using the Auto Regressive Integrated Moving Average model (ARIMA) to forecast future latency values. The same method is used to estimate battery level at runtime towards avoiding nodes crash. Once a QoS degradation is predicted, diagnostic module reasons about the cause of the degradation. This creates opportunities to repair the system. Reconfiguration actions may be introduced to ensure system survivability when QoS degradation is predicted. To achieve the mentioned goal, adaption can be accomplished by migrating a faulty broker to another node in the network ensuring a better QoS. Adaption may also be achieved by splitting an overloaded broker in order to minimize his load. An other alternative is to substitute faulty link by introducing a broker in the link.

In order to evaluate the proposed reconfiguration actions at design time, we propose to built a dynamic protocol of the architecture defining transformation rules. Grammar graph is used to describe the architecture. The Graph Matching and Transformation Engine GMTE is used to apply transformation rules to the graph describing the architecture. This enables to give a general overview describing the network after applying reconfiguration actions.

In the following of this paper, we first survey existing systems dealing with QoS aware pub/sub system. In section II, we present the decisional model predicting failure occurrence. A detailed description of the proposed reconfiguration actions is described in Section III. Thereafter, we present, in Section IV, tools used to evaluate reconfiguration actions on the monitored system. Finally, we conclude and discuss future work in Section V.

II. RELATED WORK

Work dealing with QoS aware pub/sub system when deployed on MANET do not adopt predictive schemes to fulfill their aims. We here present these systems classified in two categories.

In the first one, systems [START_REF] Cugola | REDS: a reconfigurable dispatching system[END_REF], [START_REF] Denko | A Mobility-Aware and Cross-layer Based Middleware for Mobile Ad Hoc Networks[END_REF] react after failure occurrence by proposing a reconfiguration action as a solution to link break or removal. These systems aim to ensure system survivability by maintaining connected peers. The mentioned systems ensure a better system state by substituting faulty links. However, they do not give any guarantee to have a fault free system.

The second category reacts before failure by periodically optimizing system performance taking into account some QoS parameters. For instance, Q [START_REF] Avvenuti | A Cross-Layer Approach for Publish/Subscribe in Mobile Ad Hoc Networks[END_REF] which is an infrastructure for pub/sub systems adapts the overlay network to the changing topology by means of cross-layer interaction. Indeed, to increase the efficiency of communication, the overlay network interacts with the routing layer to obtain information concerning the topology of the network. Armed by these information, it reconfigures itself in a way to keep the shortest paths between brokers. However, the shortest path is not always the optimal path ensuring a minimum of latency. In response to changes in physical network, opportunistic overlays introduced in [START_REF] Yuan | Opportunistic overlays: efficient content delivery in mobile ad hoc networks[END_REF] dynamically adapt event dissemination structures (i.e., broker overlays) with the goal of optimizing end-to-end delays in event delivery. As shortcoming, [START_REF] Yuan | Opportunistic overlays: efficient content delivery in mobile ad hoc networks[END_REF] does not provide a monitoring module to supervise the system, nor a mechanism to specify latency requirement of applications. However, in our work, we continuously monitor latency and we predict failures by computing estimated latency with adaptive threshold. Finally, the work in [START_REF] Minkyong | Efficacy of techniques for responsiveness in a wide-area publish/subscribe system[END_REF] proposes a QoS aware messaging system called Harmony. Separately, for each message topic, Harmony specifies performance requirements. Therefore, this constraint renders this approach applicable only for topic based pub/sub systems. Our approach, however, is not limited to a specific subscription language and is considered as a generic approach.

Besides, contrarily to Harmony, we use adaptive threshold according to the requirements of applications and the variation of the network. Furthermore, Harmony uses a technique consisting in sending topics across multiple paths continuously. This assumes the existence of multiple paths in the network and restricts the approach to only cyclic brokers' topology. However, our approach supports general broker overlay topologies and is not restricted to a specific topology. Although these approaches react before failure, none of them use a predictive strategy. However, all of them optimize the system periodically which can lead to wasting resources by performing unnecessary optimizations operations. Besides, the failure may occur in the system and there is no guarantee to have a free failure system. Thus, up to now, the proposed solutions do not succeed to find an efficient solution to the problem. In this paper, we propose a failure prediction approach based on statistical methods to overcome shortcomings raised by the mentioned work. We ensure system survivability and safety by avoiding failures occurrence and prohibiting QoS degradations.

III. THE PROPOSED APPROACH

The general objective of a crisis management systems (CMS) is to ensure that an abnormal or a catastrophic situation is detected and solved at time and does not go out of hand. Thus, it may require a minimum delay to manage the crisis or even minimize it by handling the situation.

The CMS is composed of a set of mobile devices that communicate via a pub/sub architecture. Thus, the CMS seem to be a relevant case study for our evaluation. CMS require QoS guarantee prohibiting communications interruption or even degradation.

In this paper, we describe our autonomic approach adding valuable self-adaptation to existing pub/sub systems without lowering their performances and characteristics. On the one hand, the proposed approach introduces a decisional model in nodes towards monitoring and analyzing nodes and link states. On the other hand, we describe new reconfiguration actions that can be applied to the system when QoS degradation are predicted. However, how to choose the appropriate reconfiguration action ? How to have a general overview on the network after applying reconfiguration actions?

To reach this goal, we propose to build a dynamic protocol that evaluate repairing actions at design time. A first level of the solution consists in giving network modeling. A second level can be reached by generating rules describing adaptation actions. Finally, the transformation engine GMTE is used to apply rules to the network and to generate the new architecture of the network after repair. This gives the opportunity to evaluate rules regarding performance parameters such network survivability.

A. The decisional module

The decisional module involves monitoring and analysis component. It introduces in each node a component responsible for monitoring and analyzing system state towards providing quality delivery to crisis management applications. The proposed analytical module monitors latency and battery level using the Auto regressive Integrated moving average formula(ARIMA). Finally, decisional module decides about node and link state and reasons about the source of degradation once predicted.

Figure 1 describes the proposed analytical module. 1) Monitoring approach: Monitoring the system aims to observe and store relevant QoS parameter values to assess the current system health. MANETs suffer from limited nodes resources and frequent nodes mobility. Thus, in order to supervise the system, we monitor node state as well as link Monitoring module tries to further monitor for total energy consumption in each node and thus conserves network lifetime. Besides, monitoring link state needs to supervise latency measurements at runtime. ⋄latency:Latency is defined as the time elapsed between the time, T1, when the first broker B1 sends the first bit of a packet to his neighbor B2 and the time, T2, when the second receives the last bit of that packet at 𝑇 1 + Δ𝑇 . Thus, the latency that takes the event between B1 and B2 is defined with the following formula:

𝑙𝑎𝑡𝑒𝑛𝑐𝑦(𝐵 2 , 𝐵 1 ) = 𝑇 2 -𝑇 1 ± (𝑜𝑓𝑓𝑠𝑒𝑡) (1)
where the clock offset refers to difference in the time reported by the two clocks. Monitoring module integrated in each broker intercepts event when translating between neighboring brokers in order to compute latency value. However, if there is no event delivery during a long time, specific messages are used in order to analysis QoS link state. Latency is calculated locally by the broker receiving the event as depicted in Algorithm 1.

⋄Battery level: One important characteristic of MANETs is that nodes are energy constrained. Nodes are battery-operated and frequent recharging or replacement of batteries may be undesirable or even impossible. This makes energy an important QoS parameter, that we should supervise. The above mentioned QoS parameter is located at the OS level and low level of the protocol stack. Thus, we explore a multi layer approach in order to provide information to the middleware layer.

The final step related to the monitoring module is saving QoS parameter values in log files.

2) Analysis approach: To deal with QoS issues, we propose a multi objective QoS model for pub/sub system on MANET.

The proposed approach enables energy efficient by proactively analyzing energy level per node towards reacting before node crush. Moreover, the analysis module analysis latency by estimating latency values using ARIMA model and comparing these values with adjustable thresholds. Adaptation is done using the Auto Regressive Integrated Moving Average Formula.

⋄latency analysis: It's based on comparing estimated latency values with adaptive thresholds.

Analyzing monitored data allow us to detect symptoms that indicate an upcoming failure. Thus, by applying basic mathematical model related to ARIMA model on the monitored latency, we can rewrite :

𝐿𝑎𝑡𝑒𝑛𝑐𝑦 𝑡 = 𝜇+ 𝑝 ∑ 𝑖=1 (𝜑 𝑖 ) * (𝐿𝑎𝑡𝑒𝑛𝑐𝑦 𝑡-𝑖 )+𝜀 𝑡 - 𝑞 ∑ 𝑖=1 (Θ 𝑖 * 𝜀 𝑡-𝑖 ) (2)
Failure prediction is performed by comparing estimated latency with appropriate threshold. In our approach, we take use of adjustable decision threshold. Adaptive thresholds are computed dynamically tacking into account the latest available information using the Exponential Moving Average Formula [START_REF] Lucas | Exponentially weighted moving average control schemes: properties and enhancements[END_REF].

𝑆 𝑡+1 = 𝛼𝐿𝑎𝑡𝑒𝑛𝑐𝑦 𝑡 + (1 -𝛼)𝑆 𝑡 ( 3
)
where 𝛼 is a constant and 𝐿𝑎𝑡𝑒𝑛𝑐𝑦 𝑡 is the measured latency value and 𝑆 𝑡 is the threshold value at present time.

⋄Energy analysis: Node crush may cause system functionality interruption. Thus, energy degradation should be predicted at runtime. In order to forecast energy measures, an algorithm is needed though, that will estimate energy values in the future, then compare these values with thresholds. Therefore, we propose a three step algorithm that first computes estimated energy values using the Auto Regressive Integrated Weighed Moving Average Formula(ARIMA) [?] by applying basic mathematical model related to the ARIMA model on the monitored energy:

𝑒𝑛𝑒𝑟𝑔𝑦 𝑡 = 𝜇 + 𝑝 ∑ 𝑖=1 (𝜑 𝑖 ) * (𝐸𝑛𝑒𝑟𝑔𝑦 𝑡-𝑖 ) + 𝜀 𝑡 - 𝑞 ∑ 𝑖=1 (Θ 𝑖 * 𝜀 𝑡-𝑖 ) (4)
Second, the proposed algorithm compares estimated latency with threshold and then trigger alarms when QoS violation is noticed. Generally, battery level is consider to be critical when it become under 20%.

3) Diagnostic approach: After predicting a QoS degradation, we execute an alternative algorithm that effectively look for the source behind latency degradation.

The key idea of source degradation identification is to study the dependency to any statistical relationship between two random variables or two sets of data. Doing this, we use the correlation method[6] towards studying the relation between latency value variation with parameters affecting latency such us hop count variation or overload variation.

Correlation method requires two steps. The first step consists in computing the Pearson correlation coefficient 𝑃 𝐶𝐶 value studying the dependency between two distributions X1 and X2 as described in the following equation.

𝑃 𝐶𝐶 = 𝑐𝑜𝑣(𝑋1, 𝑋2)/ √ (𝑣𝑎𝑟(𝑋1) -𝑣𝑎𝑟(𝑋2)
The second step consists in executing a significance test of the correlation coefficient by comparing the computed coefficient with a theoretical value issued from the table of the critical value related to 𝑃 𝐶𝐶. Thus, if 𝑃 𝐶𝐶(𝑋1, 𝑋2) is positive and greater than theoretical PCC , this means that X1 and X2 are positively and significatively correlated.

According to the significance of the computed correlation coefficient, our algorithm proceeds in different scenarios resumed in Figure 2.

B. The proposed reconfiguration actions

Once a QoS degradation is predicted, we introduce reconfiguration actions in order to repair the system and ensure its safety and survivability required by CMS. This is not the case for forwards which cannot deliver events to the substituting broker since they do not know its location. In order to solve this problem, we intend to inform forwards about the new location of the defaulting broker. So the new broker looks at the list of forwards and sends to them its address, which corresponds to the new location of the broker. During the migration time, the defaulting broker continues accepting messages from forwards until the end of the migration. This ensures that the system continues to work during migration operation. Once informed about the new location of the defaulting broker, the forwards start sending messages to the substituting broker. Consequently, we ensure that messages are never lost during the migration operation. Finally, the old broker shut down when the new broker starts its operation. Figure 3 summarizes migration operation. Our previous work [START_REF] Lahyani | QoS Aware Publish/Subscribe System Deployed on a Mobile Ad-hoc Network[END_REF] gives more detail about this adaptation action. 2) Splitting nodes: An overloaded broker becomes enable to bypass messages to neighbors. An advantageous solution is to split the broker and especially his subscription table to two others nodes in the network. 3) link substitution: When the link becomes overloaded, substitution is a beneficial solution. This may be done by intercalating a node in the defaulting link. After enumerating the proposed reconfigurations actions, we aim to explore deployment solutions for pub/sub systems on MANET. Doing this, we plan to built a simulator that give first a modeling of the pub/sub system on MANET. Second, the simulator describes rules to each reconfiguration action and finally generates the output graph by applying rules using GMTE1 as a transformation engine.

1) System modeling: Pub/sub system on MANET is modeled by a graph AG containing logical as well as physical nodes and links. Thus, in the proposed graph, there is two types of vertices as follow:

• vertices modeling logical entities, ie, brokers, characterized by : ID broker, IP node, broker state. • vertices modeling physical entities, characterized by :IP node, energy, node state. Also, there is two types of arcs. The first one indicates a physical or a logical link characterized with: type(either physical or logical),hop count between brokers, link state. For example:

• (P, -, Hop, state): indicates a physical link (1 hop),

• (-, L,Hop, state): indicates a logical link (n hops), 2) Rules description: The rule represents the reconfiguration action to apply. Rules are modeled by a rule graph RG containing 4 zones: Inv, Abs, Del, Add, where :

• Inv is the fragment that remains invariant before and after the application of the rule. • Del is the fragment that will be deleted while the rule is applied. • Abs is the fragment that should be absent from the graph in order to apply the rule. • ADD is the fragment that will be added after the application of the rule.

3) The output Graph: After applying rules on the input graph modeling network configuration, the transformation engine GMTE generates the output graph modeling the network with a new configuration.

CMS involves structured groups of robots or military personnel who cooperate for the realization of a common mission. Such application is composed of participants having different roles: a global controller of the mission, coordinators, and several groups of investigators. Each group of investigators is supervised by a coordinator. This hierarchy could be schematized as follows: We propose to elaborate communication between mobile investigators through a content based pub/sub system. In fact, the weak coupling between information producers' and consumers' identities provides P/S system the advantages of asynchronous, multi-point communication, hence freeing the communicating entities from temporal and special requirements for connection.

The communication begins when an investigator either advertises or publishes damage assessment information to find out some rescue in a particular region. Other members can selectively subscribe to the information germane to their task. Typically, the pub/sub system matches the published messages with the dynamically changing interests and routes them to the appropriate recipients.

Every investigator in the network is potentially mobile, it can join, leave the system and have access to the pub/sub service, while also acting as a broker for message forwarding.

Let consider a section of investigator. The graph AG corresponding to the section of investigators (group 3) is depicted in Figure 5.

Our decisional module predicts a QoS degradation in node (𝑁 2 in Figure 5) , Thus, applying the migration rule to the graph in Figure generates the new graph (Figure 6) Future work will evaluate the output graph in order to decide about the appropriate reconfiguration action to be applied in case of QoS degradation. Parameters to evaluate the better output graph may involve liveness of the network. Thus, the better reconfiguration action corresponds to the rule that ensure a longer lifetime to the network. In this paper, we focus on QoS guarantees for crisis management applications composed by a pub/sub system deployed on MANET. Our proposed solution monitor the system and predict failure using the ARIMA model.

In the second part of this paper, new reconfiguration actions are proposed to repair the system. In order to evaluate such reconfiguration actions, we propose a dynamic protocol at design time describing different transformation rules for reconfiguration actions on formal technics of graph grammar. GMTE is used to transform these rules towards enumerating all possible adaptation solutions. Experimentation results will evaluate the proposed adaptation rules.

Figure 1 .

 1 Figure 1. Decisional module based prediction approach

Figure 2 .

 2 Figure 2. Diagnostic approach

Figure 3

 3 Figure 3. Migration operation

Figure 4 .

 4 Figure 4. CMS architecture

Figure 5 .

 5 Figure 5. Network configuration

  Migration operation consists in three different steps: Initialization, Data Transfer and finally Informing neighbors about the new locality. 1.Initialization The first step in broker migration consists in preparing the destination node to host the mobile broker. To this end, a new broker is launched in the selected node. This broker has initially an empty subscription table. 2.Data Transfer It is the most important step in the migration operation. It consists in sending the subscription table to the new broker deployed in the last step. The latter updates its routing table. Consequently, when the new broker receives an event corresponding to one of the subscription stocked in its routing table, it notifies the corresponding subscriber. The defaulting broker sends each subscription received during the migration process to the new one. This prevents messages from being lost. During this step, the initial broker continues to work. 3.Informing Neighbors about the New Locality After the reception of the subscription table, the new broker becomes able to notify the subscribers figuring in its subscription table up on receiving an event.

1) Migration operation

Migration operation means code mobility of brokers as well as data encapsulated within it and which are necessary to its execution.

The Graph Matching and Transformation Engine (GMTE), available at http://homepages.laas.fr/khalil/GMTE


Inken Wohlers

Rumen Andonov

Gunnar W Klau

Optimal DALI protein structure alignment

Keywords: structure alignment, inter-residue distance matrix, exact algorithm, integer linear program, Lagrangian relaxation, DALI matrice de distances entre résidus, algorithme exact, programme en nombres entiers, relaxation Lagrangienne

We present a mathematical model and exact algorithm for protein structure alignment using DALI scoring, which is an NP-hard problem. DALI scoring is based on comparing the inter-residue distance matrices of proteins and is the scoring model of the widely used heuristic DALI program. Our model and algorithm extend an integer linear programming approach previously used for the related contact map overlap problem. To this end, we introduce a novel type of constraint that handles negative structure scores and relax it in a Lagrangian fashion. We also review options that allow to consider less pairs of inter-residue distances explicitly, because their large number makes it difficult to optimize DALI scoring optimally. We use our exact algorithm DALIX to compute many provably score-optimal DALI alignments for the first time, using four data sets of varying structural similarity. Further, using our exact DALIX alignments, it is for the very first time possible to qualitatively benchmark the heuristic DALI program in sound mathematical terms. The results indicate that DALI often computes optimal or close to optimal alignments, but also that in cases of aligning small proteins it tends to fail generating any significant alignment although such an alignment exists.

Introduction

Protein structure alignment, the assignment of structurally equivalent amino acids between protein structures, is an important problem in structural bioinformatics. Detecting and evaluating structural similarities is a standard task that lies at the basis of many applications like, for example, homology detection, fold recognition, protein classification, and functional annotation. For proteins with low sequence similarity, structure alignments improve upon sequence alignments and are successfully applied for detecting functional similarities or distinct evolutionary relationships. Due to the exponential increase in available protein structures, as for example in the protein data bank (PDB) [START_REF] Berman | The Protein Data Bank[END_REF], the problem will receive even more interest in the future.

The structural alignment problem is NP-hard for biologically meaningful scoring schemes used in practice [START_REF] Lathrop | The protein threading problem with sequence amino acid interaction preferences is NP-complete[END_REF]. As a result, almost all algorithms for structure alignment are heuristics; they aim at optimizing a given scoring scheme, but have no notion how much better their computed alignments could be with respect to that scoring scheme. dali (distance matrix alignment) is one of the most widely used structural alignment heuristics [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF]. It is available via the EBI structural analysis tool box, a dedicated server processes about 1500 pairwise alignment user requests a month, and the first dali paper has been cited almost 3000 times, more often than any other structure alignment program 1 .

We present an exact dali algorithm, which we call dalix. It returns either the optimal alignment according to dali scoring, if found within a predefined time limit, or an alignment together with an upper bound on the optimal score. Our algorithm uses techniques from combinatorial optimization. First, we cast the problem into an integer linear program (ILP) whose objective function maximizes the dali score and whose constraints denote that the solution represents a structure alignment. We then relax a few constraints, move them to the objective function and penalize their violation by multiplying them with Lagrangian multipliers. The relaxed problem can then be solved by double dynamic programming, which is a method that is also used in other structure alignment algorithms like ssap [START_REF] Taylor | Protein structure alignment[END_REF] or matras [START_REF] Kawabata | Protein structure comparison using the markov transition model of evolution[END_REF]. Iteratively, multipliers are adjusted and the double dynamic programming is repeated. After a specified number of iterations, the problem is split into subproblems within a branch-and-bound algorithm.

Our mathematical model and algorithm are applicable to any distance matrix-based scoring scheme, e.g. those of ssap [START_REF] Taylor | Protein structure alignment[END_REF], contact map overlap (cmo) [START_REF] Godzik | Flexible algorithm for direct multiple alignment of protein structures and sequences[END_REF], matras [START_REF] Kawabata | Protein structure comparison using the markov transition model of evolution[END_REF] or paul [START_REF] Wohlers | Towards optimal alignment of protein structure distance matrices[END_REF]; for the general framework see [START_REF] Wohlers | Algorithm engineering for optimal alignment of protein structure distance matrices[END_REF]. Here, we focus on dali because it is a popular structural alignment method that performs well in many benchmarks (for example, [START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF]. Provably maximizing the dali scoring is especially difficult, because it uses inter-residue distances between any pair of residues. Our exact algorithm thus needs to explicitly consider O(n 4 ) distance pairs. This has great influence on performance and memory requirements.

Using our algorithm dalix, we are able to compute many dali alignments for the first time to provable optimality and thereby to assess the quality of the dali heuristic. For this purpose we use (i) alignments of SCOPCath domains [START_REF] Csaba | Systematic comparison of SCOP and CATH: a new gold standard for protein structure analysis[END_REF] with lengths between 30 and 50 who share family, superfamily or fold, (ii) alignments from SKOLNICK [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF] of proteins from the same family, and alignments from (iii) the SISY and (iv) RIPC collections [START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF][START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF]. We find that dali is very reliable in returning a good alignment according to dali scoring. Although we find many cases where the dali alignment is not optimal, the difference in score between the heuristic and the optimal alignment is often negligible. When aligning short protein domains, dali's deficiency is that it misses to detect quite a few significant similarities and wrongly reports that no such similarity exists. We also evaluate the weak points of our exact algorithm dalix, which are large proteins and subtle structural similarities. In these cases, the (then suboptimal) alignment returned by dalix is often worse than the heuristic alignment returned by dali.

Mathematical Model and Algorithm

In this section we introduce the alignment graph representation of the structural alignment problem. Based on this representation, we formulate an integer linear program (ILP) that models the problem of finding the alignment of maximum dali score. We devise a Lagrangian relaxation that has been used before for contact map alignment [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] and extend it by relaxing an additional, new type of constraint that is needed for pairs of distances with negative score. We focus in the entire description of the algorithm on these novel constraints, because except for them, model and algorithm are analogous to [START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. In the last subsection we suggest algorithm engineering techniques to improve the performance of our method in practice.

Alignment Graph Representation

We denote the inter-residue distance matrix of a protein A by (A ij ). It is a symmetric square matrix of size n A ×n A , where n A is the length of the protein. A matrix entry A ij is the Euclidean distance between the C α atoms of residues i and j.

We represent the structural alignment problem using an alignment graph. For two proteins of length n A and n B , the alignment graph is a n A × n B product or grid graph as displayed in Figure 1e). Rows represent the residues of protein A and columns the residues of protein B. A node i.k in the alignment graph indicates the alignment of residue i from protein A with residue k from protein B. Directed edges (i.k, j.l) exist between any pair of nodes for which i < j and k < l. Edges are thus south-west to north-east bound, and we refer to nodes i.k and j.l as the tail and head of edge (i.k, j.l), respectively. An edge denotes the matching of distance A ij with distance B kl , see also Figure 1c) andd). An alignment of length n is represented by a set of nodes i 1 .k 1 , i 2 .k 2 , . . . , i n .k n for which i 1 < i 2 < . . . < i n and k 1 < k 2 < . . . < k n . We call such a set an increasing path. A structural alignment comprises additionally all induced edges.

Mathematical Model

We assign binary variables x ik to alignment graph nodes. They indicate whether residue i is aligned with residue k, in which case x ik = 1. An alignment graph edge between nodes i.k and j.l is described by a binary variable y ikjl which denotes whether distance A ij from protein A is aligned with distance B kl in protein B, in which case y ikjl = 1.

We define row ik (j) and col ik (l) as sets of nodes that are either tails of edges with head at i.k or heads of edges with tail at i.k and that mutually contradict because no two of them lie on an increasing path. There are many ways of constructing them, the one we use is introduced in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] and illustrated in Figure 2. The dalix ILP is then given by The black nodes are an illustration of row ik (j) (left) and col ik (l) (right) in the alignment graph. In the displayed situation, j < i and l < k, in which case the colored nodes are sets of mutually exclusive tails of contradicting edges with common head i.k. If j > i and l > k (not displayed), the colored nodes are sets of mutually exclusive heads of contradicting edges with common tail i.k.

max nA-1 i=1 nA j=i+1 nB -1 k=1 nB l=k+1 2s(A ij , B kl )y ikjl + nA i=1 nB k=1 s(A ii , B kk )x ik (1) s.t. x ik ≥ (r,s)∈row ik (j) y ikrs j ∈ [i + 1, n A ], i ∈ [1, n A -1], k ∈ [1, n B -1] (2)
x ik ≥ (r,s)∈col ik (l)

y ikrs l ∈ [k + 1, n B ], i ∈ [1, n A -1], k ∈ [1, n B -1] (3)
x ik ≥ (r,s)∈row ik (j)

y rsik j ∈ [1, i -1], i ∈ [2, n A ], k ∈ [2, n B ] (4)
x ik ≥ (r,s)∈col ik (l)

y rsik l ∈ [1, k -1], i ∈ [2, n A ], k ∈ [2, n B ] (5)
x ik ≤ (r,s)∈row ik (j) s(Ari,B sk )≤0 (y rsik -x rs ) + 1 j ∈ [1, i -1], i ∈ [2, n A ], k ∈ [2, n B ] (6) k l=1 x il + i-1 j=1 x jk ≤ 1 i ∈ [1, n A ], k ∈ [1, n B ] (7)
x, y binary. ( 8)

Here, the scoring function s(•, •) for pairs of inter-residue distances is the dali elastic similarity function [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF],

s(A i,j , B k,l ) = 0.2 - |Ai,j -B k,l | 1 2 (Ai,j +B k,l ) e -( 1 2 (Ai,j +B k,l )/20) 2 i = j and k = l 0.2 otherwise .
An alignment that maximizes the objective function of ILP ( 1)-( 8) is thus an alignment of maximum overall dali score. Based on this overall dali score S(A, B), the dali z-score Z(A, B) Inria is computed according to the formula given in [START_REF] Holm | Dictionary of recurrent domains in protein structures[END_REF],

Z(A, B) = S(A, B) -m(L) 0.5 • m(L) .
The term m(L) for L = √ n A , n B is the approximate mean score and the denominator 0.5 • m(L) estimates the average standard deviation. The z-score thus measures the significance of the detected structural similarity based on an experimentally determined background distribution of dali scores.

All constraints except ( 6) are analogous to the constraints established by [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] for the cmo model. Constraints ( 2) and ( 3) denote that an edge can only be taken if its tail node is activated and if the heads of edges with common tail i.k cannot contradict. Constraints ( 4) and ( 5) denote the reverse situation: an edge can only be taken if its head is activated and the tails of edges with common head i.k do not contradict. Activated nodes have to lie on an increasing path, which is specified by constraints (7). Different from the model for cmo, the dalix model has additional constraints (6). These describe that an edge has to be activated if its head and tail are activated. This is important since according to dali scoring, edge scores can be negative, in which case the remaining constraints allow to omit these edges. Constraints ( 6) are derived from the simple constraints

x ik + x jl -y jlik ≤ 1
for all y jlik with score less than or equal to zero. In these simple constraints, the term x jly jlik can be lifted to (r,s)∈row ik (j) (x rsy rsik ) , since each tail r.s has, according to constraints ( 2) and ( 3), no outgoing edges with contradicting heads i.k.

Lagrangian Approach

Lagrangian Relaxation.

We relax constraints ( 4), ( 5) and ( 6). This means that now an edge can be taken even if its head is not activated and even if its tail contradicts with the tail of another edge (constraints ( 4) and ( 5)) as well as that an edge can be omitted even if its head and tail are activated (constraints ( 6)). Since constraints ( 2) and (3) are not relaxed, still any edge needs to have an activated tail, and since constraints (7) dictate that activated nodes lie on an increasing path, the tails of edges can not contradict in spite of relaxing constraints ( 4) and ( 5). The solution of the relaxed problem is the following: An increasing path of activated nodes, in which each activated node picks outgoing edges of maximum overall score. The heads of these outgoing edges are not necessarily activated. See Figure 3 for a visualization. The relaxation can then straightforwardly be strengthened by constraints

l=s+1,••• ,k y rsil + j=r+1,••• ,i-1 y rsjl ≤ 1 (9) r ∈ [1, n A ], s ∈ [1, n B ], i ∈ [r + 1, n A ], k ∈ [s + 1, n B ] ,
which denote that the heads of outgoing edges picked by each node must (although still not necessarily activated) form an increasing path. The objective function of the relaxed problem is given by Figure 3: Visualization of local profit computation, the solution of the relaxed problem and the feasible solution. Left: Node 1.1 picks its best set of outgoing edges, which are the edges maximizing this node's profit. The corresponding increasing path is colored yellow. Center: The solution of the relaxed problem. It is composed of the increasing path that is the solution of the global problem, colored in blue, together with the outgoing edges that these nodes picked in their respective local problem. The relaxed solution maximizes the sum of profits. For a few edges in the solution of the relaxed problem, the heads are not activated, e.g. for edge (1.1, 3.2). Also, for nodes in the solution, the induced edge is missing, e.g. in the relaxed solution, there is no edge between nodes 3.3 and 6.6. Right: The feasible solution that can be deduced from the relaxed solution. It is composed of the nodes that are activated in the relaxed solution together with all induced edges.

LR(λ) = max nA-1 i=1 nA j=i+1 nB -1 k=1 nB l=k+1 2s(A ij , B kl )y ikjl + nA i=1 nB k=1 s(A ii , B kk )x ik + i,k j∈[1,i-1] λ h ikj   x ik - (r,s)∈row ik (j) y rsik   + i,k l∈[1,k-1] λ v ikl   x ik - (r,s)∈col ik (l) y rsik   + i,k j∈[1,i-1] λ a ikj     1 -x ik + (r,s)∈row ik (j) s(Ari,B sk )≤0 (y rsik -x rs )     . (10)
All λ are greater or equal 0. Here, λ h ikj denotes the multipliers for constraints (4), λ v ikl for constraints (5) and λ a ikj for constraints (6).

Double Dynamic Programming.

The relaxed problem can be solved in time O(n 2 A n 2 B ), where n A and n B are the protein lengths. This follows from the proof in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF] 4) and ( 5) are violated, since the head of edge (1.1, 3.2) is not activated and thus x 32 y 1132 . Lagrangian multipliers will increase the weight of node 3.2 and decrease the weight of edge (1.1, 3.2). Right: Constraint ( 6) is violated, since the edge between the activated nodes is not activated and thus x 11 + x 33y 1133

1. The multipliers for the activation constraints will decrease the weights of the nodes 1.1 and 3.3 and increase the weight of the incident edge. The relaxed problem can be solved by double dynamic programming, i.e., dynamic programming on two levels. For this purpose we solve for each node a local problem and afterwards one global problem. In the local problems, we compute for each node i.k the best set of outgoing edges with their heads on an increasing path by use of dynamic programming. We call the sum of the corresponding edge weights node profit p ik . In the global problem, we assign to each node its profit plus a node score of 0.2 and compute the increasing path of maximum overall weight, again by the use of dynamic programming.

The edge and node scores are adjusted according to the Lagrangian multiplier coefficients of x-and y-variables in the objective function (10) of the relaxed problem. How the Lagrangian multipliers redistribute score between nodes and incident edges in the case of a violated constraint is visualized in Figure 4. In the local problems, we associate to each edge (j.l, i.k) the weight

c jlik (λ) = 2s(A ij , B kl ) -λ h ikj -λ v ikl + λ a ikj if s(A ij , B kl ) ≤ 0 2s(A ij , B kl ) -λ h ikj -λ v ikl otherwise .
Once profits p ik have been computed by solving the local problem for each node i.k, we solve the global problem. In the global problem, the node weights are given by

c ik = p ik + j∈[1,i-1] λ h ikj + l∈[1,k-1] λ v ikl - j∈[1,i-1]     λ a ikj - j:r.s∈row ik (j) ∧s(Ari,B sk )≤0 λ a ikj     .
The solution of (10) subject to ( 2), ( 3) and ( 7) has objective function value LR(λ) which is an upper bound for the original problem. It comprises the set of nodes that solve the global problem, x, together with the set of edges composing the solutions of the local problems, ȳ, for those tail nodes that are in the solution of the global problem. Nodes corresponding to x together with their induced edges represent a feasible solution for the original problem, i.e., a structural alignment, with objective function value Z lb , which constitutes a lower bound. This is visualized in Figure 3.

Updating Lagrangian Multipliers.

Before each iteration t + 1, the Lagrangian multipliers are adjusted by subgradient descent,

λ t+1 = max{0, λ t -θ t g t },
with step size

θ t = α[LR(λ t ) -Z lb ] [(g h ) t ] 2 + [(g v ) t ] 2 + [(g a ) t ] 2 .
For updating the multipliers, the gradients are computed as follows

g h ikj = xik - (r,s)∈row ik (j) ȳrsik ∈ {-1, 0, 1} (11)
g v ikl = xik - (r,s)∈col ik (l) ȳrsik ∈ {-1, 0, 1} (12)
g a ikj = 1 -xik + (r,s)∈row ik (j) s(Ari,B sk )≤0 (ȳ rsik -xrs ) ∈ {-1, 0, 1} . (13)
There are only three situations in which gradient g a ikj for the activation constraints is nonzero. The first situation is g a ikj = 1, if xik = 0, (r,s)∈row ik (j) s(Ari,B sk )≤0 ȳrsik = 0;

(r,s)∈row ik (j) s(Ari,B sk )≤0 xrs = 0.

In this case λ a ikj is decreased. Since λ a ikj ≥ 0, we have to check for this situation only for indices i, k and with multiplier λ a ikj in the current iteration. The gradient g a ikj is also equal to 1, if

xik = 0, (r,s)∈row ik (j) s(Ari,B sk )≤0 ȳrsik = 1, (r,s)∈row ik (j) s(Ari,B sk )≤0 xrs = 1.
This second situation is analogous to the situation detected by identifying violated constraints (4): The tail node of an edge is activated, but its head is not. Multiplier λ a ikj is decreased, as in the first situation, which is only possible if the current multiplier λ a ikj is greater than zero.

Furthermore, in the third situation

g a ikj = -1, if xik = 1, (r,s)∈row ik (j) s(Ari,B sk )≤0 ȳrsik = 0, (r,s)∈row ik (j) s(Ari,B sk )≤0 xrs = 1.
Here, we have to identify pairs of nodes that are both activated and are connected by an edge which is not activated. Since we have to check only pairs of nodes in the solution, this can be done quickly in time O(n 2 ). After a specified number of Lagrangian iterations, the problem is split into four subproblems. The entire branch-and-bound framework is described in [START_REF] Andonov | Maximum contact map overlap revisited[END_REF].

Algorithm Engineering

In most cases preprocessing is beneficial. The idea is to exclude nodes from the alignment graph if they provably cannot be part an alignment. In order to do so, a good feasible solution is needed; it can, for example, be provided by the heuristic dali algorithm. Then, for each alignment graph node, an overestimation of the the best alignment including this node is computed. If this upper bound is less than the score of a known feasible solution, the node cannot be part of the optimal solution and will be discarded. The benefit of preprocessing is twofold: first, the computation time decreases if only a subset of nodes in the alignment graph needs to be considered. Second, less memory is needed, because together with a node also all its incoming and outgoing edges can be discarded. Since each alignment graph node can be handled one after another, there is no memory issue during preprocessing. As a result, successful filtering thus allows to compute alignments of protein pairs that otherwise would not fit into memory.

Another attempt to handle the memory requirements of especially large proteins is the exclusion of inter-residue distances between residues that belong to different domains. The orientation of domains to each other is thus entirely disregarded during alignment. As a result, alignments are scored on the base of different criteria, and thus the optimal alignments in both cases may differ. We evaluate empirically the influence of omitting inter-domain distances on the alignment accuracy.

Data Sets and Experimental Setup

The SKOLNICK data set consists of 40 proteins with length between 97 and 255 residues belonging to 5 protein families. This easy dataset has been used extensively for clustering of protein structures, see for example [START_REF] Caprara | 1001 optimal PDB structure alignments: integer programming methods for finding the maximum contact map overlap[END_REF][START_REF] Xie | A reduction-based exact algorithm for the contact map overlap problem[END_REF][START_REF] Malod-Dognin | Maximum cliques in protein structure comparison[END_REF], 2011;[START_REF] Andonov | Maximum contact map overlap revisited[END_REF]. We align only protein pairs from the same family, which amounts to 164 SKOLNICK instances.

SCOPCath [START_REF] Csaba | Systematic comparison of SCOP and CATH: a new gold standard for protein structure analysis[END_REF] is a benchmark containing 6759 domains that are consistently classified in SCOP (Murzin et al., 1995, version 1.75) andCATH (Greene et al., 2007, version 3.2.0) and that have a pairwise sequence similarity of less than 50%. We align all SCOPCath domains with 30 up to 50 residues which belong to the same family (386 pairs), to different families but to the same superfamily (151 pairs), and to different superfamilies but the same fold (926 pairs). We limited the length to maximally 50 residues to obtain alignments for which our algorithm can explore multiple branch-and-bound nodes within a few minutes.

SISY and RIPC [START_REF] Mayr | Comparative analysis of protein structure alignments[END_REF][START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF] are datasets of manually curated structural alignments assembled from the Sisyphus collection [START_REF] Andreeva | SISYPHUS-structural alignments for proteins with non-trivial relationships[END_REF], which are difficult for alignment programs because of repetitions, large indels, circular permutations, conformational variability, etc. The consolidated SISY and RIPC sets consist of 98 and 22 alignments, respectively. With consolidated we denote the subsets that have been consulted for evaluation in [START_REF] Berbalk | Accuracy analysis of multiple structure alignments[END_REF].

For comparison with dali we use DaliLite version 3.3. dali computes a number of alignments and ranks them according to z-score of protein unfolding units. Therefore, we parse all alignments returned by dali and consider the one with largest dali score. Note that maximizing the dali score will also maximize the z-score of the entire alignment. Nonetheless, in an attempt to report and rank interesting high local similarities, dali also computes z-scores for parts of the alignment, the protein unfolding units, in which case a suboptimal alignment can receive the highest z-score.

We compute pairwise alignments on cluster nodes each equipped with two quad core 2.26 GHz Intel Xeon processors and 24 GB main running 64 bit Linux. In each branchand-bound node we compute 1000 Lagrange iterations. For SISY, RIPC and SKOLNICK, a maximum running time of 30 CPU hours per instance is applied and for the short SCOPCath instances a time limit of 30 CPU minutes per instance.

Results and Discussion

We assess the capability of our algorithm to compute optimal alignments with respect to the dali scoring function on (i) 164 SKOLNICK alignments, (ii) 1463 SCOPCath alignments, (iii) 98 SISY alignments and (iv) 22 RIPC alignments. On the two later datasets, we also evaluate the alignment accuracy with respect to the manually curated reference alignment, which is defined as the percentage of correctly aligned residue pairs. We compare our alignments and their scores to those determined by the dali program. The results are summarized in Table 1 andFigure in each dataset in the first row and the number of alignments for which either dalix or dali detected an alignment with z-score greater 0 in the second row. Given that only alignments with z-score greater zero are significant and that dali reports only significant alignments, we only consider those alignments in the following table rows. "dalix optimal" denotes the number of dalix alignments that have been computed to optimality and "dali optimal" the dali alignments thereof that are also provably optimal. "dalix better" lists the number of dalix alignments with higher dali score than the dali alignment. "dali better" denotes the number of dali alignments that are better than the (not yet optimal) dalix alignment. The dalix time limit for SCOPCath alignments is 30 CPU minutes and for all other data sets 30 CPU hours.

SKOLNICK.

136 of the 164 alignments were computed to optimality (83%). For 123 alignments (75%), the heuristic dali solution was improved, but never more than 3%, and for 38 (23%) the heuristic solution was proven to be optimal. Only three dali alignments were slightly better than the corresponding dalix alignment. Results for SKOLNICK indicate that dali computes optimal or close to optimal alignments in the case of distinct structural similarities on family level. The results also demonstrate that it is feasible to compute structural alignments to optimality in the case of considerable structural similarity. Figure 5: The barplot bins the percentages of dali score improvement for the cases in which the dalix alignment has positive z-score and is better than the dali alignment. On family level, these are 278, on superfamily level 118 and on fold level 258 alignments. The improvement is computed with respect to the dalix alignment. The dalix computation time limit is 30 CPU minutes. For most alignments, the score improvement is small. The large percentage of alignments with improvement between 90 and 100% traces back to the large number of protein pairs for which dali falsely reports that there is no structural similarity.

SCOPCath.

When aligning the short SCOPCath domains, for 661 (45%) neither dali nor dalix could compute an alignment with positive z-score, especially on fold level. It is likely, but unfortunately not proven by our upper bounds, that no such alignment exists in many cases. This situation illustrates that it is difficult to design a scoring scheme and algorithm that reliably detects structural similarities on different classification levels and discriminates them from spurious similarities. An exact algorithm maximizing a "perfect" scoring function is expected to return a significant alignment for all protein pairs investigated in this paper. Given the dali score and the dali and dalix algorithms, this is not the case. We thus exclude protein pairs from the analysis for which no algorithm returns an alignment with positive z-score.

From the 359 short SCOPCath alignments of domains from the same family, 143 (40%) were solved to optimality within the time limit of 1800 seconds, and most within a few seconds. There are differences between optimal alignments and heuristic ones computed by the dali program, which are quantified in Figure 5. First, dali fails to detect 83 significant alignments with zscores up to 5 and falsely reports that there are no structural similarities. dalix computes 143 alignments (40%) to optimality. From these cases, in which the provable top-scoring alignment is known, dali returns 50 optimal and 93 close to optimal alignments. Altogether, our exact algorithm improves the heuristic solution in 287 cases (80%) by median 7%. In 16 cases, the solution returned by the exact algorithm after 1800 seconds is worse than the heuristic solution.

Computing exact alignments becomes more difficult when structural similarity gets less pronounced. For the 141 SCOPCath alignments of proteins that the same superfamily, but not the same family, only 14 (10%) are computed to optimality within 1800 seconds. Nonetheless, the exact algorithm improves the heuristic alignments in 118 cases (84%), as visualized in Figure 5. In 24 of them, dalix returns a significant alignment that is entirely missed by dali. In 14 cases, the dali heuristic alignment has larger dali score than the dalix alignment.

On fold level, only 31 of the 302 alignments were computed to optimality. 258 alignments (85%) returned by the exact algorithm are better than those returned by the heuristic algorithm, by median of 27%, see Fig. 5. Also here, the dali score for alignments produced by dali or dalix is usually very similar, but the exact algorithm detects significant alignments with z-scores up to 2.5 that are missed by the heuristic. In 123 cases (41%), dalix determines a significant alignment that is missed by dali. In 30 cases the heuristic alignment is better than the one returned by our algorithm.

Figure 4.0.2 visualizes the alignment traces of three alignments, from family, superfamily or fold level, respectively. These are the instances for which the dalix alignment improves the dali alignment with the largest relative score difference from all alignments of the respective similarity level.

SISY and RIPC.

If structural similarities are less pronounced or locally confined, determining the optimal alignment becomes inordinately more difficult. Furthermore, protein length is a problem for our algorithm, since the number of distance pairs grows quadratic with protein length. From the SISY set, whose difficulty is also confirmed by low alignment z-scores, only 62 alignments fit into memory. From these, 11 were solved to optimality in, on average, about 6000 seconds. Three of these optimal alignments were also detected by the dali heuristic, and the remaining 8 non-optimal heuristic alignments were less than 3% worse than the optimal ones. Altogether, our exact algorithm improves the heuristic alignments in 31 cases by a median of 2%. In 27 instances, the heuristic alignment performs by a median of 19% better than the the alignment returned by our algorithm. The reason is that the instances in the SISY set are large and of subtle similarity, which is both disadvantageous for our algorithm.

From the 22 RIPC alignments, 11 fit into memory. Of those, two were computed to optimality. An improvement of median 4% over the dali solution was found in 6 cases. The heuristic solution is better than the solution of our algorithm for 5 alignments, with a large median dali score difference of 100%. In 3 cases dalix thus does not detect an existing alignment with positive z-score.

On SISY and RIPC we observe that, within the given time limit, our exact algorithm fails to produce good alignments for large or remotely similar protein pairs. Furthermore, we find that improving the dali score does not necessarily implicate a higher alignment accuracy with respect to manually curated reference alignments: 13 dalix alignments from SISY have larger dali score than the dali alignment, but slightly less alignment accuracy with respect to the reference alignment. Only for three alignments an increased dali score also results in a slightly increased alignment accuracy. For the RIPC data set, in two cases an improved dali score increased the alignment accuracy, but in two other cases alignment accuracy decreased. Inspection of a few of these instances using CSA [START_REF] Wohlers | CSA: Comprehensive comparison of pairwise protein structure alignments[END_REF], our web server for comparative structural alignment, shows that in dalix alignments with larger dali score but less alignment accuracy often short, additional gaps have been inserted rendering the alignment more scattered and less intuitive than the respective dali alignment.

We evaluated the benefit of preprocessing as described in Section 2.4. Four more SISY and one more RIPC alignments fit into memory. Two of these SISY alignments wer computed to optimality. In most cases the alignment returned when additional preprocessing is used is better with respect to dali score than the one obtained without preprocessing. In those instances that were solved, the overall running time including preprocessing is in most cases significantly smaller than without preprocessing, although the preprocessing itself takes, depending on protein length, up to 13 minutes. Similar observations hold for the RIPC alignments. Here, one previously unsolved protein pair is aligned optimally when preprocessing is used, in as little as 91 seconds.

Large proteins that share only little structural similarity cause memory and performance problems for our algorithm. It was thus only possible to fit 62 SISY (63%) and 11 RIPC (50%) alignments into memory. In order to reduce memory consumption while keeping alignments very close to the ones computed by the dali program, we made another SISY and RIPC evaluation in which we excluded all inter-domain distances. Now, 92 SISY (94%) and 19 RIPC (86%) alignments fit into memory. While it allows to compute alignments for larger proteins and in shorter computing time, omitting inter-domain distances has little influence on alignment accuracy with respect to the manually curated SISY and RIPC reference alignments.

Conclusion

We presented the first exact general algorithm for distance matrix alignment and implemented and evaluated it for the dali scoring function. We used our implementation dalix to benchmark the popular dali structural alignment method. We found that computing a dali alignment to optimality is feasible if the proteins are not too large or if there is a clear structural similarity. In these cases we noticed that dali alignments are often not optimal, but nonetheless almost always very close to optimal. Further, for short protein pairs we detected that although dali rarely returns a poor alignment, it tends to entirely miss structural similarities and wrongly does not return any alignment. These findings confirm the high quality of heuristic dali alignments and identify possible improvements.

Inria

RESEARCH CENTRE RENNES -BRETAGNE ATLANTIQUE

Inria

Figure 1 :

 1 Figure 1: Different protein and alignment representations of protein A with n A = 4 residues and protein B with n B = 3 residues. a) The amino acid sequence representation. Instead of the amino acid, the corresponding residue number is given. The alignment on the right denotes which residues structurally match. The second residue of protein A is unaligned. b) The corresponding superposition. Given the alignment of residue 1 with 1, 3 with 2 and 4 with 3, protein B is translated and rotated such that the superposition minimizes RMSD. c) The residues are arranged on a horizontal line and the inter-residue distances displayed. The alignment on the right highlights three pairs of aligned inter-residue distances. d) The inter-residue distance matrices. The super-imposed and collapsed distance matrices denote the alignment and highlight

Figure 2 :

 2 Figure2: The black nodes are an illustration of row ik (j) (left) and col ik (l) (right) in the alignment graph. In the displayed situation, j < i and l < k, in which case the colored nodes are sets of mutually exclusive tails of contradicting edges with common head i.k. If j > i and l > k (not displayed), the colored nodes are sets of mutually exclusive heads of contradicting edges with common tail i.k.

Figure 4 :

 4 Figure 4: An example of the redistribution of score between nodes and edges in the case of violated constraints. Left: Constraints (4) and (5) are violated, since the head of edge (1.1, 3.2) is not activated and thus x 32 y 1132 . Lagrangian multipliers will increase the weight of node 3.2 and decrease the weight of edge (1.1, 3.2). Right: Constraint (6) is violated, since the edge between the activated nodes is not activated and thus x 11 + x 33y 1133 1. The multipliers for the activation constraints will decrease the weights of the nodes 1.1 and 3.3 and increase the weight of the incident edge.

2 nB 2 ,

 22 graph edges. In the case of the dali ILP, |V | = n A n B and |E| = nA from which the stated complexity follows.

Inria

Inria

Figure 6 :

 6 Figure 6: Three SCOPCath examples in which the dalix alignment improves the dali alignment significantly. Proteins from common family, superfamily and fold are aligned (SCOPCath IDs d2glia1 vs. d2glia4, d1b9wa2 vs. d1dx5j3, and d1pfwa3 vs. d1ryta2). Each alignment has the largest percentage of score difference of alignments from the respective level of similarity (45%, 35% and 47%). Residues aligned by dalix are colored blue, residues aligned by dali red and residues aligned by both methods gray.

  , according to which the complexity is O(|V | + |E|), with |V | the number of alignment graph nodes and |E| the number of alignment

	Inria

Table 1 :

 1 5. Comparison of dalix and dali alignments. The table lists the number alignments

		SKOLNICK		SCOPCath		SISY RIPC
			Family	Superfamily Fold
	Alignments	164	386	151	926	62	11
	Positive z-score	164	359	141	302	61	11
	dalix optimal	136	143	14	31	11	2
	dali optimal	38	50	5	5	3	0
	dalix better	123	287	118	258	31	6
	dali better	3	16	14	30	27	5

  ). Each alignment has the largest percentage of score difference of alignments from the respective level of similarity (45%, 35% and 47%). Residues aligned by dalix are colored blue, residues aligned by dali red and residues aligned by both methods gray.

RR n°7915

Acknowledgements

This work was supported partly by DFG grant KL 1390/2-1 and by an INRIA internship grant. R. Andonov is supported by BioWIC ANR-08-SEGI-005 project. Computational experiments were sponsored by the NCF for the use of supercomputer facilities with financial support from NWO.