L Thevenard

H T Zeng

D Petit

R P Cowburn

Macrospin limit and configurational anisotropy in nanoscale Permalloy triangles

Keywords: nanomagnetism, magnetic anisotropy, longitudinal Kerr effect, configurational anisotropy

In Permalloy submicron triangles, configurational anisotropy -a higher-order form of shape anisotropy -yields three equivalent easy axes, imposed by the structures' symmetry order. Supported by micromagnetic simulations, an experimental method was devised to evaluate the nanostructure dimensions for which a Stoner-Wohlfarth type of reversal could be used to describe this particular magnetic anisotropy. In this regime, a straightforward procedure using an in-plane rotating field allowed us to quantify experimentally the six-fold anisotropy fields for triangles of different thicknesses and sizes.

Introduction

Much has been accomplished in previous years to shrink the feature size of nanomagnets in view of increasing the storage density of Magnetic Random Access Memories (MRAM) and hard disk drives, but also to develop materials exhibiting a magnetic anisotropy strong enough to resist thermal fluctuations at small dimensions, while remaining switchable by accessible magnetic fields [1,2]. In the ubiquitous soft alloy Ni 81 Fe 19 (Permalloy, Py), the magnetic anisotropy can, in certain cases, be tailored by the geometry of the nanostructure. For instance, appropriately sized Py ellipses have two stable states, where the magnetization lies along the longer in-plane dimension of the element.

Another type of anisotropy exists in this material that remains to be harnessed to technological applications: configurational anisotropy (CA) [3]. This phenomenon is a direct effect of the rotational symmetry order (n) of a nanostruture on its magnetic anisotropy: triangles Email addresses: l.thevenard@imperial.ac.uk (L. Thevenard) (n=3) evidence a six-fold anisotropy, squares (n=4) a four-fold anisotropy, and pentagons (n=5) a ten-fold anisotropy [4,5,6]. Note that the necessity for easy and hard axes to present the same symmetry leads to frequency doubling for odd orders of n [7]. Configurational anisotropy relies on the fact that at small dimensions, a uniform magnetization cannot be sustained anymore in non-ellipsoidal structures, leading to a sizable deformation of the spin arrangement into an energetically more favorable state, balancing exchange and demagnetization energy costs. This effect could for instance allow the coding of multiple "bits" per nanostructure, or be used for complex toggling or switching mechanisms in MRAM type structures.

To be of interest however, it is necessary to have an excellent control of this anisotropy. Different approaches have been proposed to measure higher-order magnetic anisotropies, such as Rotating field Magneto-Optic Kerr Effect (ROTMOKE) [8] in magnetic thin films, or Modulated Field Magneto-Optical Anisometry (MFMA) [9,6] in Supermalloy (Ni 81 Fe 14 Mo 5 ) and Fe nanostructures presenting CA. In both methods, large static fields are used to impose macrospin-like behav-ior in the structure, in order to interpret the data within a Stoner-Wohlfarth (SW) model and extract a value for the strength of this anisotropy. Here we present an alternative reliable and straight-forward experimental technique to obtain the six-fold anisotropy field of submicron triangles, and we use micromagnetic simulations to define a criterion estimating the maximum triangle dimensions up to which this parameter can indeed be defined.

Samples

Equilateral Py triangles were fabricated on Silicon by a 20kV electron beam lithography, followed by thermal evaporation and a lift-off process. A typical scanning electron micrograph of these triangles was presented in Ref. [5]. Their widths and thicknesses varied between e=150-300 nm, and t=6-26 nm, and they were spaced by e into large 25 × 25µm 2 arrays. In these nanostructures, configurational anisotropy is responsible for two possible magnetization states: a 'buckle' state where the global magnetization lies parallel to the edge, and the spins buckle from one corner to the other; or the 'Ystate', where M follows a triangle bisector, the spins splaying in (or out) from two corners toward the third [9]. We have shown recently that for triangle dimensions above 100 nm wide/10 nm thick, only the buckle state can be stabilized [5], such that the 0 • [modulus 60 • ] directions are easy axes, and the 30 ] directions hard axes.

Stoner-Wohlfarth switching astroid

While nanostructures exhibiting CA cannot, by definition, exhibit a macrospin behavior in the strict sense, it is tempting to compare their reversal to a Stoner-Wohlfarth (SW) model [10,4], where the global magnetization is defined as M = (M s , Θ), and the pseudo macrospin anisotropy is described by the simple energy functional E 6 = K 6 sin 2 (3Θ), with Θ defined with respect to the base of the triangle and H 6 = 2K 6 /M s . Here K 6 is positive, giving the buckles to be easy axes as is the case in our samples. To compute a SW switching astroid, a Zeeman term was then added to this expression and the total energy E studied as a function of the applied field H = (H, φ), with φ once again defined with respect to the base of the triangle. The astroid is constructed as follows, taking the example where the field H is ramped along the hard axis φ = 30 • (Fig. 1c). Starting from the H = 0 energy functional E(Θ), a positive field is ramped up to H 1 for which a first transition is observed (120 • to 60 • ) upon a local sign inversion of d 2 E/dΘ 2 . The amplitude of this transition is 60 • . With increasing field, the 180 • to 60 • transition occurs at H 2 (amplitude of 120 • ) and finally 60 • to 0 • occurs at H 3 (once again a 60 • transition). For a field applied exactly along 30 • , the 0 • and 60 • positions are strictly speaking degenerate and the transition field then diverges; this actually leads to a computational artifact creating a slight asymmetry along the 30 • [modulus 60 • ] directions. The 60 • transition shown in Fig. 1c (field H 3 ) will however be possible as soon as the field is brought away from the bisector. Applying the field in the opposite direction (30 • + 180 • ) gives transitions of equal amplitude at identical fields. For the field applied within 15 • of the easy axes 0 • [60 • ], an initial 60 • transition is followed at much higher field by a full 180 • transition. Note that there are two different types of 60 • transitions: the "lower-field" ones with maxima along 0 • [60 • ], and the "higher-field ones", with maxima along 30 • [60 • ]. For each field direction φ, the transition fields H i are evaluated in this way, and their locus plotted on the complete astroid, where transitions of different amplitudes 60, 120 or 180 • have been labeled (Fig. 1a).

Very much like the SW astroid of an ellipse presents maxima along its long and short axes directions, Fig. 1a presents equal maxima along the easy and hard axes. The envelope of the full astroid is therefore a twelvepointed star. The sole effect of increasing the anisotropy field is to expand the astroid, keeping all features identical, as shown in Fig. 1b where three SW astroids were calculated for values of K 6 between 2000 and 12000 J/m 3 (for clarity, only the envelope of the astroid is shown).

The rotating field method

Micromagnetic simulations

Following these calculations, a possibility to extract a value of K 6 from experimental data would be to use an experimental switching astroid. The latter has been measured [5], but this method tends to yield fairly large error bars on the transition fields, as it requires to isolate specific transitions as small kinks in the hysteresis loops. Another solution consists of forcing 60 • rotations of the magnetization in the triangles, for instance by using a large rotating field. By measuring the field amplitude required to observe a magnetization flip from one buckle to the next, the 60 • transition fields can be obtained and compared to the equivalent feature in the SW model.

In order to test the validity of this approach, numerical simulations (OOMMF package [11]) were first 2 done on 310 nm wide, 10 nm thick triangles, using a 5 nm cell size, M s = 800 kA/m, A = 13 pJ/m and α = 0.5. A rounding of ⊘ = 60 nm was moreover included to take into account the physical rounding of our structures [5], and the triangular mesh was titled by 15 • to spread pixelation effects. A counterclockwise field of amplitude H rot was then applied, starting from the positive x direction with the magnetization initialized in a 0 • buckle. Unsurprisingly, for low amplitude rotating fields (H rot ≤ 150 Oe), no rotation of the global magnetization M was observed at all, whereas for H rot ≥ 1000 Oe, M followed exactly the field. In the intermediate regime, the temporal evolution of the magnetization and the corresponding spin configuration are as shown in Fig. 2a,b (for clarity, the triangles were not shown tilted by 15 • ). Six abrupt 60 • rotations of M are observed during a full field cycle, with the magnetization flipping from one buckle to the next. The amplitude of the transitions in Fig. 2a corresponds to the projec- tion of the global magnetization drawn schematically by large white arrows in Fig. 2b. The simulations therefore seem to show that it is possible to address a particular transition, the 60 • rotation, by using a rotating field of carefully chosen amplitude.

Experimental methods

The effect of a 1 Hz rotating field was then studied experimentally using Magneto-Optical Kerr Effect (MOKE) in the longitudinal geometry (x axis in Fig. 2b). The laser was focused into a 5 µm spot on each array, and the field created by a combination of inplane fields (H x ,H y ) from a quadrupole electromagnet. A typical hysteresis loop is presented in Fig. 2c, where the longitudinal magnetization M x is plotted against the x component of a counterclockwise rotating field H rot = 324 Oe. The loop exhibits 6 very clear and sudden changes in M x , in two sets of 3 symmetrical transitions (longitudinal field sweeping down or up). Comparing with the simulations, these can be identified as successive 60 • rotations of the magnetization occurring at specific projections along x of the rotating field: h x1 = ±270 Oe, h x2 = ±17 Oe and h x3 = ±291 Oe. Because the amplitude of the field remains constant during a full cycle, it is immediate to calculate the angle of the field at the moment of the transition, using φ i = cos -1 (h xi /H rot ), or φ i = cos -1 (h yi /h xi ). This yields for the loop shown:

φ 1 = 34 • [180 • ], φ 2 = 93 • [180 • ] and φ 3 = 154 • [180 • ].
For a triangle in a 0 • buckle, a rotating field of 324 Oe will therefore induce a rotation toward the 60 • buckle just as it crosses the φ 1 = 34 • direction. Faraday effects and other perturbations tend to introduce an arbitrary slope to the loops, which we have subtracted for clarity; the amplitude of the transitions are therefore not directly comparable with the simulations .

Results

For each array, M(H) loops were measured under increasing values of H rot , starting from the lowest one giving a 6-stepped loop. The data was then represented as a polar plot representing the rotating field amplitude as a function of the transition angles φ i (Fig. 3). Note that both +60 • and -60 • transitions are represented. They were obtained by applying the field counterclockwise or clockwise. The experimental H rot (φ) data evidences the same characteristics as the computed higher-field 60 • transitions either side of φ = 30 • [60 • ] in the SW astroid. Mainly, if the magnetization lies along the bottom edge of the triangle (0 • buckle), it can only flip by ±60 • to the following buckle under the torque of a field applied along ±30 • at the least. When this angle is increased, the field required to obtain a transition strongly decreases, and ceases to be measurable at about φ = ±45 • . This is naturally reminiscent of the magnetization reversal mechanism in ellipses, where the additional torque from an H y field perpendicular to H x greatly reduces the switching field.

We then attempted to fit the experimental H rot (φ) data to the higher-field 60 • transitions calculated in the SW astroid (Fig. 1), the only adjustable parameter being K 6 . For the plot shown in Fig. 3 for instance (300 nm wide, 9 nm thick triangles), the best fit is thus obtained for K 6 = 10560 ± 300 J/m 3 , or expressed in term of anisotropy field, H 6 = 264 ± 8 Oe (M s = 800 kA/m).

Establishing a macrospin limit

This procedure was applied to triangles of different widths (e =150-300 nm), and thicknesses (t =6-26 nm). For all of these arrays, 6-steps hysteresis loops were easily obtained under a rotating field, and a H rot (φ) polar plot of 60 • transitions readily extracted. Before retrieving H 6 values from these data however, the validity limits of this method needs to be questioned. Indeed, while the 6-stepped hysteresis loops seem to warrant a macrospin-like behavior, and therefore legitimize the approach described above, it can be argued that the existence of a 60 • transition is not a good criteria for a macrospin behavior in our nanostrutures, being an easy, small amplitude 'sliding' of spins from one buckle to the next. A full 180 • reversal, or any procedure allowing the magnetization to relax would be a more reliable control test. The experimental switching astroid now becomes a useful tool, as it is measured by applying a saturating field along a direction -φ, letting M relax while reversing the field direction, and observing the configurational changes of the spins leading to a full reversal of the magnetization along +φ.

A micromagnetic simulation of the switching astroid was then done following the field sequence described above, and the fields inducing abrupt changes of the magnetization plotted versus the direction of the applied field. The dimensions (thickness, width) investigated were identical to the experimental ones. For thicknesses t ≤ 10 nm at all triangle widths, as well as the t=15 nm, e=150 nm structure, a unique shape was obtained for the astroid, an example of which is shown in Fig. 4a (first quadrant of the simulation for the t=10 nm, e=300 nm triangle). Having been described in detail in Ref. [5], only the main differences and similarities with a macrospin type of reversal (Fig. 1a) will be highlighted. For fields reversed along directions close to the easy axes 0 • [60 • ], the magnetization does a complete 180 • reversal, going directly from state (i) to (iii) in Fig. 4c. For fields close to the hard axes 30 • [60 • ], it reverses in two steps, a 120 • transition followed by a 60 • transition (closed symbols in Fig. 4a). The latter is the exact same transition observed under a rotating field, or calculated using the SW model. Contrary to the SW astroid, the peaks along 0 • [60 • ] are lower than the ones at 30 • [60 • ]. This is a direct consequence of the complex and different micromagnetic configurations in the easy and hard axes states. In both easy and hard-axis types of reversal, note that there is one less transition compared to the SW model: for all field directions, the first 60 • transition does not occur in the micromagnetic simulations. Moreover, the locus of the 120 • transition fields has a quite different shape in both cases. On the one hand, the macrospin model identifies local energy minima, without favoring any initial or final configuration: the first 120 • to 60 • transition H 1 in Fig. 1c can only occur if the initial configuration of M is along 120 • for instance. In the micromagnetic simulations on the other hand, the triangle is first fully saturated in one direction, thereby fixing the initial magnetic configuration, before an opposite field is applied.

When the triangle width e is increased above 150 nm (e ≥200 nm) for t ≥ 15 nm, the envelope of the astroid remains the same, a 12-pointed star, but a new transition appears at low fields, as shown in Fig. 4b by the dashed line. Indeed, the simulations evidence for any applied field direction a passage by an intermediary vortex state (ii) shown in Fig. 4c. This vortex is eventually expelled as a 60 • transition does occur at higher fields. Obtaining a 6-stepped hysteresis loop under rotating field (series of 60 • transitions) is therefore not a good 'macrospin' criterion. Following the simulations however, a simple remanence loop along the 0 • direction of the triangle should evidence the passage by a vortex state, seen as the presence or absence of a second step in the loop, and determine whether it is legitimate to extract an H 6 value from the data obtained under a rotating field.

Remanence loops along the 0 • direction were measured for the different arrays. A clear cut-off dimension appears below which the reversal is a single step 180 • reversal (thickness 15 nm, width 300 nm), as shown in Fig. 5, open symbols. Above this value (thickness 18 nm, width 150 nm, full line in Fig. 5), a double step transition evidences the formation of a vortex, and therefore the boundary of applicability of our SWbased analysis of the six-fold anisotropy field. The cutoff dimension between macrospin and non-macrospinlike behaviors are therefore quite similar between simulations and experiments in the explored size range: 15nm thick, 200 nm wide for the first, 18 nm thick, 150 nm for the second. It may be argued that the zero- temperature nature of these computations lets the system enter shallow energy wells, such as the vortex configuration, whereas thermal fluctuations in the experiment leave them unnoticed. Within the explored dimensions, the volume of the nanostructure at the boundary can however be considered identical between simulations and experiment.

Discussion and Conclusions

Following this measurement, six-fold anisotropy fields were extracted from the MOKE data for triangles e=150-300 nm and t=6-15 nm, as summarized in Table 1. H 6 increases with thickness, ranging from 128 Oe to 706 Oe. Where they can be compared, the values agree well with the anisotropy fields obtained in thinner triangles using MFMA [4]; the thinnest sample presented in the present study for instance (6 nm thick, 300 nm wide) has an anisotropy field of H 6 = 128 ± 10 Oe, very comparable to the 107 Oe obtained by MFMA in 5 nm thick, 270 nm wide triangles. The method presented here however, gives an estimation of the limit of a macrospin description of configurational anisotropy. Finally, these anisotropy fields are comparable to experimental uniaxial anisotropy fields H un obtained in Permalloy ellipses [12], with the major difference that they are obtained in non-elongated structures via configurational anisotropy.

In first approximation, but mainly on the demagnetizing field H demag of the structure. Increasing t creates more charges at the structure edges, and decreasing the width (keeping a/b constant) brings these poles closer together, both leading to a larger H demag [13,14]. A monotonous increase of H un as a function of t/a is then expected. This is in part what is observed in the triangles: at fixed triangle width, the anisotropy field increases with t (Table 1). At fixed triangle thickness however, H 6 first increases when the width decreases, and peaks at around 200 nm before decreasing. Interestingly, this latter trend with width is not observed for anisotropy fields extracted from zero-temperature micromagnetic simulations (Fig. 2) following the procedure described in section 4.3: the H 6 fields are in this case found to increase monotonously with increasing thickness and decreasing structure size. For the smaller structures (e.g 9 nm thick, 150nm wide) the anisotropy energy U = 2M S VH 6 /n 2 deduced from the experimental anisotropy fields amounts to about 10k B T (volume V of the triangle, symmetry order n of the structure [4]). Thermal fluctuations then become sufficient to overcome the anisotropy barrier, and H 6 starts to decrease as the structure is brought closer to the superparamagnetic regime. Note that configurational anisotropy is weak enough to let us observe the ferromagnetism-superparamagnetism transition at larger dimensions than in Permalloy ellipses [15].

In conclusion, we have devised a simple experimental method to extract the six-fold anisotropy fields of Permalloy triangles presenting configurational anisotropy. Using a large rotating field to impose a particular 60 • transition, we created conditions to compare our data to a Stoner-Wohlfarth model. Moreover, using micromagnetic simulations, we proposed a straightforward test to assess the validity of this approach for triangles of any given dimensions. Finally, we highlighted the specificity of configurational anisotropy in which the usual thickness over width dependence of the anisotropy field is not fully observed due to the proximity of the superparamagnetic regime. Equipped with this fundamental parameter characterizing the magnetic anisotropy, further work will focus on assessing the stray field interaction between these structures, in order to compare them to typical values of nanomagnets currently used in data storage technologies.

Figure 1 :

 1 Figure 1: (a) Stoner-Wohlfarth switching astroid calculated for a sixfold anisotropy term K 6 = 12000 J/m 3 , with the transitions ±60 • , ±120 • or 180 • . (b) Evolution of the envelope of the astroid for K 6 varying from 2000 to 12000 J/m 3 . (c) Energy landscape for fields of increasing amplitude applied along 30 • , leading to -60 • or -120 • transitions.

Figure 2 :

 2 Figure 2: (a) and (b) [color on-line] Numerical simulations of the effect of a counterclockwise rotating field of constant amplitude on a 310nm wide triangle. The projection along x of the magnetization is plotted with the corresponding field direction (thick red arrow) and global magnetization orientation (white arrow). (c) Experimental longitudinal Kerr magnetization plotted against the x component of a counterclockwise rotating field H rot = 324 Oe (sweeping direction given by the black arrows), for an array of 300 nm wide, 9 nm thick triangles.

Figure 3 :

 3 Figure 3: [color on-line] Polar plot of the rotating field amplitude versus the angle at which occurs a 60 • rotation of the global magnetization in 300 nm wide, 9 nm thick triangles. Clockwise (resp. counterclockwise) fields yield the -60 • transitions, empty symbols (resp. +60 • , closed symbols). The solid line is a fit to the Stoner-Wohlfarth calculation of the 0 -60 • transition, with K 6 = 10560 J/m 3 .

Figure 4 :

 4 Figure 4: [color on-line] Micromagnetic simulation of the switching astroids. The field H is saturated along -φ, and ramped back up along +φ (large empty red arrow), with the micromagnetic configuration evolving from (i) to (iii). Configuration (ii) is only observed for t=15 nm. The switching fields are plotted versus φ. (a) 300 nm wide, 10 nm thick triangle: abrupt switching of M for H reversed along φ = 0 • . (b) 200nm wide, 15nm thick triangle: a new, low-field transition appears when increasing the thickness. (c) Spin configurations corresponding to the φ = 0 • reversal for (a) and (b), with the global magnetization direction given by the white arrow, and the field direction by the full red arrow.

Figure 5 :

 5 Figure 5: [color on-line] Longitudinal Kerr remanence loops for the field applied along the base of the triangles: 15 nm thick, 300 nm wide (open symbols), or 18 nm thick, 150 nm wide (full line). The double-step in the latter evidences the passage through a vortex state, and therefore the boundary of the macrospin behavior.

Table 1 :

 1 elliptical or rectangular platelets (area a * b, thickness t), H un does not depend on exchange energy in Six-fold anisotropy fields H 6 extracted from rotating field loops, using H 6 = 2K 6 /M s .

	Thickness (nm) Width (nm) H 6 (Oe)
	6	300	128
		150	161
	9	200	317
		300	264
		150	561
	15	200	706
		300	670


Guillaume Mercier
email: guillaume.mercier@labri.fr

Jérôme Clet-Ortega
email: jerome.clet-ortega@labri.fr

Towards an Efficient Process Placement Policy for MPI Applications in Multicore Environments

Keywords: Message-Passing, Multicore architectures, Process placement

This paper presents a method to efficiently place MPI processes on multicore machines. Since MPI implementations often feature efficient supports for both shared-memory and network communication, an adequate placement policy is a crucial step to improve applications performance. As a case study, we show the results obtained for several NAS computing kernels and explain how the policy influences overall performance. In particular, we found out that a policy merely increasing the intranode communication ratio is not enough and that cache utilization is also an influential factor. A more sophisticated policy (eg. one taking into account the architecture's memory structure) is required to observe performance improvements.

Introduction

In the last decade, parallel computer architectures have evolved dramatically. Clusters have shifted from an assembly of uniprocessor machines interconnected by a single network to a complex and highly hierarchical structure. Nodes are now composed of several multicore processors sharing memory banks physically scattered across the node. Major CPU manufacturers like AMD (HyperTransport) and Intel (Quick Path Interconnect) follow this trend. The memory access time depends on the location of both the core and the memory bank. This is often designed as the Non-Uniform Memory Access (NUMA) effect. The memory hierarchy is also more complex due to the increase of cache levels. This sharing of memory resources depends on the CPU architecture and differs from a manufacturer to another. As for the interconnection network, multirail systems where several high-speed NICS are connected to a node are sometimes also encountered.

A real challenge for a parallel applications is to exploit such architectures at their full potential. In order to achieve the best performance, many factors must be taken into consideration and studied. The first one is to make use of an implementation of the MPI specifications [START_REF]Message Passing Interface Forum: MPI-2: Extensions to the message-passing interface[END_REF] able to efficiently take advantage of a multicore environment. Whilst the MPI standard is architecture-independent, it is an implementation's task to bridge the gap between the hardware's performance and the application's. Indeed, recent MPI-2 implementations such as Open MPI [START_REF] Gabriel | Open MPI: Goals, concept, and design of a next generation MPI implementation[END_REF] or MPICH2 [3] fulfill this purpose and offer a very satisfactory performance level on multicore architectures. However, in order for an MPI implementation to fully exploit the underlying hardware, the MPI application processes have to be placed carefully on the cores of the target architecture. This placement policy has to be defined by both the application's communication pattern and the hardware's characteristics. For instance, if some application processes communicate more frequently than others, they should be regrouped and placed on the same multicore node. By doing so, the amount of intranode communication increases and the application global performance will improve since intranode communication (shared memory) is faster than internode communication (network).

In this paper we expose the method and software tools we employed to allow an MPI application to better take advantage of a multicore environment. We will show a performance improvement not due to modifications of the MPI implementation itself but rather due to a relevent process placement. The rest of this paper is organized as follows: Section 2 describes how process placement is determined and on which set of tools and algorithms it relies. Experimentals results are presented in Section 3 with performance figures for some of the NAS computing kernels. Section 4 lists related works in this area of MPI process placement and discusses some issues raised by information gathering. Section 5 concludes this paper and opens future perspectives.

Computation of a relevant MPI process placement

In order to place the MPI processes in a relevant fashion, we have to gather information about the target architecture and the application's communication pattern. Once both are available we analyze them and determine the best possible placement. The criterion choice should follow a user-defined strategy. However, the current scope of this work does not encompass all MPI applications.

Hypotheses about MPI applications and their execution environment

In the rest of this paper we consider static MPI applications. By static, we mean that the application does not use any of the dynamic processes features offered by MPI-2. We also exclude hybrid MPI applications that rely on multithreading features (such as OpenMP directives). The number of computing entities (threads and processes) is therefore guaranteed to remain constant during an application's execution. We also consider the target machine to be fully dedicated to the MPI application. All cores are usable by the MPI processes with the restriction that only a single MPI process runs on a given core. As a consequence of these points, a static mapping between the MPI processes and the CPU cores can be computed before launching the application. This placement will not need to be modified during the application's execution.

Gathering the hardware's information

As previously explained, clusters of NUMA nodes are hierarchically structured. For instance, figure 1 shows the architecture of an AMD Opteron-based compute node. This compute node is composed of four dies with two cores each. Each die possesses a set of main memory banks attached to it. A core features its own Level 1 cache (not shown) and Level 2 cache, not shared with the other core on the same die. A core located on a die can access the main memory of any other die but the access time increases as the physical distance between the core and the memory bank lengthens. A network interface card (NIC) can be attached to a bus connected to Die #0 and Die #1 . Since the cores located on these dies are physically closer to the I/O bus, one might expect faster network transactions for processes mapped on these cores. Some tools can provide us with the needed information (e.g libtopology or Portable Linux Processor Affinity [START_REF] Squyres | Portable Linux Processor Affinity[END_REF] available for Linux). However, such tools are not portable and accurate enough over a wide spectrum of operating systems. To gather the hardware's information, we used a topology discovery mechanism implemented in the PM 2 runtime system [START_REF] Namyst | Utilisation des processus légers pour le calcul parallèle distribué : l'approche PM2. Calculateurs Parallèles[END_REF]. This feature is employed to help the scheduling of threads on multicore nodes [START_REF] Thibault | Building Portable Thread Schedulers for Hierarchical Multiprocessors: the BubbleSched Framework[END_REF]. Currently, this mechanism only deals with CPUs and does not deliver information about networks and I/O buses. For instance, when applied to the Opteron node as depicted by figure 1 Fig. 2. Hardware information generated by PM 2 's topology discovery mechanism.

Mem Mem

We generate from this ouptut a data structure that other tools will use to compute the placement. This data structure is a graph with vertices representing the CPU cores and weighted edges. This graph is complete and non-oriented. The weight affected to each edge increases as more elements of the memory hierarchy are shared between cores. It is also affected by the NUMA effect.

Figure 3 shows the various weight values chosen for the Opteron compute node (figure 1). We decided to put the largest weight for cores on the same die because they directly share memory banks. The next largest weight value corresponds to NUMA effects: for instance core #0 is able to access the memory attached to Die #2 and Die #1 faster than the memory attached to Die #3. When the machine is made of several compute nodes, we build a larger graph and affect a larger weight for cores located within the same compute node. The smallest weight values thus correspond to communication using the network.

Collecting the application's communication pattern data

The second piece of information deals with the application's communication pattern. Each application possesses its own pattern influenced by the number of participating processes. Our chosen characterization criterion for this pattern is the amount of MPI data exchanged between processes. Therefore, we need to compute this amount for each pair of processes in the application.

Several sophisticated tools are provided for MPI application tracing and analysis, such as the MPI Parallel Environment (MPE). However, they do not provide all the necessary information. For instance, MPE is able to trace all calls to MPI routines made by the application. By analyzing such a trace and focusing on the point-to-point calls, we can have hints about the communication pattern. This approach is simple and require to configure the MPI implementation with MPE support and to link the application with the appropriate libraries. This is limited by two factors: first, the trace generated can be potentially very large and second, the amount of data exchanged in collective operations is not taken into account.

As a consequence we found simpler and more accurate to modify an MPI implementation to collect the desired information. By modifying directly an implementation, we reduce drastically the trace size and take into account collective communication operations. In order to get a generic (that is, not implementationspecific) information, we trace only the size of MPI user data exchanged between processes. All costs induced by the implementation's internal protocols are not counted. As for the hardware's data, we represent this communication pattern with a complete, non-oriented graph with weighted edges. In this case, the weight value increases as the amount of data exchanged between MPI processes grows. An example is depicted by figure 4 (i,j) represents the amount of communication between processes #i and #j).

Mapping a MPI process rank to a CPU core number

The final step is to extract an embedding of the application's graph from the target machine's graph. We use the Scotch software [START_REF] Pellegrini | Scotch and LibScotch 5.1 User's Guide[END_REF] to solve this NP graph problem. Scotch applies graph theory, with a divide and conquer approach, to scientific computing problems such as graph and mesh partitioning, static mapping, and sparse matrix ordering. In our case, we use the ability to construct a static mapping. Scotch implements dual recursive bipartitioning algorithms to perform this task [START_REF] Pellegrini | Static Mapping by Dual Recursive Bipartitioning of Process and Architecture Graphs[END_REF] and computes static mappings for graphs larger than 2 32 vertices. This ensures that we can create a mapping for all MPI applications, regardless of their size. For instance, table 1 gives the static mapping computed by Scotch in the case of the NAS lu.B.8 benchmark (figure 4) launched on the Opteron compute node (figure 3). Using this static mapping, we finally generate a specific command line fully customized for each (Application,Target Machine) couple. Practically, each MPI process is affected to its dedicated core with the numactl command.

Evaluation: a case study with NAS computing kernels

In this section, we present some results obtained by applying the method previously described on several NAS computing kernels. We first describe our experimental environment, then show the results and finally comment on them.

Experimental environment and NAS benchmarks choice

We carried out experiments on a 10-nodes cluster called Borderline, part of the Grid5000 testbed [START_REF] Bolze | Grid'5000: a large scale and highly reconfigurable experimental Grid testbed[END_REF]. Borderline nodes are similar to the exemple depicted by figure 1: each node features four dies (a 2.6 GHz AMD Opteron 2218) with two cores each. A core possesses its own Level 2 cache (1 MBytes) not shared with the other core on the same die. The total amount of memory is 32 GBytes per node (8 GBytes per die). A Myrinet 10G NIC is attached on a bus to Die #0. The Linux kernel installed is 2.6.22.1 (SMP version).

We benchmark some of the NAS computing kernels in order to assess the relevance of our placement method and policy. We first make a run of all NAS benchmarks to gather their respective communication pattern data, as described in section 2.3. Since our placement technique is currently based only on the global amount of data exchanged, we only run tests where this amount is significant: 64-processes jobs of classes C and D. Likewise, we test only applications with irregular communication patterns. Indeed, they are likely to be the most influenced by the placement of processes in a multicore environment. Among all NAS kernels, only BT, CG, LU, MG and SP have an irregular communication pattern. In the case of CG, some processes do not communicate with others. But when communication occurs, the global amount of data exchanged for each pair of processes is of same magnitude. In the case of BT, LU, MG and SP these amounts differ dramatically from one pair to the other.

NAS performance results for two MPI implementations

In order to confirm that the method used and the data gathered are implementation-independent, we expose performance comparisons for two different MPI implementations. This first one is MPICH2-Nemesis, which relies on a very efficient intranode communication system using shared-memory [START_REF] Buntinas | Implementation and Evaluation of Shared-Memory Communication and Synchronization Operations in MPICH2 using the Nemesis Communication Subsystem[END_REF]. We configured Nemesis to use its most recent MX support (available in the MPICH2 1.1 release). The other MPI implementation is MPICH2-MX [START_REF] Myricom | MPICH2-MX[END_REF] designed and implemented by Myricom. This implementation also features an efficient intranode communication system based on an in-kernel mechanism. Table 2, figures 5 and 6 show the NAS performance for two different placement policies. The first placement policy is a simple Round-Robin-type policy where the process number i is executed on the node number n where n ≡ i mod 8 (in our case each node features eight cores). This type of placement is commonly used when launching MPI applications without any knowledge of their communication patterns. With this scheme, the operating system chooses on which core a process is executed. The other policy, called Placed is the one determined thanks to the method described in section 2. In the case of MPICH2-Nemesis:MX, there is a significant performance gap between the Placed policy and the Round-Robin policy. This result applies to both C and D classes. As our results indicate, this improvement for Classes C and D kernels is roughly of 25%. In the case of Myricom's MPICH2-MX, this gap is even larger: execution times are 34% faster when our Placed policy is enforced. One might note that in the case of the CG kernel (class C), performance is much better with our custom placement than with the Round-Robin one.

NAS performance and placement policies analysis

The Placed policy regroups MPI processes on the same node as much as possible.

The ratio of intranode communication (using shared-memory) versus internode communication (using the Myrinet network) increases as Table 3 shows. However, questions remain: is this performance gap solely due to this increase of intranode communication ratio in the application? Do others factors -such as memory hierarchy and structure -influence performance too? In order to answer

BT CG LU MG SP C D C D C D C D C D Round-
Robin 32% 33% 0% 0% 46% 46% 39% 50% 31% 33% Placed 53% 54% 78% 78% 70% 70% 55% 66% 52% 54%

Table 3. Overall intranode communication ratios in NAS kernels for two different placement policies.

these questions, we run the NAS kernels with three other placement policies. The first new policy is the same as the Round Robin policy, except that each process is bound to a particular core with the numactl command. In this way, the operating system cannot change a process location when scheduling occurs. The Level 2 cache utilization is therefore better. In this case, the intranode communication ratio is not improved. The second new policy regroups the MPI processes on the nodes as much as possible (like in the Placed policy) but instead of binding each MPI process on its dedicated core, we let the operating system choose the placement (like in the regular Round-Robin policy). By doing so, the overall ratio of intranode communication versus internode communication is increased but we do not take anymore into account factors such as NUMA effects, memory hierarchy or die sharing by processes. Cache utilization will be negatively impacted when processes are placed according to this policy. The last policy also regroups process on the compute nodes but we force the intranode placement to be suboptimal: in this case we place the processes that communicate the most on opposite dies in the compute node. This last policy and the Regroup policy share a common point in that process placement within the node is very poor. The difference between the two is that with the Placed Core Reverse policy, cache utilization is better since processes are pinned to a core. Figure 7 shows the results: both Round-Robin and Placed policies are the same as in section 3.2, while Round-Robin Core Binding, Regroup and Placed Core Reverse represent the new placement policies described above. The results indicate that merely improving the intranode communication ratio is not enough to deliver better performance. Increasing this ratio without taking into account cache utilization leads to suboptimal results. Globally, our Placed policy delivers slightly better results than the others policies. The CG kernel is the benchmark for which the performance improvement is the most noticeable. These results suggest that cache utilization has a great influence on NAS performance. and advocates for a placement policy that takes architectural factors into consideration. However, we think that the Opteron compute nodes we used have a too small NUMA effect for it to impact applications performance. Also, the results obtained with the Placed Core Reverse policy show that despite the increase of traffic on the memory bus (induced by this policy), the impact on performance is limited on the compute nodes we used. Last, Level 2 cache is not shared between cores on the same die as explained in section 2.2. Those three aspects lead us to think that our placement policy could perform even better on nodes featuring more cores, with substantial NUMA effects, and where cache is effectively shared between cores. Finally, since we do not know the theoretical achievable performance of this particular cluster, we cannot assess how close to the optimal placement we actually are.

Related Works and Discussion

An adequate process placement is a necessity to fully take advantage of multicore environments. The MPI standard itself offers a set of routines that allow the applications developers to create and manipulate topologies. Using such topologies, MPI processes could be placed according application's communication pattern. However, only a subset of existing MPI applications makes use of these topology creation and management features. Also, as [START_REF] Larsson | Implementing the MPI process topology mechanism[END_REF] points out, not all MPI implementations support these efficiently. At last, this mapping between the virtual topology created and the physical topology (the issue on which we focus in this paper) is outside the scope of the MPI standard. Our approach -deploying the MPI processes according to a matching between the application's communication pattern and the machine's architecture -is more generic and can benefit to any MPI application. Actually, this approach is not MPI-dependent, but rather message-passing-dependent. It can be applied even to non-MPI applications as long as the necessary information is collectable.

Placing MPI processes according to the underlying target architecture using graph theory has already been explored. Several vendor MPI implementations, such as the ones provided by Hewlett-Packard [START_REF] Solt | A profile based approach for topology aware MPI rank placement[END_REF] or by IBM (according to [START_REF] Duesterwald | Method and System for Optimizing Communication in MPI Programs for an Execution Environment[END_REF]) make use of such mechanism. [START_REF] Larsson | Implementing the MPI process topology mechanism[END_REF] also formalizes the problem with graphs. In these cases, however, the algorithm computing the final mapping is MPI implementation-specific. In our framework, we rely on an external piece of software (Scotch) fully tailored for graph computations. This ensures both performance and scalability. Scotch is indeed able to work on very large classes of graphs, larger that the maximum number of MPI processes present in any MPI application.

As explained in section 2.2, architecture's information is gathered thanks to PM 2 's topology discovery mechanism. But this feature is fully embedded in this software stack and rather unconvenient to use. Such a topology discovery feature would be very useful for many other applications, especially MPI process managers. Hybrid MPI+OpenMP applications could also take advantage from it [START_REF] Rabenseifner | Hybrid MPI/OpenMP Parallel Programming on Clusters of Multi-Core SMP Nodes[END_REF]. There is a clear need for a portable and accurate tool dedicated to topology discovery. We already started to work on this specific point.

The last issue we would like to address regards the gathering of applications communication patterns. In order to get the needed information to create the mapping, we have to execute a prior run of the applications compiled with a modified MPI implementation fitting our needs. Other works in this field use the same coarse scheme. They also emphasize that tracing tools could be enhanced to address the specific issue of getting the amount of data exchanged between two given processes, as HP's Light Weight Instrumentation [START_REF] Solt | A profile based approach for topology aware MPI rank placement[END_REF] does. The necessary information could be an estimate of the communication flows that would allow us to rank pairs of processes accordingly. However, the ability to perform this prior run does not systematically exists. Other solutions should be investigated. For instance, relying on an application's programmer's knowledge is also possible but once again far from always possible. Could this kind of information be computed at compile time? Would it be possible to pass an option to the mpicc compiler that would automatically generate the customized command line? These are issues we would like to address in the near future.

Conclusion and Future Work

In this paper, we expose a method that leads an MPI application to better exploit its target architecture, especially complex and hierarchical multicore environments. This method relies on gathering information about the underlying hardware and the application's communication pattern. This information is then used to create a mapping between MPI process ranks and each node's core numbers. Finally an application-specific command line is generated. Our method uses free, open-source software and is not tightly integrated within a particular MPI implementation. Being able to place the processes according to the marchine's topology increases performance. We analyzed several placement policies and found out that merely increasing the intranode communication ratio in an application is not enough to deliver more performance: an adequate cache utilization is also mandatory to enhance performance. The current results are mitigated: our sophisticated placement policy improves MPI performance applications but should yield even better results on more complex architectures. We look forward to run our experiments on a different (and more complex) class of hardware.

Possible future works are numerous: as we previously stated, we began to work on a software tool that would allow us to easily extract hardware's information. Integrating information about I/O buses or even GPUs should also be considered. We would like to relax some constraints about the type of MPI applications falling into the scope of this work. We do consider static MPI applications, that is, applications where the amount of computing entities remains constant throughout the execution. This excludes both applications spawning new MPI processes and multithreaded ones. Since hybrid programming, mixing message-passing and multithreading, is considered as a possible way to better program multicore architectures, we plan to address the issue of MPI processes placement when OpenMP parallel sections appear in MPI processes. Another interesting direction would be to refine the application data regarding its communication pattern. For now, we did only consider a spatial pattern. But what about the temporal pattern? Indeed, the data exchanges occuring between a given pair of processes may vary during the application's execution. In order to take this phenomenon into account we will have to isolate application time slices and remap the MPI processes during such time slices. What granularity for slices would be the most beneficial to performance? Also, we would have to modify the mapping during execution. For intranode communication this task is easy but for internode communication we would have to migrate processes from one node to the other. Using virtual machines in this context might be a way to implement this. Also, some other influential factors such as contention on the memory and I/O buses could be taken more into consideration. This would lead to an even more refined placement policy but requires to comprehend thoroughly the target application. Generating MPI-implementation dependent information could also be considered in order to increase the policy's accuracy. Last, we plan to investigate the feasability of gathering the application's information at compile time.

Fig. 1 .

 1 Fig. 1. An exemple of hierarchical compute node: an 8-cores Opteron

Fig. 3 .

 3 Fig. 3. Matrix representation of the hardware's information from figure 2. (Machine(i,j) represents the potential of communication between cores #i and #j).

Fig. 4 .

 4 Fig. 4. Matrix representation of NAS LU (class B, 8 processes) communication pattern. (Application(i,j) represents the amount of communication between processes #i and #j).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Process placements comparison for MPICH2-Nemesis:MX.

Fig. 7 .

 7 Fig. 7. Process placements comparison for MPICH2-Nemesis:MX.

  , the PM 2 topology discovery mechanism outputs the information displayed by figure2.

	Machine:
	NUMANode + Die: Node#0(8GB) Die#0
	L2Cache + Core + L1Cache + SMTproc : L2#0(1MB) Core#0 L1#0(64kB) CPU#0
	L2Cache + Core + L1Cache + SMTproc : L2#4(1MB) Core#1 L1#4(64kB) CPU#4
	NUMANode + Die: Node#1(8GB) Die#1
	L2Cache + Core + L1Cache + SMTproc : L2#1(1MB) Core#0 L1#1(64kB) CPU#1
	L2Cache + Core + L1Cache + SMTproc : L2#5(1MB) Core#1 L1#5(64kB) CPU#5
	NUMANode + Die: Node#2(8GB) Die#2
	L2Cache + Core + L1Cache + SMTproc : L2#2(1MB) Core#0 L1#2(64kB) CPU#2
	L2Cache + Core + L1Cache + SMTproc : L2#6(1MB) Core#1 L1#6(64kB) CPU#6
	NUMANode + Die: Node#3(8GB) Die#3
	L2Cache + Core + L1Cache + SMTproc : L2#3(1MB) Core#0 L1#3(64kB) CPU#3
	L2Cache + Core + L1Cache + SMTproc : L2#7(1MB) Core#1 L1#7(64kB) CPU#7

Table 1 .

 1 Resulting mapping for application NAS lu.B.8 on the Opteron compute node.

Table 2 .

 2 MPICH2-Nemesis:MX execution times in seconds for two different placement policies.

		BT		CG		LU			MG	SP
		C	D	C	D	C	D	C	D	C	D
	Round-Robin 51.6 1038.6 23.9 1177.5 45.5 1356.5 5.6 119.2 78.6 2015.63
	Placed	45.6 851.7 15.6 848.4 33.6 938.3 3.7	85.3	60.2 1386.8

Acknowledgments

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr). The authors wish to thank Sébastien Fourestier, our Scotch Guru.


Roger Botton

Dominique Cosserat

Souhila Poncin

Gabriel Wild

   A SIMPLE GAS-LIQUID MASS TRANSFER JET SYSTEM, 

Keywords: Gas-liquid reactors, impinging jet, energy efficiency, bubble columns, stirred tank. 

An original gas-liquid contacting system is proposed, consisting of a pump, an orifice, a vertical tube coaxial to the orifice and an impinging plate. The pump generates a downward vertical liquid jet through the orifice situated above the gas-liquid dispersion level. The two phase jet is directed towards an impinging plate near the bottom of the tank and dispersed in the volume of the liquid. Liquid is withdrawn below the impinging plate and recycled. This reactor may be used for gas-liquid reactions (ie hydrogenations) and also to mix liquids, to disperse particles, to oxygenate waste water etc…. Performances and design rules of this equipment are proposed. Then, the results are compared to performances of bubble columns, stirred tanks, and other academic and industrial jet systems. It is shown that, at a given energy dissipation, this system yields much higher mass transfer densities than a classical stirred tank provided with a Rushton turbine. Finally some suggestions about mass transfer mechanisms and efficiency of dissipated power are given.

INTRODUCTION

 Chemical and biochemical industries need simple high performance gas-liquid systems to achieve high productivity and selectivity reactors. The purpose of this article is to propose a very simple reactor achieving this goal. It consists only of an orifice, a tube, a plate and a pump. In the present system, mechanical energy is provided by the liquid pump exclusively, while in bubble columns it is provided by the gas compressor, and in stirred tank reactors by the stirrer engine and the gas compressor. The aim of this article is to provide simple and reliable design rules for the jet device presented here, and to compare the energy efficiency in terms of k L a values with that of other devices. 

EXPERIMENTAL SET-UP AND OPERATING CONDITIONS

2.1 Description of the impinging jet device used here.  Fig. 1 shows an example of an experimental set-up. It consists of a tank, a pump, a nozzle with an orifice diameter d o , a perforated tube D (coaxial to the orifice) above and below the dispersion surface, and of an impinging plate P. The gas is introduced in the tank by the top. It is entrained downwards in tube D by the liquid jet provided by the nozzle; further gas and gas-liquid mixture is sucked in through the holes in the wall of D and is entrained downwards. The gas-liquid jet is then carried through tube D and impinges on plate P; the impact causes the gas to be dispersed in the whole the tank. Above the plate, the dispersion of bubbles has a uniform aspect, while below the plate, with an appropriate design, no bubbles are observed. The liquid is recycled from below the plate through the pump. The dimensions indicated in Fig. 1 correspond to optimal results according to experiments made in industrial scale equipment. If necessary, several tanks can be used to obtain the gas and liquid flow rates needed by the reaction. If, for safety reasons (e.g. hydrogenations), the exit of liquid at the bottom of tank is not wanted, it is possible to place it at the top, using a tube parallel or coaxial to tube D and the nozzle d 0 . The advantages, compared with a stirred tank are the absence of moving parts (eliminating sealing problems) and the possibility, if a heat exchanger is required, to install it outside of the tank on the liquid recycling line.

 

Operating conditions and measuring techniques using the impinging jet device.



The gas used in the experiments was air at room temperature and pressure, the liquid was either tap water or a water-sulfite solution in case of mass transfer or interfacial area measurements. Several jet systems were tested with air and water using a nozzle diameter of d o =0.010m and tanks of diameters T=0.25 and 0.48m where dispersion levels were adjusted by overflow. Tests have consisted of measuring the total bubble volume and mass transfer rate generated by a given energy input. In the case of a coalescing gasliquid system, under the conditions investigated, no gas entrainment below the plate was observed, so that an effective recycle of the liquid was possible. The effect of the orifice diameter, (always with air and water) was investigated using the following two plants. The first one, with a nozzle orifice diameter d o =0.012m, was called "optimized system N°1" and the second, with d o =0.024m, "optimized system N°2". The other dimensions are those indicated in Fig. 1. During all these studies, the following variables were generally tested: jet velocity V o (4 to 12 m.s -1 ), nozzle orifice diameter d o (0.009, 0.010, 0.012, 0.024 m), tube diameter D (0.025, 0.033, 0.045, 0,057, 0.067 m with or without immersed holes to recycle dispersion), tank diameter T (0.25, 0.48 m), dispersion height H E fixed by t overflow (0.5, 0.85, 1.85, 2.35 m). Some experiments were made with a sodium sulfite solution to investigate the effect of a coalescence inhibiting system and to determine the volumetric mass transfer coefficient k L a or the gas-liquid interfacial area by the chemical technique.

The measuring techniques are quite classical: the gas volume in the tank was determined by stopping the liquid flow and by measuring the liquid level before and after gassing. The liquid flow rate was determined using a rotameter, the induced gas flow rate by a Pitot tube at the gas entrance and the pressure in the nozzle by a manometer. The global liquid side mass transfer capacity Q T (Q T =k L a E V E, =k L a L V L ) (m 3 .s -1 ) was determined by using the classical dynamic reoxygenation technique or by the slow oxidation of sufite reaction technique. The interfacial area by using the pseudo-mth order rapid reaction technique (oxidation of sulfite catalyzed by cobalt sulfate).

Experimental results

Contrarily to the classical gas-liquid mixed tanks and bubble columns, the induced gas flowrate depends on the liquid velocity at the outlet of the nozzle (ie the liquid flowrate and orifice diameter d o of the nozzle). Only some hydrodynamics and mass transfer results obtained with the optimized reactors will be presented here. The effect of operating conditions on hydrodynamic (ie gas induced flow rate and gas holdup) and mass transfer parameters (interfacial area and volumetric liquid side mass transfer k L a) has been studied in order to establish correlations between design parameters and operating conditions for design purpose.

Hydrodynamic parameters :

The effect of the liquid velocity at the nozzle outlet on the induced gas flowrate Q G obtained with both (N°1 and N°2) optimized jet impinging reactors is represented on figure 2. One can observe the logical increase of Q G with V o due to the decrease of pressure as V o increases at the outlet of the nozzle. Furthermore, no effect of reactor geometry is obtained, which is probably due to the fact that the optimized reactors are almost geometrically similar. The evolution of the gas holdup with V o (figure not shown) is similar to that of Q G .

Mass transfer parameters:

Results of the volumetric gas-liquid mass transfer coefficient k L a =Q T /V L obtained from the dynamic method using air/water system are reported on figure 3 as a function of the gas to the liquid volume ratio. As expected, considering the increase of the gas induced flowrate and the gas holdup with the liquid velocity V o at the nozzle outlet, the values of k L a also increase with V B /V L regardless of the geometry of the optimized impinging jet reactor as in the case of the gas induced flow rate results shown before.

Other results of k L a have also been obtained for different emulsion heights and from the slow sulfite oxidation reaction technique. The volumetric mass transfer coefficient k L a generally increases with a decrease of the emulsion height due to a the increase of the gas induced flowrate. In presence of coalescence inhibiting system which is the case of the air/sulfite system, k L a values are much higher than in presence of air/water system due to higher gas suction rates and gas holdups (not represented here). This result is commonly reported in literature for other types of gas-liquid reactors (stirred tanks, bubble columns etc…). [START_REF] Van 't Riet | Mass Transfer in Fermentation[END_REF], [START_REF] Midoux | Réacteurs gaz liquide à cuves agitées[END_REF], [START_REF] Botton | Operating zone and scale up of mechanically stirred gasliquid reactors[END_REF]; (5) Impingement of two aerated jets. [START_REF] Gaddis | Hydrodynamics and mass transfer in a impinging-stream loop reactor aerated from the top[END_REF]; (6) Plume dispersion confined into a tube. [START_REF] Evans | Performance of confined plunging liquid jet bubble column as a gas-liquid reactor[END_REF] The volumetric gas-liquid mass transfer coefficients (k L a = Q T /V L ) obtained with the impinging jet reactor was finally compared to those determined in the studied bubble columns and stirred tanks and also to those published for other reactor geometries. For this purpose, k L a values have been represented in figure 4 as a function of the specific power input P W /V L . Figure 4 shows that the impinging jet system is more efficient than bubble columns and stirred tanks for coalescence inhibiting liquids. For coalescent liquids, it is more efficient than the classical gas-liquid stirred tanks provided with a Rushton turbine and avoids the sealing problems commonly encountered in stirred tanks; it is almost as efficient as bubble columns but easier to scale up and to simulate than bubble columns. Also indicated are values obtained with different other contacting devices: higher volumetric mass transfer coefficients can be achieved (e.g. by the impingement of two aerated jets, [START_REF] Gaddis | Hydrodynamics and mass transfer in a impinging-stream loop reactor aerated from the top[END_REF], but at the expense of higher mechanical energy dissipation.

Fig. 1 .

 1 Fig. 1. Description of an optimized gas-liquid impinging jet system

Fig. 2 .

 2 Fig.2. Gas rate induced by the liquid jet for both optimized reactors :N°1 (d o =12 mm): open symbols; N°2 (d o =24 mm) black symbols. Results obtained with Air/water system.

Fig. 3 .

 3 Fig. 3. Evolution of k L a (= Q T /V L ) with the gas to liquid volume ratio for both optimized reactors: N°1 (d o =12 mm): open symbols; N°2 (d o =24 mm): black symbols. Air/ water system.Comparison of the impinging jet reactor with other classical reactors:Hydrodynamic and mass transfer experiments have also been conducted in bubble columns and in stirred tanks provided with a Rushton turbine in order to compare k L a values at a given specific dissipated power. The chosen reactors had quite similar heights and diameters as those of the impinging jet reactors.

Fig. 4 .

 4 Fig. 4. Liquid-side mass transfer versus specific dissipated power. Continuous lines: air/water system; doted lines: air/sodium sulphite solution; (1A, 1B) .Studied impinging reactors for power at entrance and exit of orifice; (2) Bubble columns (This work); (3h, 3b) Plume dispersions with down and up aerated jets. Tojo et al (1982); (4) Mechanically stirred tanks[START_REF] Van 't Riet | Mass Transfer in Fermentation[END_REF],[START_REF] Midoux | Réacteurs gaz liquide à cuves agitées[END_REF],[START_REF] Botton | Operating zone and scale up of mechanically stirred gasliquid reactors[END_REF]; (5) Impingement of two aerated jets.[START_REF] Gaddis | Hydrodynamics and mass transfer in a impinging-stream loop reactor aerated from the top[END_REF]; (6) Plume dispersion confined into a tube.[START_REF] Evans | Performance of confined plunging liquid jet bubble column as a gas-liquid reactor[END_REF]


Guillaume Collet

Rumen Andonov

Nicola Yanev

Jean-Franois Gibrat

Keywords: Mixed integer programming, combinatorial optimization, protein threading problem, protein structure alignment

During the last decade, significant progresses have been made in solving the Protein Threading Problem (PTP). However, all previous approaches to PTP only perform global sequence-structure alignment. This obvious limitation is in clear contrast with the "world of sequences", where local sequence-sequence alignments are widely used to find functionally important regions in families of proteins. This paper presents a novel approach to PTP which allows to align a part of a protein structure onto a protein sequence in order to detect local similarities. We show experimentally that such local sequencestructure alignments improve the quality of the prediction. Our approach is based on Mixed Integer Programming (MIP) which has been shown to be very successful in this domain. We describe five MIP models for local sequencestructure alignments, compare and analyze their performances by using ILOG CPLEX 10 solver on a benchmark of proteins.

Local Protein Threading by Mixed Integer Programming

Résumé : Au cours des dernières décennies, le problème de la reconnaissance de repliements des protéines ou Protein Threading Problem (PTP) a connu de grandes avancées. Néanmoins, toutes les méthodes développées ne proposent que des alignements séquence-structure globaux. Or, les alignements séquenceséquence locaux ont montrés qu'ils permettaient de détecter des régions fonctionnelles dans des familles de protéines. Ce rapport présente une nouvelle approche de la reconnaissance des repliements qui permet d'aligner une partie d'une structure de protéine avec une partie d'une séquence de protéine afin de détecter des similarités locales. Nous montrons que cette approche, basée sur la programmation mixte en nombre entiers (MIP), améliore la qualité de la reconnaissance. Nous avons modélisé les alignements séquence-structure locaux par cinq modèles MIP que nous comparons et analysons grâce au logiciel CPLEX 10.0 sur un jeu de test de protéines.

Mots-clés :

Programmation mixte en nombre entiers, Optimisation combinatoire, Reconnaissance des repliements, Alignements séquence-structure de protéines

Introduction

To exploit the amount of new genomic data, the most important in silico methods are based on the concept of homology. The principle of homology-based analysis is to identify a homology relationship between a new protein and a protein whose function is known. Detecting homology relationships when the sequences are sufficiently similar is relatively straightforward. There exist efficient O(N 2 ) algorithms based on dynamic programming techniques to align protein sequences [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF][START_REF] Smith | Identification of Common Molecular subsequences[END_REF]. Suitable modifications of the fundamental algorithm enable the user to perform global, semi-global and local alignments (see Fig. 1).

Global alignment

Semi-global alignment

Local alignment

Figure 1: Types of alignments. When aligning protein sequences, it is important to be able to carry out different types of alignments according to the situation encountered. Top: When aligning proteins belonging to the same family, global alignments are used, i.e. gaps (illustrated here by dashed lines) at both ends of the alignments are penalized. Middle: On the other hand, when one looks for a protein domain in a longer sequence, semi-global alignments that allow the shorter sequence to be aligned (entirely) along the longer one are required. Gaps at both ends are not penalized in semi-global alignments. Bottom: The most general type of alignment is the local alignment, where only substrings of both sequences can be aligned (this might correspond to a common domain found in two proteins that are otherwise different).

For remote homologs, i.e. homologous proteins for which the sequences are no longer similar, sequence alignment methods fail. One, then, must resort to the good conservation of the 3-dimensional (3D) structures. In such a case, to figure out a homology relationship between 2 proteins, one aligns the sequence of the new protein with the 3D structures of known proteins. Such methods are called fold recognition methods or threading methods since this can be thought as "threading" the sequence through the 3D structure.

For a long time, threading methods using non local parameters (see section 2.1 for the definition of this term) suffered from the lack of a rigorous method capable of determining the sequence-structure alignment with the optimal score. They relied on heuristic techniques, e.g. stochastic techniques such as a Gibbs Monte Carlo [3]. R. Lathrop showed that, in the most general case, when variable length alignment gaps are allowed and pairwise amino acid interactions are considered in the score function (i.e. non local parameters), the problem of aligning a sequence onto a 3D structure is NP-hard [4]. The algorithm proposed by Lathrop & Smith was the first for finding the global alignment with the optimal score [START_REF] Lathrop | Global optimum protein threading with gapped alignment and empirical pair potentials[END_REF]. Since then, other methods have been developed that improved the efficiency of the global sequence-structure alignment algorithms [START_REF] Xu | RAPTOR: optimal protein threading by linear programming[END_REF][START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF][START_REF] Yanev | Lagrangian approaches for a class of matching problems in computational biology[END_REF][START_REF] Andonov | Recent Advances in Solving the Protein Threading Problem[END_REF].

This study focusses on the local sequence-structure alignment problem. We propose an approach based on integer programming which expands upon our previous results [START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF][START_REF] Andonov | Recent Advances in Solving the Protein Threading Problem[END_REF]. To the best of our knowledge, no attempt has been done previously for solving this problem. We show experimentally that such local sequence-structure alignments improve the quality of the prediction. This allows threading methods to cover the whole spectrum of alignment types needed to analyze homologous proteins.

Outline of the Protein Threading Problem

In this section we use classical notations and definitions concerning PTP as given in [START_REF] Lathrop | Global optimum protein threading with gapped alignment and empirical pair potentials[END_REF].

Definition of alignments

Query Sequence A query sequence is a string of length N over the 20-letter amino acid alphabet. This is the amino acid sequence of an unknown protein which must be aligned to structure templates from a given database.

Structure Template All current threading methods replace the 3D coordinates of the known structure by an abstract template description in terms of blocks, neighbor relationships, distances, environments, etc (see Fig 2). Blocks correspond to the most conserved parts of the structure, usually the secondary structure elements (SSEs : α-helices and β-strands). We consider that a structure template is an ordered set M of m segments or blocks. Block k has a fixed length of L k amino acids. Alignments The alignment of a sequence with a 3D structure can be described as positioning the blocks along the sequence. Such an alignment is called feasible if blocks preserve their original order and do not overlap. An alignment Figure 3: Example of alignment. Alignment of query sequence of length N = 24 and template containing m = 3 blocks of lengths L 1 = 6, L 2 = 4 and L 3 = 8. An arrow is drawn between 2 blocks if there is at least one position in each block which are in interaction in the 3D structure. The corresponding generalized contact map graph is given by G = (M, I), where M = 1, 2, 3 and [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF]3), (2, 3)} is completely determined by the starting positions of all blocks along the sequence. In fact, in the classical PTP, it is more convenient to use relative positions instead of absolute positions. If block k starts at the jth query character, its relative position is r k = j -k-1 i=1 L i . In this way the possible (relative) positions of each segment are between 1 and n = N -m k=1 L k +1. The set of possible alignments (feasible threadings) is

I = {(1, 2),
T = {(r 1 , . . . , r m ) | 1 ≤ r 1 ≤ • • • ≤ r m ≤ n}.
The cardinality of this set (the search space size of PTP) is given by

|T | = m+n-1 m
, which is a huge number even for small instances (for example, if m = 20 and n = 100 then |T | ≈ 2.5 × 10 22 ).

With this definition of the PTP, gaps are not allowed within blocks. They are confined to loops joining the blocks.

Network flow formulation

In order to develop an appropriate mathematical model, PTP has been reformulated in [START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF][START_REF] Yanev | Parallel Divide&Conquer Approach for the Protein Threading Problem[END_REF] as a network optimization problem. Let G = (V, E) be a m-partite graph with vertex set V and edge set E. The vertex set V is organized in n × m grid as shown in Fig 4 . The vertex (i, k) (where i stands for row, while k stands for column) corresponds to the relative position of the block k along the sequence. The graph G is called an alignment graph. To preserve blocks order, each edge ((i, k), (j, l)) ∈ E satisfies the condition i < j and k < l (i.e. all edges are southwest-northeast oriented). For this reason, the alignment graph can be also thought as a directed graph, and since we search for a southwest-northeast oriented path in this graph, there is no ambiguity when we sometimes indicate the edges of a given vertex as input or output edges.

A feasible path F in the alignment graph G is a set of exactly m vertices A cost C ik is associated to each vertex (i, k) as defined by the scoring function. C ik is a local score that depends only on the location of the block k at position i along the sequence. To each edge ((i, k), (j, l)) we associate a cost C ikjl as defined by the scoring function. C ikjl is a non local score that depends on the location of blocks k and l on positions i and j along the sequence. We say that an edge is activated by a feasible path if both ends are on the path. The subgraph induced by the activated edges of a feasible path is called an augmented path. The main result from [START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF] claims that solving PTP is equivalent to finding the optimal (in our case the longest) augmented path in the alignment graph G.

F = {v i1,1 , v i2,2 , . . . , v im,m } such that i k ≤ i k+1 , ∀k,

Mixed Integer Programming Formulation

Let y ik be binary variables associated with vertices in the previous network. Then y ik is 1 if block k is at position i and 0 otherwise (vertex (i, k) is activated or not). To take into account the interaction costs, we introduce a second set of variables 0 ≤ z ikjl ≤ 1, with (k, l) ∈ I and 1 ≤ i ≤ j ≤ n. The variable z ikjl is set to 1 if the corresponding edge is activated. Finding the optimal augmented path in graph G (i.e. solving PTP) is then equivalent to maximizing the following function :

m k=1 n i=1 C ik y ik + ((i,k),(j,l))∈E C ikjl z ikjl (1)

INRIA

This objective function is subject to the following constraints:

n i=1 y ik = 1 k ∈ [1, m] (2)
i j=1

y j(k+1) - i j=1 y jk ≤ 0 k ∈ [1, m -1], i ∈ [1, n -1] (3)
y ik ∈ {0, 1} k ∈ [1, m], i ∈ [1, n] (4)
y jl = j i=1 z ikjl (k, l) ∈ I, j ∈ [1, n] (5)
y ik = n j=i z ikjl (k, l) ∈ I, i ∈ [1, n] (6) 0 ≤ z ikjl ≤ 1 ((i, k), (j, l)) ∈ E (7)
This model, known as MYZ and first introduced in [START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF], has been shown to outperform the MIP model used in the RAPTOR package [START_REF] Xu | RAPTOR: optimal protein threading by linear programming[END_REF] for large PTP instances. Both models (MYZ and RAPTOR) are solved by first relaxing constraints (4), i.e., letting y ik be real variables such that 0 ≤ y ik ≤ 1. This allows the use of a Linear Programming (LP) technique. The solution of the LP technique is then used as a lower bound in a subsequent branch & bound algorithm that finds the integer solution.

Local alignments : towards better PTP models

Constraints (2) from the above model force each block (i.e SSE) to be aligned with the query sequence. However, proteins belonging to the same family do not always have the same number of SSEs. Forcing all blocks to be aligned with the sequence results in spurious alignments with bad scores. Such alignments prevent the method to detect remote homologs. To tackle this issue, we develop new models that allow a block to be omitted (not aligned) if its score with all other interacting blocks is negative. Such an approach realizes local sequencestructure alignments. Towards this goal, we slightly modify the definition of a feasible path by accepting that any non-decreasing set of vertices is a feasible path of the alignment graph G (i.e. the cardinality of a feasible path could be now less than m). Obviously, there is a one-to-one correspondence between the set of local alignments and the set of feasible paths in G.

In a local alignment, each block can be aligned with the whole sequence because the number of omitted blocks is not known in advance. It follows that relative positions are no longer meaningful. Thus, each block k takes n k = N -L k + 1 absolute positions along the sequence. This results in columns having different heights in the network flow formulation and leads to the need of using an "offset" to move from one column to the next (this offset is illustrated in Fig 5).

In order to implement the mechanism for omitting blocks during an alignment, we propose two schemes : (i) we modify constraints (2) to allow the sum of a column to be zero, (ii) N + 1 dummy vertices are added in every column.

If a dummy vertex is activated in a column, the corresponding block is omitted.

In the second scheme, two types of vertices exist : real (R) and dummy (D). This generates 4 types of subpaths: R to R , R to D, D to R, and D to D (as shown in fig 6). According to the above two schemes, we implemented 5 models : Compact model (CM) for the first one and Extended models (EM1, EM2, EM3, EM4) that make use of the dummy nodes. These models are described below.

INRIA

Mathematical models

In the following models, we use N , m, L k , I, C ik , C ikjl , y ik , z ikjl notations from section 2.3 and the following :

n k = N -L k + 1
is the number of possible positions of block k along the query sequence.

E = {((i, k), (j, l))|(k, l) ∈ I, 1 ≤ i ≤ n k , i + L k ≤ j ≤ n l } is the set of edges.
All models share the same objective function:

max m k=1 n k i=1 C ik y ik + ((i,k),(j,l))∈E C ikjl z ikjl (8)

Compact model (CM)

The objective function ( 8) is subject to the following constraints:

y ik ∈ {0, 1}, k ∈ M, i ∈ [1, n k ] (9) 0 ≤ z ikjl ≤ 1, ((i, k), (j, l)) ∈ E (10)
n k i=1 y ik ≤ 1, k ∈ M (11)
n l j=i+L k z ikjl -y ik ≤ 0, (k, l) ∈ I, i ∈ [1, n k ] (12)
j-L k i=1 z ikjl -y jl ≤ 0, (k, l) ∈ I, j ∈ [1, n l ] (13)
y ik + min(n l ,i+L k -1) j=1 y jl ≤ 1, 1 ≤ k ≤ l ≤ m, i ∈ [1, n k ] ( 14
)
n k i=1 y ik + n l i=1 y il - l j=L k +1 j-L k i=1 z ikjl ≤ 1, (k, l) ∈ I (15)
Constraints [START_REF] Marin | FROST: A Filter Based Fold Recognition Method[END_REF] correspond to constraints (2) in MYZ model and they impose a block to have one position or zero (i.e to be omitted). Constraints ( 12) and ( 13) force node activation if an outpout (resp. input) edge is activated. Constraints [START_REF] Lupyan | A new progressive-iterative algorithm for multiple structure alignment[END_REF] keep blocks order and impose that they do not overlap. Finaly, constraints (15) impose the activation of an edge if its ends are activated.

Extended model 1 (EM1)

Based on the compact model, each extended model is an attempt to find a more similar model with MYZ [START_REF] Andonov | Protein Threading Problem: From Mathematical Models to Parallel Implementations[END_REF] which was shown to be efficient for global alignments. EM1 extends the compact model by adding dummy positions for each column. An active dummy position in a column means that the corresponding block is deleted. d ik ∈ {0, 1} are dummy variables added for all extended models. The objective function ( 8) is subject to constraints ( 9), ( 10), ( 12), ( 13) and (15) from CM. We add the following constraints:

d ik ∈ {0, 1}, k ∈ M, i ∈ [1, N + 1] (16) n k i=1 y ik + N +1 i=1 d ik = 1, k ∈ M (17) j i=1 d ik + min(j,n k ) i=1 y ik - j i=1 d ik-1 - j-L k-1 i=1 y ik-1 ≤ 0, k ∈ [2, m], j ∈ N + 1 (18)
Constraints (17) replace constraints [START_REF] Marin | FROST: A Filter Based Fold Recognition Method[END_REF] with the use of dummy nodes. Our goal was to transform these inequalities in equalities. The use of dummy nodes also change constraints ( 14) into constraints (18). One can notice that, m -1, the number of constraints ( 18) is smaller than, m * (m-1)

2

, the number of constraints (14).

Extended model 2 (EM2)

This model extends EM 1 by adding edges between dummy vertices and original vertices (or real vertices). The aim of this model is to transform constraints ( 12) and ( 13) in equalities.

rd ikjl ∈ R, (k, l) ∈ I, (i, j) ∈ RD are dummy variables representing an edge from a dummy vertex to a real vertex(with

RD = {(i, j) |(k, l) ∈ I 1 ≤ i ≤ n k , i + L k ≤ j ≤ N + 1}).
dr ikjl ∈ R, (k, l) ∈ I, (i, j) ∈ DR are dummy variables representing an edge from a real vertex to a dummy vertex(with

DR = {(i, j) |(k, l) ∈ I, 1 ≤ i ≤ j ≤ n l }).
The objective function ( 8) is subject to constraints ( 9), ( 10), ( 15), ( 16), ( 17), (18) and to the following constraints:

0 ≤ rr ikjl ≤ 1, (k, l) ∈ I, (i, j) ∈ RD (19) 0 ≤ dr ikjl ≤ 1, (k, l) ∈ I, (i, j) ∈ DR (20) n l j=i+L k z ikjl + N +1 j=i+L k rd ikjl -y ik = 0, (k, l) ∈ I, i ∈ [1, n k ] (21) j-L k i=1 z ikjl + j i=1 dr ikjl -y jl = 0, (k, l) ∈ I, j ∈ [1, n l ] (22)
Constraints ( 12) and ( 13) are replaced by constraints (21) and ( 22) which are equalities.

INRIA

Extended model 3 (EM3)

This model extends EM 2 by adding edges between dummy vertices.

dd ikjl ∈ R, (k, l) ∈ I, (i, j) ∈ DD are dummy variables representing an edge from a dummy vertex to a dummy vertex(with

DD = {(i, j) | 1 ≤ i ≤ j ≤ N + 1}).
The goal of this model is to have equality constraints like (21), ( 22) applied to dummy vertices. The objective function ( 8) is subject to constraints ( 9), ( 10), ( 16), ( 18), ( 19), ( 20), ( 21), ( 22) and to the following:

0 ≤ dd ikjl ≤ 1, (k, l) ∈ I, (i, j) ∈ DD (23) N +1 j=i dd ijkl + n l j=i dr ijkl -d ik = 0, (k, l) ∈ I, i ∈ [1, N + 1] ( 24
) j i=1 dd ijkl + j-L k i=1 rd ijkl -d jl = 0, (k, l) ∈ I, j ∈ [1, N + 1] (25)
Constraints (24) (resp (25)) impose the activation of one and only one edge going out (resp. in) a dummy node. These constraints replace activation constraints (15).

Extended model 4 (EM4)

This model is based on EM 2. In this model, we try to delete constraints (15). The objective function ( 8) is subject to constraints ( 9), ( 10), ( 16), ( 17), ( 18), ( 19) , ( 20) , ( 21), ( 22) and to the following :

n l j=i dr ijkl -d ik ≤ 0, (k, l) ∈ I, i ∈ [1, N + 1] (26) j-L k i=1 rd ijkl -d jl ≤ 0, (k, l) ∈ I, j ∈ [1, N + 1] (27)
In this model we replace constraints (15) by constraints ( 26) and ( 27).

Results

MIP models comparison

Benchmark

These five MIP models have been implemented with Ilog CPLEX 10.0 Library in our protein threading package FROST [START_REF] Marin | FROST: A Filter Based Fold Recognition Method[END_REF]. We also randomly created a benchmark of 100 alignments to compare their performances. Each alignment can be solved in less than 1000 seconds with CPLEX MIP solver for all models. Alignments have been processed on a cluster of computers with 2.5GHz AMD Opteron biprocessors and 4MB of memory. Comparisons of models have been carried out using four measures : the number of variables, the number of constraints, the relative gap, and the computation times.

Number of variables and constraints

CM obviously has the smallest number of variables because there is only real nodes and real edges. EM1 is larger than CM because dummy nodes are added. Then EM2 and EM3 follow due to the addition of edges. EM2 and EM4 have the same number of variables. In summary, based on the number of variables, the ranking of the models is: CM < EM1 < EM2 = EM4 < EM3.

EM1 is the model with the smallest number of constraints. Although EM1 has more variables, dummy nodes create less constraints (18) compared to constraints [START_REF] Lupyan | A new progressive-iterative algorithm for multiple structure alignment[END_REF] in CM for the same effect. EM1 and EM2 have similar structure, but EM1 has also less constraints than EM2. Finally, EM3 and EM4 have the largest number of constraints, about twice the number of constraints in EM1.

Relative Gap

The relative gap (RG) is the relative difference between the solution of the relaxed problem (LP ) and the optimal solution (OP T ):

RG = LP -OP T OP T
. RG is a good indicator of the efficiency of the model since the smaller RG, the easier for the branch & bound algorithm to find the solution. EM3 relative gap values are on the thick line. One can notice that EM3 and EM4 are the tightest models, and that EM1 always give an equal relative gaps with EM2. Overall, the tightest model is EM3 because it gives a tighter (or equal) relative gap than EM4 for all instances.

Figure 7 shows that EM3 and EM4 give tighter relative gaps compared to other models. In fact, EM3 always gives a tighter relative gap than any other model. Moreover, EM3 finds the optimal solution by solving the relaxed problem for 77% of alignments. This rate if greater than for CM (37%), EM1 and EM2(64%), and EM4(75%).

Computation Times

Figure 8 presents statistics on computation times for the five models. EM1 and EM2 seem to be the fastest models compared to CM, EM3 and EM4.

Actually, EM1 computation times are always smaller than the others.

Biological considerations

In order to measure the quality of local alignments, we need accurate alignments between proteins. These alignments are provided by structural alignments of proteins. We used the TopMatch server [START_REF] Sippl | A note on difficult structure alignment problems[END_REF][START_REF] Sippl | On distance and similarity in fold space[END_REF] for alignments of two proteins and the Mammoth server [START_REF] Lupyan | A new progressive-iterative algorithm for multiple structure alignment[END_REF] for multiple structure alignments.

Local similarities

A first improvement of a local alignment approach compared to a global alignment approach is to align a long template with a small sequence. Such alignments are not allowed by global alignments because all blocks must be aligned. For example, the TopMatch alignment of proteins 9gaaA (length 180) and 1apyB (length 141) show that 1apyB is a sub domain of 9gaaA (85% structure similarity and 40% sequence identity). Thus, a global alignment between template 9gaaA and sequence 1apyB is impossible. However, with a local alignment, a part of 9gaaA can be deleted and the subdomain can be aligned with 1apyB as illustrated in figure [START_REF] Andonov | Recent Advances in Solving the Protein Threading Problem[END_REF].

Better alignments

In section 3 we presented the problem of having proteins with similar structures but not the same number of SSEs. Such a case is illustrated in figure (10.A). Template 1qddA is aligned with sequence 1rjhA. These two proteins share a common domain (76% structure similarity and 21% sequence identity) but the first three blocks of 1qddA (in gray) are not in this domain. Because the global alignment approach imposes these three blocks to be aligned with the sequence, it results in a decrease of the overall score and in a spurious alignment. GNLV---TQAEGAFVASLIKESG----NVW--LHD-------------------WHW------------PQARIS----TNA----CYYFN----TWVDADLYCQNM---GNLV----TQAEGAFVASLIKESG-----NVW---LHD----WHW-PQARIS----TNA----CYYFN----TWVDADLYCQNM---GNLV----TQAEGAFVASLIKESG-----NVW---LHD----WHW-PQARIS----TNA----CYYFN----TWVDADLYCQNM---GNLV----TQAEGAFVASLIKESG-----NVW---LHD----WHW-PQARIS----TNA----CYYFN----TWVDADLYCQNM---GNLV----TQAEGAFVASLIKESG-----NVW---LHD----WHW [START_REF] Lupyan | A new progressive-iterative algorithm for multiple structure alignment[END_REF]. Frost alignment is in bold. A) Gray blocks are aligned although they are not in the structural domain shared by 1qddA and 1rjhA. This results in an alignment which is totaly different from the best alignment. B) Because two gray blocks have been omitted, local alignment approach gives a better score and an alignment more coherent with the best alignment (frames correspond to exact matchings between Frost and Mammoth alignments).

A local alignment permits to omit blocks when needed (i.e. based on score function) as illustrated in figure (10.B). The deletion of two gray blocks results in a better score and an alignment which is more consistent with the best alignment found by Mammoth.

Discussion

This paper deals with the description and the comparison of five MIP models for local alignments. We do not explore the sensitivity and specitivity of the approach from a biological point of view. Such statistics about recognition rates and quality of alignments need a huge amount of alignments and are beyond the scope of this paper.

Based on the compact model, we develop four extended models. The aim of these extended models is to find a similar model to MYZ. MYZ has been proven to be very efficient and to give an optimal solution for 95% of cases by solving the relaxed problem. Our assumption is that a close model to it would be as efficient for local alignments as MYZ is for global alignments. EM3 is the closest model to MYZ because it uses very similar equality constraints. By always giving the tighter relative gap, EM3 confirms our hypothesis.

Aligning difficult instances (i.e. long sequences with long templates) with the CPLEX solver is too slow for a massive usage. For this reason, we need INRIA to develop a fast dedicated algorithm based on one of our models. Usually, we use the model with the tightest relative gap but in this case, the tightest model (EM3) is also the slowest (section 4.1.4) and the largest (section 4.1.2). Thus, choosing the best model for the development of a dedicated solver is still an open research problem on which we are currently working.

Although the local alignment approach seems to behave as we expected, it completely relies on the score function. Indeed, a block is omitted if its score decreases the overall alignment score, i.e. its score is negative. But we do not know if 0 is the right boundary to decide to omit a block. We are currently looking for an automatic learning method to evaluate the correct boundary for a given score function.

Conclusion

This paper describes a method allowing local alignments of protein sequences and structures. As far as we know, this is the first attempt to propose alignments that give the possibility of omitting parts of the 3D structure that might not be conserved in remote homolog proteins. Results show that our algorithm is indeed capable of omitting blocks when required. The accuracy of the resulting alignment strongly depends on the quality of the score function. It is likely that score functions adapted to the local alignments will contribute for the success of the approach.

Local alignments have been modelled and tested with the CPLEX 10 solver. This is a general purpose solver for integer programming problems that is very convenient for rapidly testing new ideas. However, our computation experiments show that it is not enough rapid for our particular problem. We show that EM1 and EM3 are respectively the fastest and the tightest model. The next step is to choose one of these models to create a dedicated algorithm in order to solve difficult instances and to evaluate the recognition rate of the local alignment approach.

Figure 2 :

 2 Figure 2: Proteins as contact maps. A) "ribbon" representation of a protein with 3 α helices, B) only the secondary structure elements, here the 3 helices, are taken into consideration, C) dotted lines represent interactions between positions in the helices, i.e., residues in contact, D) contact map representation of the protein, double headed arrows correspond to the dotted line in C. In the blocks respectively 6, 4 and 8 amino acids can be aligned.

Figure 4 :

 4 Figure 4: Example of alignment graph. The thick line path corresponds to the threading in which the relative positions of blocks are 1,2,2,3,4,4. Dashed line edges represent non local interactions where the set of block interactions is I = {(1, 3), (2, 5), (3, 5)}.

Figure 5 :Figure 6 :

 56 Figure 5: Example of flexible alignment graph. The path in thick lines corresponds to the threading in which: blocks 1, 2, and 4 are omitted, and blocks 3, 5, and 6 are on positions 1, 3 and 5. All edges in E are represented. Because of absolute positions, notice that a offset of size L k is needed to go from column k to the followings (blocks sizes are : 2,3,2,3,2,2).

Figure 7 :

 7 Figure7: Relative gap of CM, EM1, EM2 and EM4 compared to relative gap of EM3 for 100 alignments. EM3 relative gap values are on the thick line. One can notice that EM3 and EM4 are the tightest models, and that EM1 always give an equal relative gaps with EM2. Overall, the tightest model is EM3 because it gives a tighter (or equal) relative gap than EM4 for all instances.

Figure 8 :

 8 Figure 8: Statistics on computation times of each model for 100 alignments. Times are in seconds on a logarithmic scale. A box represents 50% of instances, the thick line inside is the median. The ends of lines represent maximum and minimum values. EM1 and EM2 are the fastest models.

Figure 9 :

 9 Figure 9: Alignment of proteins 9gaaA and 1apyB by TopMatch and FROST. Although 9gaaA is larger than 1apyB, a local alignment is possible. Moreover, FROST local approach produces a coherent alignment with the Top-Match structural alignment (90% of exact matching).

FTQTKTFHEASEDCISRGGTLSTPQTGSENDALYEYLRQSVGNEAEIWLGLNDMAAEGTWVDMTGARIAYKNWETEITAQPDGGKTE

Figure 10 :

 10 Figure10: Frost Global and Local Alignments of core 1qddA with query 1rjhA. "Best Align" is the alignment obtained by Mammoth structural alignment tool[START_REF] Lupyan | A new progressive-iterative algorithm for multiple structure alignment[END_REF]. Frost alignment is in bold. A) Gray blocks are aligned although they are not in the structural domain shared by 1qddA and 1rjhA. This results in an alignment which is totaly different from the best alignment. B) Because two gray blocks have been omitted, local alignment approach gives a better score and an alignment more coherent with the best alignment (frames correspond to exact matchings between Frost and Mammoth alignments).

RR n 7122

Collet & al.

Acknowledgements

This work is supported by ANR project PROTEUS "ANR-06-CIS6-008". N. Yanev is supported by DVU/01/197 project from SWU "Neofit Rilski" and by DO 02359/2008 grant. All computations were done on the Biogenouest bioinformatics platform 1 .