Jérôme Malick
email: jerome.malick@inria.fr

Frédéric Roupin
email: frederic.roupin@lipn.univ-paris13.fr

Solving k-cluster problems to optimality with semidefinite programming

Keywords: combinatorial optimization, k-cluster, exact resolution, semidefinite programming, Lagrangian relaxation, branch-and-bound

This paper deals with the computation of exact solutions of a classical NP-hard problem in combinatorial optimization, the k-cluster problem. This problem consists in finding a heaviest subgraph with k nodes in an edge weighted graph. We present a branch-and-bound algorithm that applies a novel bounding procedure, based on recent semidefinite programming techniques. We use new semidefinite bounds that are less tight than the standard semidefinite bounds, but cheaper to get. The experiments show that this approach is competitive with the best existing ones.

and split graphs [START_REF] Corneil | Clustering and domination in perfect graphs[END_REF]. Strong negative results have been obtained about the approximability of the clique problem, but there is no direct way to use them for the k-cluster problem: it is open whether an approximation algorithm with a fixed ratio exists in general. The best general approximation ratio is actually O(n 1 3 ) [START_REF] Feige | The dense k-subgraph problem[END_REF], and several approximation algorithms with a better ratio have been designed for special graphs (see e.g. [START_REF] Liazi | A constant approximation algorithm for the densest k-subgraph problem on chordal graphs[END_REF] for recent results and references).

These results give an idea of how tough the k-cluster problem is to be solved to optimality. Even if the machinery of linear-optimization-based branch-and-bound algorithms has reached a high degree of development (see for example the recent [START_REF] Achterberg | Branching rules revisited[END_REF]), only medium-size instances of the k-cluster problem can be solved in general with these solvers (see [START_REF] Erkut | The discrete p-dispersion problem[END_REF], [START_REF] Billionnet | Different formulations for solving the heaviest k-subgraph problem[END_REF], [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF]). Linear relaxations seem to be not tight enough for k-cluster in order to be used efficiently within an algorithm for exact resolution. Therefore, the best exact and approximate resolution methods devoted to this problem are based on nonlinear approaches, such as convex quadratic programming [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] or semidefinite programming [START_REF] Han | An improved rounded method and semidefinite programming relaxation for graph partition[END_REF], [START_REF] Roupin | Form linear to semidefinite programming: an algorithm to obtain semidefinite relaxations for bivalent quadratic problems[END_REF], [START_REF] Jäger | Improved approximation algorithms for maximum graph partitioning problems[END_REF].

We present in this paper an approach for solving k-cluster problems to optimality by using standard branch-and-bound techniques together with new semidefinite programming bounds. The bounding procedure trades computing time for a (small) deterioration of the quality of the bound, so that it fits well within a branch-and-bound algorithm. To get these semidefinite bounds, we follow the scheme sketched in [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF] for general binary quadratic program: we simplify and specify it for the k-cluster problem -and we push the development to show numerically that this approach is efficient. The main contribution of this paper is to prove numerically that these new semidefinite bounds are interesting in view of solving hard combinatorial optimization problems to optimality. For large-scale k-cluster problems indeed, our approach is competitive with the best known approach [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF], [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] that uses CPLEX.

The paper is organized as follows. We formulate the k-cluster problem in Section 2 in a way to introduce the new semidefinite bounds in Section 3. We compare these bounds with the usual semidefinite bound in Section 4 to argue that they have nice features in view of exact resolution. Finally we embed them within a simple branchand-bound algorithm in Section 5 and compare the overall efficiency of this approach with the best known solver for exact resolution of k-cluster problems.

Equivalent formulations

This section introduces the notation used in the sequel, and presents a sequence of equivalent formulations of the k-cluster problem: namely, the initial formulation (1), simplified (2), augmented with additional constraints (3), the formulation as a binary quadratic programming problem (4) and the formulation as semidefinite programming problems (6), (8). All these formulations are standard, except the last one. We emphasize: all of them are exact reformulations; we start relaxing in the next section.

Initial modeling as a {0, 1}-quadratic programming problem. Consider an undirected weighted graph G = (V, E) with n vertices {v 1 , . . . , vn} and with nonnegative weights w ij on edges (v i , v j ). For an integer k in {2, . . . , n -2}, the so-called k-cluster problem consists in determining a subset S of k vertices such that the total edge weight of the subgraph induced by S is maximized. To select subgraphs, assign a decision variable z i ∈ {0, 1} for each node (z i = 1 if the node is taken, and z i = 0 if the node is not). The weight of the subgraph given by z is P (i,j)∈E w ij z i z j . Thus the k-cluster problem can be phrased as the 0-1 quadratic problem

 max 1 2 P n i=1 P n j=1 w ij z i z j P n i=1 z i = k, z ∈ {0, 1} n .
(1)

We rewrite this problem synthetically as

( max 1 2 z ⊤ W z e ⊤ z = k, z ∈ {0, 1} n .
(2)

with the vector of all ones e = (1, . . . , 1) ∈ R n and the weight-matrix W = (w ij ) ij of the graph G (which is a symmetric n × n matrix).

Reformulation 1: Enforcement of redundant constraints. In view of relaxing, we add redundant constraints to have a better duality gap. We introduce the standard productconstraints as follows. Observe that a feasible z (that is, such that e ⊤ z = k and z ∈ {0, 1} n ) also satisfies for all j ∈ {1, . . . , n}

n X i=1 z i z j = kz j .
The left-hand side of this equality is quadratic in z; so we introduce the symmetric n × n-matrix C j such that P n i=1 z i z j = 1 2 z ⊤ C j z. Adding these n product constraints, we come up with the following equivalent formulation of k-cluster

8 > < > : max 1 2 z ⊤ W z e ⊤ z = k, z ∈ {0, 1} n z ⊤ C j z = 2k z j , j ∈ {1, . . . , n}.
(3) These product constraints are the only one we add, because they are tight in a certain way; this follows indeed from the general results that we recall briefly here. For a {0, 1}-quadratic problem, it is well-known (see e.g. [START_REF] Lemaréchal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]) that the best bound by Lagrangian duality is obtained when dualizing only the {0, 1}-constraints (and not the linear constraints). Though this bound does not lead directly to a SDP problem, the semidefinite relaxation of the general problem reinforced by the product constraints (e.g. (3) for k-cluster) is equivalent to this best bound [START_REF] Faye | Partial lagrangian for general quadratic programming. 4'OR[END_REF]. Adding other redundant quadratic equality constraints would lead to the same bound, so we stick with (3).

We also recall that adding to (2) the single constraint (e ⊤ z -k) 2 = 0 instead of the n product constraints is an equivalent approach: the two formulations lead to two SDP relaxations which give the same bound -but on which solvers behave differently. The numerical comparison between the two SDP formulations is made in [START_REF] Malick | Numerical study of SDP bounds for the k-cluster problem[END_REF] showing that (3) is preferable for our developments. We will come back to this in Section 4.2.

Reformulation 2: Change of variables. Dealing with quadratic problems with {-1, 1} constraints and purely quadratic problems will ease the next developments. So we introduce the new variable x = (x 0 , . . . , xn) ∈ R n+1 performing two operations:

1. The change of binary variable x i = 2z i -1 ∈ {-1, 1} for i ∈ {1, . . . , n}.

The "homogenization" of the quadratic forms with the extra-variable x

0 ∈ {-1, 1}.
See more precisely the role of this homogenization in Lemma 3 in appendix. A graph interpretation is as follows: we add an isolated vertex x 0 to the graph, and we want to compute a (k + 1)-cluster containing this vertex; in other words, x 0 is like a "flag" that indicates which vertices are in the k-cluster.

The change of variable z → x gives the following equivalent formulation of k-cluster

8 > < > : max x ⊤ Q x x ⊤ Q j x = 4k -2n, j ∈ {0, . . . , n} x ∈ {-1, 1} n+1 (4)
where the n + 2 symmetric matrices Q and Q j (for j ∈ {0, . . . , n}) are defined by

Q := 1 4 » e ⊤ W e e ⊤ W W e W - , Q 0 := » 0 e ⊤ e 0 - and Q j := » 0 ẽ⊤ j ẽj C j - (5)
with ẽj = e + (n -2k)e j the vector of R n made up from e and e j the j-th element of the canonical basis of R n . The change of variable and the reformulation are detailed in Lemma 3 in appendix.

Reformulation 3: Lifting in matrix space. We now lift the problem (4) up to the matrix space S n+1 in order to transform the quadratic forms with respect to x ∈ R n to linear forms with respect to a matrix variable X. We follow the classical pattern (e.g. [START_REF] Lovász | On the Shannon capacity of a graph[END_REF], [START_REF] Goemans | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF]) and we use classical notation, recalled below. The natural inner product in matrix space S n+1 is defined for any X, Y ∈ S n+1 by

X, Y = n+1 X i,j=1 X ij Y ij = trace(XY ).
This inner product is very convenient for our purposes through the relation:

∀ x ∈ R n+1 , ∀ A ∈ S n+1 , x ⊤ A x = A, xx ⊤ .
We denote • the associated norm; it is the same notation as the norm in R n , but no confusion should be possible since the matrices are represented by capital letters. The idea of the standard lifting is to introduce the symmetric matrix of rank one X = xx ⊤ . With the help of X, we express the binary constraints x i ∈ {-1, 1} (that is x i 2 = 1) as X ii = 1, and the quadratic constraint x ⊤ A x = c as A, X = c. So we get the following equivalent formulation of (4)

8 > > < > > : max Q, X Q j , X = 4k -2n, j ∈ {0, . . . , n} X ii = 1, i ∈ {0, . . . , n} X = xx ⊤ .
Notice then that X = xx ⊤ is a rank-one symmetric semidefinite matrix, and that conversely any rank-one symmetric semidefinite matrix can be written this way. So we can cast the above problem as an SDP linear problem with rank-one constraint

8 > > < > > : max Q, X Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n} rank(X) = 1, X 0. ( 6
)
where E i is the matrix with zeros entries except in position (i, i) where there is a one.

Reformulation 4: Introducing the spherical constraint. The hard constraint in (6) is the rank-one constraint, which is moreover difficult to handle. So we go one step further to the standard SDP lifting by rewriting this rank-one constraint as a norm, along the lines of [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF]. The (curious and) useful property is the following (see [Mal07, Theorem 1]): for all X 0 satisfying X ii = 1, we have X ≤ n + 1 and moreover

X = n + 1 ⇐⇒ rank X = 1. ( 7
)
Therefore the idea is to replace the rank-one constraint in the SDP formulation of the k-cluster by the constraint X 2 = (n + 1) 2 , called the "spherical constraint". Thus we have the following formulation of k-cluster that we formalize in the next proposition; this is the formulation we will use in this paper.

Lemma 1 With the notation of this section, the optimal value of the k-cluster problem is equal to the optimal value of

8 > > > > < > > > > : max Q, X Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n} X 2 = (n + 1) 2 X 0. ( 8
)
Proof The proof is obvious in view of the previous developments: by Lemma 3, the optimal value of (1) is equal to the one of (4) which is in turn equivalent to (6), and then to (8) by (7).

⊓ ⊔

We have transformed a quadratic problem with quadratic constraint in R n+1 to a linear SDP problem with one single quadratic constraint (the spherical constraint X 2 = (n + 1) 2 ). This quadratic constraint carries the nonconvexity and the combinatorial difficulty of the initial problem. The idea is now to treat it by Lagrangian dualization, as explained in the Section 3.2.

Semidefinite relaxations

Standard semidefinite bounds

The previous section explains the formulation of the k-cluster problem as SDP problems (6) and (8). All the combinatorial difficulty of the original problem is now located in the rank-one constraint in (6), or in the spherical constraint in (8). The standard SDP relaxation then consists in enlarging the feasible set in (6) by dropping the rank constraint, to derive a convex problem. In our situation, the standard SDP relaxation of the k-cluster problem (6) is indeed

8 > > < > > : max Q, X Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n} X 0. (9)
We note that this SDP problem is equivalent to the one obtained from the {0, 1}formulation (3) (used by [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] in particular). As mentioned in the discussion following (3), it is also equivalent to the SDP problem coming from (2) augmented by the squared constraint (e ⊤ x -k) 2 = 0. More precisely, by following the same lines to go from (3) to (9), we get the equivalent SDP problem

8 > > < > > : max Q, X Q S , X = -(2k -n) 2 E i , X = 1, i ∈ {0, . . . , n} X 0 (10)
with a particular Q S ∈ S n+1 . In practice, this bound is known to be tight and several approximation algorithms have been proposed using it (see e.g. [START_REF] Han | An improved rounded method and semidefinite programming relaxation for graph partition[END_REF]).

Those SDP problems (9) and ( 10) have linear objective functions and as well as linear constraints together with the conic constraint. Over the last decade, several algorithms and associated solvers have been developed to solve problems of this type. We refer for instance to the introduction of the recent [START_REF] Malick | Regularization methods for semidefinite programming[END_REF] for references and for some explanations about the different approaches; here we just point out:

-Interior points. The most known and used methods for SDP; they have nice complexity theory (see e.g. the review [Tod01]), and reliable numerical performances.

For the forthcoming numerical experiments of Section 4.2, we use the software CSDP [START_REF] Helmberg | An interior point method for semidefinite programming[END_REF], [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF]. -Spectral bundle. The bundle methods of convex optimization [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms[END_REF] have been specialized successfully to semidefinite programming, and especially to applications of semidefinite programming in combinatorial optimization [START_REF] Helmberg | A spectral bundle method for semidefinite programming[END_REF]. For the numerical experiments, we use the software SB [START_REF] Helmberg | A C++ implementation of the Spectral Bundle Method[END_REF]. -Regularization. Recently other approaches based on regularization (augmented Lagrangian and proximal method) have been developed for solving large-scale SDP problems, as for example the relaxation of max-stable problem [START_REF] Malick | Regularization methods for semidefinite programming[END_REF]. Let us also mention the penalized augmented Lagrangian method of PENNON [START_REF] Kočvara | On the solution of large-scale sdp problems by the modified barrier method using iterative solvers[END_REF]. -First-order methods. To tackle ever larger SDP problems, several first-order methods have been developed: we mention for example low-rank methods (see [START_REF] Burer | A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization[END_REF] and references therein) that have recently drawn renewed interest.

Though there exist several types of methods and several efficient solvers, computing the SDP bound (9) or (10) is still expensive, and this may prevent its direct use as bounding procedure within a branch-and-bound for solving k-cluster to optimality. Note in particular that the method of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] uses the SDP bound only once at the beginning of the search tree (as an initial calibration of the quadratic convex relaxation used later for bounding, see more in Section 5). We are not aware of research for kcluster in the line of [START_REF] Rendl | Solving max-cut to optimality by intersecting semidefinite and polyedral relaxations[END_REF] for max-cut. Here, we consider different bounds that have an SDP-quality but that are less tight that (9) and (10). In the next sections, we argue that those new bounds are easier to compute and then well-adapted to be embedded within a branch-and-bound.

New family of semidefinite bounds

We investigate now a way to keep a SDP-like quality of bound without paying the full computational price to get it. The idea is to keep the rank-one constraint, to write it with the help of the spherical constraint as in (8) and then to dualize it. This approach is sketched in [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF] in a general setting; here, we specialize the study for the k-cluster to push to the end and to get tools for exact resolution.

For the real parameter α ∈ R, we consider the Lagrangian function

L(X; α) := Q, X - α 2 ( X 2 -(n + 1) 2 )
and the associated dual function

Θ(α) := 8 > > < > > : max L(X; α) Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n} X 0. ( 11
)
By weak duality, each value Θ(α) is an upper bound for the optimal value of k-cluster, since we can write: for all feasible X for (8) (and then in (11)),

Q, X = L(X, α) ≤ Θ(α) (12)
so that Θ(α) is upper bound for (8) indeed.

Hence we have a new family of SDP bounds Θ(α) (parameterized by α ∈ R). In a way, these bounds generalize the standard SDP bound (9): observe indeed that for α = 0, (11) is exactly (9). They have moreover interesting properties: the important theoretical properties are gathered in the following proposition; we discuss the numerical properties in the next section. In what follows, val( * ) denotes the optimal value of an optimization problem numbered by ( * ).

Proposition 1 Function Θ : R → R defined by (11) satisfies the following properties:

-For any α ∈ R, we have a bound for k-cluster Θ(α) ≥ val (8) = val (1).

-If α > 0, this bound is weaker than the SDP bound -If α > 0, we may compute Θ(α) by solving a semidefinite least-squares problem:

Θ(α) ≥ val (9) ≥ val (8) = val (1),
Θ(α) = " α 2 (n + 1) 2 + 1 2α Q 2 " -α 8 > > < > > : min 1 2 X -Q/α 2 Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n} X 0 (13)
Proof The first property comes from (12). The second property about the comparison between Θ(α) and Θ(0) = val(9) follows from Theorem 3 of [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF]: the function Θ : R → R is convex (so continuous at 0) and nondecreasing (so Θ(α) ≥ Θ(0) if α > 0). To prove the third point, we write for α > 0

L(X, α) = α 2 (n + 1) 2 - α 2 ( X 2 -2 X, Q/α ) = α 2 (n + 1) 2 - α 2 " w w w w X - Q α w w w w 2 - Q 2 α 2 " = " α 2 (n + 1) 2 + 1 2α Q 2 " - α 2 w w w w X - Q α w w w w 2 ,
which yields that (11) corresponds to (13).

⊓ ⊔

The point is thus the following: for α > 0, the new SDP bound Θ(α) is always less tight than the usual SDP bound Θ(0), but it can be reduced arbitrary close to it by decreasing α. The key question is now how easier is Θ(α) to compute; this is addressed in the next section.

Relaxed resolution: computation of bounds

The new SDP bound Θ(α) looks more complicated than the usual SDP bound Θ(0). It turns out that nice geometrical properties make it cheaper to compute. Section 4.1 explains the computation, and then Section 4.2 presents a numerical comparison.

Computing the bound by projection

By (13), computing the new SDP bound Θ(α) when α > 0 comes down to solving 8 < :

min 1 2 X -Q/α 2 Q j , X = 4k -2n, j ∈ {0, . . . , n} E i , X = 1, i ∈ {0, . . . , n}, X 0. ( 14
)
This quadratic SDP problem has a simple geometric interpretation: Consider indeed the affine subspace Aff in S n+1 defined by the affine constraints of this problem

Aff := n X ∈ S n+1 : Q j , X = 4k -2n, E i , X = 1, for i, j ∈ {0, . . . , n} o .
Then the problem ( 14) is

 min X -Q/α 2 X ∈ Aff ∩ S + n+1
and thus consists in projecting the matrix Q/α onto the intersection Aff ∩ S + n+1 of the affine subspace Aff and the cone S + n+1 . This problem is an instance of the so-called semidefinite least-squares [START_REF] Malick | A dual approach to semidefinite least-squares problems[END_REF]. Efficient algorithms have been recently developed to solve these problems, based on three paradigms: (1) alternating projection methods [START_REF] Higham | Computing a nearest symmetric correlation matrix -a problem from finance[END_REF], (2) interior-point methods [START_REF] Tutuncu | Inexact primal-dual path-following algorithms for a special class of convex quadratic SDP and related problems[END_REF], (3) duality [START_REF] Malick | A dual approach to semidefinite least-squares problems[END_REF] (see also [START_REF] Qi | Quadratic convergence and numerical experiments of Newton's method for computing the nearest correlation matrix[END_REF] and [START_REF] Borsdorf | A preconditionned newton algorithm for the nearest correlation matrix[END_REF] for developments on an important special case).

We have developed our own semidefinite least-squares solver to solve (14). This solver follows the general method of [START_REF] Malick | A dual approach to semidefinite least-squares problems[END_REF] while being specific to be more efficient in exploiting the particular features of (14). Interesting properties are summarized in the next proposition. theorem. Consider the linear operator A : S n+1 → R 2n+2 defined by A(X) := ` Q 0 , X , . . . , Qn, X , E 0 , X , . . . , En, X ´, and the vector b ∈ R 2n+2 defined by b := (4k -2n, . . . , 4k -2n, 1, . . . , 1), such that the affine subspace Aff admits AX = b as an equation. Notice that the adjoint of A, denoted by A * : R 2n+2 → S n+1 , and defined by A * (y), X = y ⊤ A(X) for all y ∈ R 2n+2 and X ∈ S n+1 , is

A * (y) = n X i=0 y i Q i + n X i=0 y i+n+1 E i .
For a given y ∈ R 2n+2 , we also define the matrix X(y) := P S + n+1 `Q/α + A * (y) ẃhere

P S + n+1
is the projection onto S + n+1 . Note that the dependence with respect to α is dropped for simplicity.

Proposition 2 The function f : R 2n+2 → R defined by

f (y) := 1 2 X(y) 2 -y ⊤ b
is convex and differentiable with gradient ∇f (y) = AX(y) -b. Besides, we have

val (14) ≥ Q 2 2α 2 -f (y) ( 15
)
for all y ∈ R 2n+2 . Assume furthermore that there exists a vector y ⋆ ∈ R 2n+2 that attains the minimum f ⋆ of f ; then equality holds in (15) and X ⋆ = X(y ⋆ ) is the unique solution of (14).

Proof Up to a change of sign and to the constant

Q 2
2α 2 , the function f corresponds to the function θ in [START_REF] Malick | A dual approach to semidefinite least-squares problems[END_REF]. We just apply results of this paper: the convexity and differentiability of f comes from Theorems 3.1 and 3.2; (15) is equation (4.2); and the rest comes from Theorem 4.1.

⊓ ⊔

The differentiability of f allows us to use standard algorithms for unconstrained differentiable optimization for solving (14) thus to compute the bound Θ(α). We can cite, among others: gradient, quasi-Newton, Newton-like methods as truncated generalized Newton (see [START_REF] Nocedal | Numerical Optimization[END_REF], [START_REF] Bonnans | Numerical Optimization[END_REF]). For its simplicity and robustness, we choose the limited memory quasi-Newton algorithm [START_REF] Gilbert | Some numerical experiments with variablestorage quasi-Newton algorithms[END_REF], which has proved to be efficient in many (academic and industrial) applications.

Numerical comparisons of semidefinite bounds

This section gives numerical comparisons of the new SDP bound Θ(α) and the usual SDP bound Θ(0) in terms of tightness, computation time and balance between both.

Solvers and settings. To compute the SDP bound Θ(0), that is solving (9), we use:

-SB [Hel04], a SDP solver based on the spectral bundle method [START_REF] Helmberg | A spectral bundle method for semidefinite programming[END_REF]. This software can handle large-scale problems and is known to be efficient in the context of combinatorial optimization (see e.g. the recent [START_REF] Rendl | Solving max-cut to optimality by intersecting semidefinite and polyedral relaxations[END_REF]). We use it with default settings, except that we add an initial scaling of the constraints and we set the stopping criterion to 1e-4 (instead of 1e-5).

-CSDP [START_REF] Borchers | CSDP, a C library for semidefinite programming[END_REF], a robust and efficient interior point solver. We use it with default settings, except that we activate "Fastmode" and we set the stopping criterion "objtol" to 1e-5 (instead of 1e-8).

To compute the SDP bound Θ(α), that is essentially solving (14), we use:

-SDLS, a home-made solver implementing a version of the algorithm of [START_REF] Malick | A dual approach to semidefinite least-squares problems[END_REF] adapted for solving (14). We use the quasi-Newton algorithm [START_REF] Gilbert | Some numerical experiments with variablestorage quasi-Newton algorithms[END_REF] to minimize the function f of Proposition 2. We set the stopping criterion to 1e-7, and, as for SB, the constraints of ( 14) are scaled so that the constraints matrices are of norm 1. For sake of simplicity, we also decide to keep α constant. In view of preliminary experiments, we then fix α to 1e-4: our strategy is indeed to have a SDP-like bound, and Θ(10 -4 ) is almost as tight as Θ(0), while its computing time is reasonable. We illustrate this with Figure 1 which plots (for an instance with n = 300) the decrease of Θ(α) to Θ(0) when α → 0 and the increase of the corresponding computing time. For example, Θ(α) for α = 10 -4 (the leftmost point) is much more tight than for α = 0.002 with less than twice the computing price.

Different formulations of the same bound. Computing times depend obviously on the solvers, and in turn the performances of the solvers depend on the formulation of the bounds. In our case, the SDP bound can be computed by solving (9) or (10), and each of both turns out to advantage one solver.

To have a fair comparison between the three solvers, we have done preliminary tests (reported in [START_REF] Malick | Numerical study of SDP bounds for the k-cluster problem[END_REF]) to choose the best semidefinite formulation for each solver. As expected, (10) is better suited for CSDP (when n = 100 about five times faster than using (9)), since here the computing time directly depends on the number of variables and constraints. On the other hand, the formulation (9) provides the best results for SB and SDLS (when n = 100 about four times faster than using (10)). In the numerical experiments below, we thus have: SB solves (9), CSDP solves (10), and SDLS solves (14). Test problems. We consider instances of k-cluster, randomly generated as follows.

Given a density d ∈ [0, 1], a random number ρ ∈ [0, 1] is generated for any pair of indexes i < j; if ρ > d, then w ij is set to 0, otherwise it is set to 1. The parameters of these instances are: the size of the graph n, the graph density d and the value of k. In this section, we take n = 80, 100, 300, d = 25%, 50%, 75% and k = n/4, n/2, 3n/4.

We have 5 instances of k-cluster problems for each set of parameters. Each forthcoming numerical result is thus an averaged of 5 results. The instances with 80 and 100 vertices have already been used in several others papers, as [START_REF] Billionnet | Different formulations for solving the heaviest k-subgraph problem[END_REF] and especially [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] which is the reference method for the comparison of Section 5.2. Note that these problems have unweighted graphs, which are at least as hard as the ones with weighted graphs. Numerical results. The numerical experiments have been carried out on a Pentium IV 2.2 GHz with 1 Go of RAM under Linux. Table 1 reports the computation times (in seconds) of the three solvers to compute the SDP bounds. Figure 2 illustrates the run on two characteristic instances. The results show that Θ(α) provides a SDP-quality bound which is easier to get than the usual SDP bound Θ(0).

Since the solvers do not compute the same bounds (remember that SB and CSDP compute Θ(0) and SDLS computes Θ(α)), we report:

gap SDP -the relative difference between Θ(α) and Θ(0), i.e., (Θ(α)-Θ(0))/Θ(0), -time SDLS -the computing time for SB and CSDP to achieve the value Θ(α). Some comments are in order. Although less tight, the new bounds Θ(α) are very close to Θ(0): the relative gap always lower than 0.25%. Remember moreover that the gap will be further reduced within the branch-and-bound by rounding down the bound to the integer value. We notice also that the gap is almost constant in absolute value so that its relative value decreases mechanically, when k or d increase (since the optimal value of the instance increases in this case).

The crucial point is that our solver provides bounds faster than SB and CSDP. As expected, CSDP is very robust and efficient for the medium-size instances (n = 80); but compared to SB, its performance deteriorates when the size of the problem gets bigger. In this case, the interest of SDLS becomes clear: the gap between the bounds does not increase, whereas the running times increase slower than those of SB and CSDP.

Furthermore, our solver provides these bounds faster whenever the solving process is interrupted. This is illustrated in Figure 2 for example: the convergence curve is always below the one of CSDP and also the one of SB (excepted at the very beginning). Note also that none of SB and CSDP dominates the other.

We insist finally on the three following points.

1. Robustness. As expected, the running times CSDP are almost constant for the instances of the same size (see for example the two examples of Figure 2). Our solver has also a similar behavior, which is a desirable property in the context of branch-and-bound. Note that SB does not have this behavior (it was faster than CSDP on some instances, but slower on others). For example, let us give the standard computing time deviations for n = 100. More precisely, the mean standard deviations for the five instances at given triplet (n, k, d) are 0.02 for SDLS, 0.11 for CSDP, and 2.32 for SB. 2. Fast initial decrease. Our solver SDLS (as well as SB) has a very fast initial convergence, as illustrated on Figure 2. This is also a highly desirable property in the context of branch-and-bound. In particular, the improvement of the bound computed by SDLS is big in the first iterations. 3. Low-memory. We also mention that our solver uses little memory: for instance, it requires about 4 MB to solve (14) when n = 100.

The conclusion of these first experiments is that the new SDP bounds have an advantageous trade-off between tightness and speed of computation, together with interesting features in view of exact resolution.

Exact resolution: solving k-cluster to optimality

We use the solver to compute Θ(α) as the bounding procedure of a branch-and-bound for solving k-cluster problems to optimality. As presented in Section 5.1, our branchand-bound algorithm is very basic; the novelty is essentially the use of Θ(α). In Section 5.2, we compare it with [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] the best known method to solve k-cluster (that mixes nicely CPLEX and quadratic relaxations parameterized with SDP).

Branch-and-bound for k-cluster

We summarize in this section the main ingredients of our branch-and-bound algorithm.

Heuristics for initial feasible point. Good feasible points give good lower bounds at the beginning of the branch-and-bound, and this yields better pruning. Following [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF], we apply the following two-step heuristic that gives very good feasible points:

1. Apply the natural greedy algorithm (see e.g. [START_REF] Asahiro | Greedily finding a dense subgraph[END_REF]) consisting in removing n-k vertices from G by choosing successively the vertex with the minimum degree in the reduced graph. 2. Enhance the feasible point by a local search which consists in swapping two nodes (one in the k-cluster and not the other) until no improvement is possible.

Branching strategy. For simplicity, the branching scheme is based on a basic depthfirst search. Moreover, we have to choose, at each node of the search tree, a vertex to add to the cluster or to exclude from it (in other words, we want to fix one of the {-1, 1}-variables). Again for simplicity, we use a branching order fixed in advance, as follows. As in [START_REF] Pisinger | Upper bounds and exact algorithms for p-dispersion problems[END_REF], we estimate of the expected contribution of each node, by solving for each i

D i = max 8 < : X j∈V \{i} w ij z i j : X j∈V \{i} z i j = p -1, z i j ∈ {0, 1}, j ∈ V \{i} 9 = ;
. This problem asks to choose p-1 largest values among the edges connected to i -which can be done efficiently by a median search algorithm. Then we order the variables with decreasing D i . Numerical testing has showed us that this choice is competitive with other approaches (such as branching on the most fractional variable).

Branching step. Branching on the index i 0 consists in decomposing (14) into two smaller instances, as follows. To obtain the new objective function, we remove the i 0 th line from Q, and add it (or subtract if the vertex is excluded) to the first one. We do the same for the i 0 th column. We also remove two constraints: E i0 , X = 1 and Q i0 , X = 4k -2n. Finally we use the previous optimal y to make a warm-start for the two reduced problems: we keep the current values associated to the remaining constraints. We observed that this speeds up the solving process of a mean factor of 3.

Early termination. Within branch-and-bound algorithms, the early termination of the bounding procedure often provides a substantial overall speed-up (see e.g. [START_REF] Faye | A cutting planes algorithm based upon a semidefinite relaxation for the quadratic assignment problem[END_REF] and [START_REF] Rendl | Solving max-cut to optimality by intersecting semidefinite and polyedral relaxations[END_REF]). In our case, the next observation is important, so we formalize it as a lemma.

Lemma 2 For any variable y ∈ R 2n+2 and any α > 0, the value α((n + 1) 2 /2 + f (y)) is a bound for the k-cluster problem.

Proof Just combine ( 13) and (15) to get that val (1) ≤ Θ(α) ≤ α " (n + 1) 2 2 + f (y)

" .

Thus we get a bound of the optimal value of (1).

⊓ ⊔

As any differentiable optimization algorithm, the algorithm to minimize f is stopped when an approximate solution is computed; more precisely when ∇f (y k ) = AX(y k ) -b ≤ ε with the fixed tolerance ε. Though the computed X(y k ) does not satisfies perfectly the constraints, the lemma guarantees that we still get a bound for our problem with y k . In fact, we can stop the algorithm anytime before convergence, and we still get a bound. This is very useful within the branch-and-bound: we can stop the bounding procedure when the bound gets lower than the previous threshold. More precisely, we use the following early termination rule. Let β be the current best lower bound given by a k-subgraph; then we stop the run of the algorithm when (n + 1) 2 /2 + f (y) < (β + 1)/α.

Fixing α and ε. For simplicity, we keep the level-parameter α and the stopping criterion ε constant. We chose these parameters empirically: we did numerical tests on some instances (as those of Section 4, see also Figure 1), we observed that ε = 10 -7 and α = 10 -4 seem to be good values, and we take those values for all the tests. Dynamic adjustement of those parameters might provide additional speed-ups, but we postpone this technical point for further research.

Experiments on solving k-cluster to optimality: comparison

The best approach for k-cluster. As far as we are aware, the currently strongest results are obtained by the convex quadratic relaxation procedure of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] (see also [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF], and [START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF]), that we briefly recall here. Start with the initial formulation (2), and observe that for any u, γ, the problem is equivalent to

8 > < > : max z ⊤ Q(u, γ)z e ⊤ z = k z ∈ {0, 1} n (16)
of SDP-quality everywhere, without paying the whole price in computing time. Thus our method is particularly well-adapted for the instances with small d and k; these are the most difficult instances, where good pruning is essential. Moreover, the increase of number of nodes with increase of dimension is less dramatic for our method than for QCR (passing from n = 80 to n = 100, the number is multiplied by 2 instead of 3). Finally, we point out that the computing time of our method is not strongly correlated to the number of the evaluated nodes in the search tree. This is a consequence of our solver ability to be interrupted during the solving process.

Experiments on solving k-cluster to optimality: larger problems

We have conducted numerical experiments on problems beyond the actual state of the art. Table 3 reports the results for instances with n = 120. For those large graphs, we change the experimental protocol, as follows. We do not compare with QCR (the largest tests of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] and [START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF] are with n = 100). So we do need anymore to use the same computer as their: we run this third experiment on a Dell precision T7500 Intel Xeon 2.80GHz with 4GB RAM under linux, using singlethread (single core). We fix a computing time limit of 60000 seconds, and we report the number of instances solved for each setting (in the column "# pbs solved").

As far as we know, this is the first time that k-cluster problems with unrestricted graphs of size larger than 100 are solved to optimality. On the other hand, we see that our solver runs into some trouble, especially for the most difficult problems (those with low density). Two reasons could explain those difficuties, and open ways to improve the approach. First, recall that, except for the bounding procedure, our branch-andbound is very basic. Including more sophisticated strategies on each point of Section 5.1 would probably bring much performance improvement. Second, the bound Θ(α) is not tight enough to avoid losing oneself in the branch-and bound tree (see the gap at the root). Considering stronger bounds of the same kind (adding cuts for example) would therefore be interesting. Those two points are subject to current investigations.

Conclusions, perspectives

This paper develops a branch-and-bound algorithm using a novel bounding procedure to solve k-cluster problems to optimality. The key is to use new semidefinite bounds for k-cluster that trade some quality of bound for a speed-up of computation time. This is the first successful attempt to solve this NP-hard optimization to optimality with a semidefinite programming approach competitive with state-of-the-art methods.

The numerical experiments of Sections 4.2 and 5.2 show that the new semidefinite bounds θ(α) have a practical interest. They have indeed good balance between tightness and computing time; they provide SDP-quality bounds while being faster to compute. The dedicated solver for computing θ(α) also combines advantages of the SDP solvers SB and CSDP: like SB, it gives guaranteed upper bounds, has a fast initial convergence, and allows to be interrupted; and like CSDP it is reliable, in the sense that we observe only small variations in computational times. The branch-and-bound using θ(α) takes advantage of SDP-like bounds (all way long) to prune well. Its performance is comparable with the best method for this problem. In practice, our method works particularly fine on the most difficult instances of k-cluster (with a large number of vertices, small density and small k).

The exact resolution scheme presented here could be adapted to other combinatorial problems: the semidefinite bounds upon which the approach is based are indeed introduced in [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF] for general binary quadratic problems. The first step toward a generalization would be to extend the numerical study of Section 4.2 (and [START_REF] Malick | Numerical study of SDP bounds for the k-cluster problem[END_REF]); this is what proposes our recent work [START_REF] Malick | On the bridge between combinatorial optimization and nonlinear optimization: a family of semidefinite bounds leading to newton-like methods[END_REF] (which also presents a different derivation of the SDP bounds).

Lemma 3 With the notation of this section, the k-cluster problem (1) is equivalent to the quadratic problem in dimension n + 1 8 > < > :

max x ⊤ Q x x ⊤ Q j x = 4k -2n, j ∈ {0, . . . , n} x ∈ {-1, 1} n+1 (17)
where the symmetric (n + 1) × (n + 1)-matrices Q and Q j (j ∈ {0, . . . , n}) are defined by (5). More precisely, this equivalence means that the optimal values of (1) and (17) are the same and that the optimal solutions coincide as follows:

if z is a solution of (1) then x = (1, 2z -e) is a solution of (17); -if x is a solution of (17), then z = ((x 0 x1 , . . . , x0 xn) + e)/2 is a solution of (1).

Proof

The formulation of (3) is equivalent to (18) in the sense that the optimal values are the same and the solutions are in one-to-one correspondance with x = 2z -e.

We consider now that the following purely quadratic problem with the additional variable x 0 ∈ {-1, 1} 8 > > < > > : max 1 4 (x ⊤ W x + 2x ⊤ W e x 0 + e ⊤ W e x 0 2 ) e ⊤ x x 0 = 2k -n x ⊤ C j x + 2x ⊤ ẽj x 0 = 4k -2n, j ∈ {1, . . . , n} x i ∈ {-1, 1} i ∈ {0, . . . , n}.

We observe that this problem is equivalent to (17) in view of the definitions of the matrices in (5). So we just have to establish that (18) is equivalent to (19); we do so in two steps.

  but we get arbitrarily close to the SDP bound when α tends to 0 lim α→0,α>0 Θ(α) = val (9).

Fig. 1

 1 Fig.1Illustration of the features of Θ(α) on a graph with n = 300, k = 75 and d = 25%: when α → 0, we get more tightness (Θ(α) → Θ(0)) but for an increasing computing cost. The leftmost point corresponds to α = 10 -4 . .

Fig. 2

 2 Fig. 2 Comparison of the CPU time of the three solvers for two instances of the k-cluster problem with n = 300, d = 0.5, and k = 150 (left side) k = 225 (right side)

  Operate first the change of variable x = 2z -e, and express the objective and the constraints with respect to x = (x 1 , . . . , xn) ∈ R n . Just develop the objectivez ⊤ W z = (x + e) W x + 2x ⊤ W e + e ⊤ W e),and similarly transform e ⊤ z = k as e ⊤ x = 2k -n, and for all j = 1, . . . , nz ⊤ C j z = 2kz j ⇐⇒ x ⊤ C j x+2x ⊤ (C j e-2ke j ) = 4k-e ⊤ C j e ⇐⇒ x ⊤ C j x+2x ⊤ ẽj = 4k-2nthe last equivalence coming from e ⊤ C j e = 2n and ẽj = C j e-2ke j . So we get quadratic problem in x ∈ {-1, 1} n ⊤ W x + 2x ⊤ W e + e ⊤ W e) e ⊤ x = 2k -n x ⊤ C j x + 2x ⊤ ẽj = 4k -2n, j ∈ {1, . . . , n} x i ∈ {-1, 1} i ∈ {1, . . . , n}.

Table 2

 2 Computing times (in seconds) for exact resolution: comparison on the test-problems of[START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] and[START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF].

	n	k	d(%)		B&B using Θ(α)		B&B of [BEP09]
				time	# nodes	# pbs	gap(%)	time	# nodes
				(sec.)		faster	(root)	(sec.)
	80	20	25	186.6	14869	1/5	9.43	72.4	207562
			50	674.2	43702	0/5	8.19	325.6	527300
			75	1809.4	132455	0/5	6.60	1322.6	2170499
		40	25	28.5	1426	1/5	2.76	19.6	31442
			50	53.6	6774	0/5	1.86	24.8	36544
			75	413.4	19784	0/5	1.39	176.9	276198
		60	25	3.69	200	0/5	0.90	0.92	1699
			50	6.91	330	0/5	0.62	2.09	4358
			75	20.44	1098	0/5	0.43	3.60	7592
			mean	355.2	24515			216.5	362577
	100	25	25	1443 (3)	155658	2/3	9.45	1483 (3)	2499727
			50	1802 (2)	172912	2/2	8.10	1978 (1)	3510428
			75	n.a.	n.a.	n.a.	5.75	n.a.	n.a.
		50	25	314.2	13532	4/5	2.63	473.6	504644
			50	349 (1)	12020	1/1	2.23	392 (1)	396780
			75	704.2	29503	2/5	1.09	548.4	580512
		80	25	16.08	565	0/5	0.86	9.00	12996
			50	72.09	2547	0/5	0.68	50.81	74196
			75	40.49	2662	0/5	0.36	14.32	21045
			mean	592.6	48675			668.6	950041

Table 3

 3 Computing times (in seconds) for exact resolution of large-scale instances.

	n	k	d(%)		B&B using Θ(α)
				time	# nodes in	# pbs	gap (%)
				(sec.)	search tree	solved	(root)
	120	30	25	57156	964968	1/5	10.00
			50	n.a.	n.a.	0/5	n.a.
			75	n.a.	n.a.	0/5	n.a.
		60	25	16230	223187	3/5	3.26
			50	3716	56101	3/5	2.61
			75	9020	176355	5/5	1.47
		90	25	124	2952	4/5	1.07
			50	12818	179349	5/5	1.19
			75	1845	24549	5/5	1.03

  If (x 1 , . . . , xn) is feasible in (18), then obviously (1, x) and (-1, -x) are both feasible in (19), with same objective value. It follows that val (19) ≤ val (18). Conversely if (x 0 , . . . , xn) is feasible in (19), then x 0 = ±1 and (x 0 x 1 , . . . , x 0 xn) is feasible in (18) with the same objective value. It follows that val (18) ≤ (19), so that we have the equality in fact. The relation between the argmins then becomes clear. ⊓ ⊔

Acknowledgement

This work was supported by CNRS (through "GdR Recherche Opérationnelle") and Grenoble University (Université Joseph Fourier, through "Pôle Math-STIC"). We thank Quentin Monnet and Lise-Marie Veillon (Master students of ENSIIE, Evry, France) for their valuable help in developing parts of the solver. We also thank the two anonymous referees for numerous suggestions about a preliminary version of this article.

with the quadratic form

Relaxing the integrality constraints gives the upper bound

If u and γ are chosen such that Q(u, γ) is negative semidefinite, the bound R(u, γ) is obtained by solving a convex quadratic problem. In practice, the variables (u * , γ * ) are chosen such that they minimize the upper bound

which turns out to be a linear SDP problem, more precisely, the dual of (9). Thus a solution of the problem ( 16) with (u * , γ * ) is computed by a standard solver for convex quadratic 0-1 problem, more precisely a branch-and-bound algorithm using R(u * , γ * ) as bounding procedure. The interest is then that they can advantageously use CPLEX implementing the state-of-the-art for mixed-integer quadratic solver.

Numerical results. We solve to optimality, on the same computer, the same testproblems as in [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] (see also [START_REF] Plateau | Reformulations quadratiques convexes pour la programmation quadratique en variables 0-1[END_REF]). Table 2 presents the numerical results: we can solve all the instances with comparable times; since we have implemented a very basic branch-and-bound (simple branching strategy and depth-first search), whereas the results of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] are obtained with the mixed-integer programming solver of CPLEX, this shows the interest of our approach. Table 2 reports for each method the computing time and the number of nodes. Recall that each entry of the table is the mean over the 5 instances with the same parameters settings. To go beyond average numbers, the table also shows the number of problems (out of 5) solved faster by our method (in the column "# pbs faster").

Let us comment more precisely the comparison of Table 2. We observe that QCR is faster in average than our approach for graphs with n = 80 vertices but it is slower when n = 100. As in [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF], the computation is stopped after one hour: we indicate 'n.a.' if the corresponding problems were not solved within the time limit, and we indicate '(x)' next to the computing time when we were able to solve only 'x' instances (out of 5) in less within the time limit.

We did not include the results for small graphs (n = 40) that are in [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF]: for those graphs, QCR clearly outperforms our method (by a factor between 5 and 10), but both methods converge very quickly (less than few seconds). In that case, the tightness of the bound is actually not crucial, and it seems like there is no need to spend time in computing semidefinite bounds.

Though the two methods are comparable with respect to computing times, they have opposite strategies. QCR does many (cheap) bound evaluations, whereas our algorithm computes more expensive bounds and prunes better: we enumerate 15-20 times fewer nodes on average. The bound of [START_REF] Billionnet | Improving the performance of standard solvers for quadratic 0-1 programs by a tight convex reformulation: The QCR method[END_REF] is of SDP-quality at the root of the branch-and-bound tree only, but deteriorates down in the tree, while Θ(α) is Appendix: Formulation of the k-cluster as a quadratic {-1, 1}-problem

The spherical constraints appears more easily on purely quadratic {-1, 1}-optimization problems (see [START_REF] Malick | Spherical constraint in Boolean quadratic programming[END_REF]). The second reformulation in Section 2 considers the transformation of the k-cluster problem, from the natural modeling as a quadratic {0, 1}-problem (with linear and quadratic constraints), to a purely quadratic {-1, 1}-problem. We specify here this transformation, that uses standard techniques.


Emidio ; Scarpellini

Kathleen Blondeau

Veerle Boecxstaens

Rita Vos

Antonio Gasbarrini

Ricard Farre

M.D Jan Tack

R Farré

EFFECT OF RIMONABANT ON ESOPHAGEAL MOTOR FUNCTION IN MAN

Keywords: GERD or GORD < Disease-based, Oesophagus < Organ-based, Motility < Topics, Acidity (oesophageal) rimonabant, cannabinoid-1 receptors, transient lower esophageal sphincter relaxations, gastroesophageal reflux disease

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Gastroesophageal reflux disease (GERD) is defined as the presence of symptoms or lesions that can be attributed to the reflux of gastric contents into the esophagus. [START_REF] Vakil | The Montreal definition and classification of gastroesophageal reflux disease: a global evidence-based consensus[END_REF][START_REF] Tack | Recent developments in the pathophysiology and therapy of gastroesophageal reflux disease and non-erosive reflux disease[END_REF] The pathophysiology of GERD is multifactorial and involves several mechanisms such as failure of the anti-reflux barrier, impaired esophageal clearance, the presence of offensive factors in the refluxate and defective esophageal tissue resistance [START_REF] Koek | Multivariate analysis of the association of acid and duodeno-gastro-oesophageal reflux exposure with the presence of oesophagitis, the severity of oesophagitis and Barrett's oesophagus[END_REF]. Among the dysfunctions of the antireflux barrier, transient lower oesophageal sphincter relaxations (TLESRs) are the major mechanism underlying gastro-oesophageal reflux events in normal subjects, and in most GERD patients [START_REF] Dent | Mechanism of gastroesophageal reflux in recumbent symptomatic subjects[END_REF][START_REF] Dodds | Mechanisms of gastroesophageal reflux in patients with reflux esophagitis[END_REF][START_REF] Dent | Mechanisms of lower esophageal sphincter incompetence in patients with symptomatic gastroesophageal reflux[END_REF]. TLESRs are a vago-vagal reflex, triggered by activation of stretch receptors of the proximal stomach, and organized in the brain stem. [START_REF] Dodds | Mechanisms of gastroesophageal reflux in patients with reflux esophagitis[END_REF][START_REF] Dent | Mechanisms of lower esophageal sphincter incompetence in patients with symptomatic gastroesophageal reflux[END_REF][START_REF] Mittal | Characteristics and frequency of transient relaxations of the lower esophageal sphincter in patients with reflux esophagitis[END_REF] Pharmacological restoration of the anti-reflux barrier is considered a relevant target for the control of GERD [START_REF] Achem | A prokinetic approach to treatment of gastroesophageal reflux disease[END_REF]. Several gastroprokinetic agents which also enhance lower esophageal sphincter (LES) pressure and esophageal motility have been developed and used in the treatment of GERD [START_REF] Tack | Recent developments in the pathophysiology and therapy of gastroesophageal reflux disease and non-erosive reflux disease[END_REF][START_REF] Achem | A prokinetic approach to treatment of gastroesophageal reflux disease[END_REF] Because of side-effects or lack of efficacy, none of these drugs has led to a currently established treatment for GERD. Drugs inhibiting TLESRs are considered a novel and attractive target for GERD therapy. [START_REF] Blackshaw | Inhibition of transient LES relaxations and reflux in ferrets by GABA receptor agonists[END_REF][START_REF] Lehmann | Activation of the GABAB receptor inhibits transient lower esophageal sphincter relaxations in the dog[END_REF][START_REF] Lidums | Control of transient lower esophageal sphincter relaxations and reflux by the GABAB agonist baclofen in normal subjects[END_REF] Baclofen, an agonist at the gamma-aminobutyric acid receptor type B (GABA-B), has been demonstrated to inhibit considerably the occurrence of TLESRs and gastro-oesophageal reflux, with an acceptable side-effect profile [START_REF] Lidums | Control of transient lower esophageal sphincter relaxations and reflux by the GABAB agonist baclofen in normal subjects[END_REF][START_REF] Zhang | Control of transient lower esophageal sphincter relaxations and reflux by the GABAB agonist baclofen in patients with gastroesophageal reflux disease[END_REF][START_REF] Van Herwaarden | The effect of baclofen on gastro-oesophageal reflux, lower oesophageal sphincter function and reflux symptoms in patients with reflux disease[END_REF][START_REF] Lee | Differential effects of baclofen on lower oesophageal sphincter pressure and proximal gastric motility in humans[END_REF]. It is currently thought that the antireflux effect of baclofen depends on peripheral inhibition of mechanosensitive gastric vagal afferents, thereby raising the threshold for action potential firing, and partly on central inhibition of vagal afferents leading to reducedtransmitter release [START_REF] Partosoedarso | GABA(B) receptors on vagal afferent pathways: peripheral and central inhibition[END_REF]. Cannabinoid receptors (CBRs) have many similarities with GABA-B receptors. For example agonism of both GABA-B receptors and CBRs inhibits transmitter release presynaptically and postsynaptically hyperpolarizes neurons in the central nervous system 5CNS) [START_REF] Bowery | Gamma-Aminobutyric acid B receptors: first of the functional metabotropic heterodimers[END_REF][START_REF] Pertwee | Pharmacology of cannabinoid CB1 and CB2 receptors[END_REF]. Common CBR-mediated actions in vivo include skeletal muscle relaxation, hypothermia, and antinociception [START_REF] Baker | Cannabinoids control spasticity and tremor in a multiple sclerosis model[END_REF]. CBRs have an important role in inhibition of emesis triggered by stimuli in the gastrointestinal tract or CNS [START_REF] Simoneau | The cannabinoid agonist WIN55,212-2 suppresses opiod-induced emesis in ferrets[END_REF][START_REF] Van Sickle | Cannabinoids inhibit emesis through CB1 receptors in the brainstem of the ferret[END_REF][START_REF] Van Diepen | Prejunctional and peripheral effects of the cannabinoid CB(1) receptor inverse agonist rimonabant (SR 141716)[END_REF]. They have been localized to specific points along the emetic pathway in the central nervous system that receive input from gastrointestinal vagal afferents. There is considerable anatomical overlap of Rimonabant and esophageal function Scarpellini et al.,[START_REF] Dent | Mechanism of gastroesophageal reflux in recumbent symptomatic subjects[END_REF] pathways involved in emesis and triggering of TLESRs [START_REF] Hornby | Central neurocircuitry associated with emesis[END_REF][START_REF] Mittal | Transient lower esophageal sphincter relaxation[END_REF]. CBRs, like GABA-B receptors, belong to the G protein-coupled receptor superfamily and are subdivided into 2 groups, CBR1 and CBR2 [START_REF] Howlett | International Scarpellini et al., 15 Union of Pharmacology. XXVII. Classification of cannabinoid receptors[END_REF]. Delta-9-tetrahydrocannabinol (delta(9)-THC) is the principal psychoactive compound of marijuana and is well established as an agonist at CBRs. More recently, endogenous ligands (endocannabinoids) have been isolated, the most important of which are arachidonylethanolamide (anandamide) and 2-arachidonyl glycerol [START_REF] Howlett | International Scarpellini et al., 15 Union of Pharmacology. XXVII. Classification of cannabinoid receptors[END_REF][START_REF] Diana | Short-term retrogradeinhibition of GABAergic synaptic currents in rat Purkinje cells is mediated by endogenous cannabinoids[END_REF][START_REF] Wilson | Endogenous cannabinoids mediate retrograde signalling at hippocampal synapses[END_REF]. The endocannabinoids appear to act as retrograde messengers, i.e. they can be released by postsynaptic neurons and diffuse to nerve terminals where they reduce transmitter release [START_REF] Diana | Short-term retrogradeinhibition of GABAergic synaptic currents in rat Purkinje cells is mediated by endogenous cannabinoids[END_REF][START_REF] Wilson | Endogenous cannabinoids mediate retrograde signalling at hippocampal synapses[END_REF]. A number of synthetic ligands for CBRs have been made, and whereas CBR agonists are moderately subtype selective, antagonists show definitive selectivity [START_REF] Wilson | Endogenous cannabinoids mediate retrograde signalling at hippocampal synapses[END_REF]. The CB1 antagonist rimonabant was commercially available as a drug in the treatment of obesity, until being expelled from the market because of an increased prevalence of depression. [START_REF] Currioni | Rimonabant for overweight or obesity. Cochrane database of systematic reviews[END_REF][START_REF] Akbas | A critical review of the cannabinoid receptor as a drug target for obesity management[END_REF].

In dogs and in humans, delta(9)-THC significantly reduced resting LES pressure and the occurrence of TLESRs after a standard meal [START_REF] Beaumont | Effect of delta9-tetrahydrocannabinol, a cannabinoid receptor agonist, on the triggering of transient lower oesophageal sphincter relaxations in dogs and humans[END_REF]. A study in dogs established that CB1 agonists inhibit the occurrence of TLESRs while rimonabant reversed this effect [START_REF] Lehmann | Cannabinoid receptor agonism inhibits transient lower esophageal sphincter relaxations and reflux in dogs[END_REF]. The effects of rimonabant on LES function in man has not been addressed so far. [START_REF] Scarpellini | Effect of itopride one esophageal motitiy and lower esophageal sphincter function in man[END_REF].

RECORDING METHODS

Following an overnight fast, an oesophageal manometric catheter fitted with a 6-cm Dent Sleeve was introduced through the mouth. Subsequently, the oesophageal catheter was positioned so that pressures could be recorded from the fundus (side hole 2 cm below the sleeve), the lower oesophageal sphincter (sleeve), oesophageal body (side holes 4, 7 and 10 cm proximal to the sleeve) and pharynx (side hole 28 cm proximal to the sleeve, to detect swallows). The oesophageal catheter was infused at a flow rate of 0.5 mL/min with distilled water using a low-compliance pneumo-hydraulic capillary infusion system (Arndorfer Medical Specialties, Milwaukee, WI, USA). The infusion system was connected to external pressure transducers, and signals were recorded on a polygraph (Synectics Medical, Stockholm, Sweden).

The oesophageal pH was measured with an antimony pH electrode (Synectics Medical, Stockholm, Sweden) positioned 5 cm above the proximal margin of the sleeve. The pH electrode was calibrated in buffers of pH 1 and pH 7 before and after each study. During the study period, the oesophageal pH was recorded continuously using an ambulatory data-logger (MicroDigitrapper, Synectics Medical, Stockholm, Sweden).

DATA ANALYSIS

Lower oesophageal motility

The tracings were analysed in a blinded fashion by one of the authors (ES) and checked by a second readere (JT). The basal lower oesophageal sphincter pressure was measured at end expiration relative to the end-expiratory intragastric pressure. The basal lower oesophageal sphincter pressure was visually determined every 3 min and averaged over

Oesophageal pH

The percentage of time with an oesophageal pH < 4 and the number of acid reflux episodes were calculated. Acid reflux episodes were defined as a decrease in oesophageal pH to a value below pH 4 for at least 4 s or, if the basal oesophageal pH was already below pH 4, as a rapid further drop in pH of at least 1 pH unit.

Statistical analysis

Data are presented as the mean ± standard error of the mean (S.E.M.). The paired t-test was used for the comparison of the mean values between the periods or between the rimonabant and placebo studies. The changes in the basal lower oesophageal sphincter pressure, was evaluated using analysis of variance for repeated measures. The frequency and duration of transient lower oesophageal sphincter relaxations, the number of acid reflux episodes and the symptom scores were compared using the Wilcoxon signed rank test. A p-value <0.05 was considered to be statistically significant. Based on previous studies, the study had a 85% power to detect 30% difference in TLESR rate at 5% significance level [START_REF] Lee | Differential effects of baclofen on lower oesophageal sphincter pressure and proximal gastric motility in humans[END_REF][START_REF] Scarpellini | Effect of itopride one esophageal motitiy and lower esophageal sphincter function in man[END_REF].

RESULTS

Conduct of the study

The positioning of oesophageal manometric catheter and pH probe were well tolerated, and all subjects completed both the sessions of studies. No adverse events were reported

Lower esophagal sphincter pressure

Prior to drug administration, and in the preprandial period prior to the meal, LES resting pressure and swallow-induced relaxations were similar for both conditions (Table 1,2). In the placebo studies, ingestion of the meal was associated with a significant decrease in the LES pressure during both the first and the second postprandial hour (Figure 2). In contrast, in the rimonabant studies, no significant decrease in LES pressure occurred after the meal (Table 1).

Esophageal motility

The amplitude of peristaltic contractions and swallow-induced relaxations were not significantly altered by rimonabant in the pre-prandial and post-prandial periods. However, rimonabant significantly increased the duration of peristaltic contractions in the preprandial and post-prandial periods (Table 1). Rimonabant did not alter the amplitude of swallow-induced relaxations, but the duration of contractions was significantly increased by rimonabant (Table 1).

Transient lower esophageal sphincter relaxations.

The numbers of transient lower oesophageal sphincter relaxations after the administration of placebo and rimonabant are summarized in Figure 2

Esophageal pH

The percentage of time pH < 4 in the oesophagus did not differ between the rimonabant and placebo studies in the pre-prandial (0.1± 0.05 % vs. 0.3±0.1%, NS) but tended to be lower after rimonabant in the postprandial period (0.1± 0.05 vs. 0.4±0.1%; p= 0.06). The number of acid reflux episodes increased significantly after the meal in the placebo studies but not in the rimonabant studies (Figure 4). The number of acid reflux episodes during the post-prandial period after rimonabant was significantly lower than after placebo (Figure 4).

Symptoms

No significant differences in symptom scores (calculated as area under the curve, AUC)

during both the pre-prandial and post-prandial periods were found between rimonabant and placebo studies (details not shown).

DISCUSSION

The aim of the present study was to investigate how suppression of endocannabinoid action by rimonabant would influence esophageal motility, the occurrence of TLESRs elicitation and gastro-esophageal reflux events in healthy volunteers. We found that rimonabant inhibited the decrease in LES pressure and the increase in the rate of TLESRS and reflux events after a meal, leading to a significant lower number of acid reflux events.

In addition, the duration of peristaltic contractions in the distal esophagus was prolonged after rimonabant.

The observation that rimonabant increased LES pressure after a meal is in line with previous studies which demonstrated that the cannabinoid agonist delta(9)-THC reduced resting LES pressure after a standard meal in man [START_REF] Beaumont | Effect of delta9-tetrahydrocannabinol, a cannabinoid receptor agonist, on the triggering of transient lower oesophageal sphincter relaxations in dogs and humans[END_REF]. The inhibition of TLESRs and reflux events we observed is more surprising, as delta( 9)-THC administration inhibited the occurrence of TLESRs after a standard meal in healthy volunteers [START_REF] Beaumont | Effect of delta9-tetrahydrocannabinol, a cannabinoid receptor agonist, on the triggering of transient lower oesophageal sphincter relaxations in dogs and humans[END_REF], and in the dog rimonabant enhanced the rate of TLESRs and reflux events induced by an acidified meal and intra-gastric air insufflation [START_REF] Lehmann | Cannabinoid receptor agonism inhibits transient lower esophageal sphincter relaxations and reflux in dogs[END_REF]. These observations which indicate ongoing suppression of TLESR by a CB1 receptor in the dog [START_REF] Lehmann | Cannabinoid receptor agonism inhibits transient lower esophageal sphincter relaxations and reflux in dogs[END_REF], were not confirmed in the present human study.

The reason for the apparent species difference between dog and man, and the mechanism underlying the inhibition of TLESRs and reflux events by rimonabant remains to be established. The dosing regimen differed strongly between the dog studies and the present studies: dogs received a unique dose of 0.22 umol/kg while the volunteers were treated for 3 days of medication with 20 mg/day. Furthermore, binding studies show that rimonabant has an affinity for the CB1 receptor in the low nanomolar range, but in functional experiments, rimonabant is not a neutral antagonist but, has rather been found to be an inverse agonist at the CB1 receptor [START_REF] Meschler | Inverse agonist properties of N-(piperidin-1-yl)-5-(4-chlorophenyl)-1-(2, 4-dichlorophenyl)-4-methyl-1H-pyrazole-3-carboxamide HCl (SR141716A) and 1-(2-chlorophenyl)-4-cyano-5-(4methoxyphenyl)-1H-pyrazole-3-carboxyl ic acid phenylamide (CP-272871) for the CB(1) cannabinoid receptor[END_REF][START_REF] Thomas | Cannabidiol displays unexpectedly high potency as an antagonist of CB1 and CB2 receptor agonists in vitro[END_REF]. The effects of rimonabant in the present study therefore can reflect antagonist properties, inverse agonist effects, or both.

TLESRs are controlled by a vago-vagal reflex pathway which is triggered by gastric distention, integrated in the brainstem and which induces release of nitric oxide from intrinsic nerves at the LES [START_REF] Tack | Anti-relaxation therapy in GORD[END_REF][START_REF] Tack | Recent developments in the pathophysiology and therapy of gastroesophageal reflux disease and nonerosive reflux disease[END_REF][START_REF] Boeckxstaens | Reflux inhibitors: a new approach for GERD?[END_REF] afferents, in the brain stem, on interneurons in the wall of the gastrointestinal tract and on postganglionic fibers innervating smooth muscle in the gastrointestinal tract [START_REF] Pertwee | Cannabinoids and the gastrointestinal tract[END_REF][START_REF] Duncan | Review article: endocannabinoids and their receptors in the enteric nervous system[END_REF][START_REF] Hornby | Involvement of cannabinoid receptors in gut motility and visceral perception[END_REF][START_REF] Boesmans | Cannabinoid receptor 1 signalling dampens activity and mitochondrial transport in networks of enteric neurones[END_REF]. In the dog, rimonabant only increased the number of TLESRs while their characteristics and esophageal motility were unchanged, and this was considered an argument in favor of a central site of action [START_REF] Lehmann | Cannabinoid receptor agonism inhibits transient lower esophageal sphincter relaxations and reflux in dogs[END_REF]. The observation that fasting LES pressure was not affected by rimonabant in the present study, while the postprandial drop in LES pressure and the occurrence of TLESRs were inhibited could also argue in favor for a central action, for instance on vago-vagal reflex pathways. Rimonabant also has the potential to inhibit the occurrence of TLESRs at the level of triggering through distention of the proximal stomach.

In a gastric barostat study, we demonstrated a significant decrease in postprandial volume of the proximal stomach in healthy volunteers after rimonabant pretreatment [START_REF] Ameloot | Endocannabinoid control of gastric sensorimotor function in man[END_REF]. In the present study in man, the duration of distal esophageal contractions was prolonged and the duration of TLESRs was shortened by rimonabant. Another important site of action for rimonabant, therefore, to consider is the enteric nervous system. Indeed, we and others previously provided evidence for a continuous endocannabinoid tone acting on CB1 to submit neurotransmitter release in the enteric nervous system [START_REF] Pertwee | Cannabinoids and the gastrointestinal tract[END_REF][START_REF] Duncan | Review article: endocannabinoids and their receptors in the enteric nervous system[END_REF][START_REF] Hornby | Involvement of cannabinoid receptors in gut motility and visceral perception[END_REF][START_REF] Boesmans | Cannabinoid receptor 1 signalling dampens activity and mitochondrial transport in networks of enteric neurones[END_REF]. Rimonabant may thus enhance the release of for instance acetylcholine from intrinsic motor neurons, which could contribute to the observed longer duration of peristaltic contractions and higher postprandial LES pressure. The absence of any change in symptoms from the gastrointestinal tract could also be viewed as an argument against a primarily central site of action. Further clarifying the site of action where rimonabant inhibits TLESRs and reflux would require studies with a peripherally acting CB1 antagonist, and no such drug is presently available for human studies.

Inhibition of TLESRs is a now well-established therapeutic target in GERD, and several drugs are being investigated for their ability to inhibit TLESRs, including GABA-B agonists, metabotropic glutamate-5 receptor antagonists and the mixed dopamine-2 antagonist/cholinesterase inhibitor itopride (9-15, 32, 44-46). The magnitude of inhibition of TLESRs and reflux events observed with rimonabant is at least in the same order of magnitude as observed with these agents. On the other hand, long-term administration of rimonabant has been associated with an increased occurrence of depression [START_REF] Moreira | Central side-effects of therapies based on CB1 cannabinoid receptor agonists and antagonists: focus on anxiety and depression[END_REF] be considered for this and other indications, provided that the site of action in TLESR inhibition is outside the blood-brain barrier [START_REF] Pavón | Central versus peripheral antagonism of cannabinoid CB1 receptor in obesity: effects of LH-21, a peripherally acting neutral cannabinoid receptor antagonist, in Zucker rats[END_REF][START_REF] Gómez | A peripheral mechanism for CB1 cannabinoid receptor-dependent modulation of feeding[END_REF].

The present study has a number of limitations. First, only a single dose of rimonabant was studied, as this was the only one available for the treatment of obesity in man. Second, the number of TLESRs and reflux events was low, which may be due to the fact that we studied asymptomatic healthy volunteers, and due to the relatively small size of the test meal. Finally, reflux events were only assessed by pH monitoring. It is conceivable that the use of esophageal impedance monitoring might have increased the number of detectable and quantifiable reflux events, especially after the meal.

In summary, in this placebo-controlled, double blind, randomized, crossover study we demonstrated for the first time that the CB1 receptor antagonist, rimonabant, inhibits the meal-induced increase in TLESRs, increases postprandial LES pressure leading to a lower number of acid reflux events. In addition, rimonabant increased the duration of distal esophageal peristaltic waves. Further studies will need to address whether this effect of rimonabant has a central or a peripheral site of action. and pH measurement studies after 3 days premedication with rimonabant 20 mg or placebo. After placement of the assembly ten wet swallows of 5 ml of water will be administered, followed by ingestion of the medication. After 60 minutes a standardized meal will be administered and measurements will continue for another 120 minutes. At 30minute intervals, 10 wet swallows will be administered. Throughout the study, at 15-minute intervals the intensity of 8 epigastric symptoms will be scored on visual analogue scales. We included a total of 4 figures, some of which are composite figures. We added individual graphs for the most relevant of these, see figures 2B and 3B.

F

Some missing references:

The references (41 and 45, now 32 and 43) have been updated. The authors have performed a small study in healthy volunteers demonstrating that rimonabant does have some potential effects on esophageal physiology.

It is unclear how the sample size was arrived at as the statistical analysis describes no power calculation. It is also unclear whether the tracings were read in blinded fashion, though I suspect that they were.

We made theses additions to the text:

1) Based on previous studies, the study had 85% power to detect 30% difference in TLESR rate at 5% significance level. We apologise for an error in the manuscript: there were 12 subjects, not 11. This has been corrected.

2) The tracings were read in blinded fashion by two of the authors.

The authors do well to point out the study limitations and the issues with the compound in general as well as pointing out that the magnitude of inhibition of TLESRs and reflux events similar to other agents being investigated for antireflux therapy. Early data from the reflux inhibitors cited have been underwhelming. Is there anything about rimonabant that makes it more attractive?

Together with the safety profile of the drug, the findings of the present study would not justify therapeutic development for GERD at this point, in our opinion. We have clarified this in the manuscript.

Table 2 could be eliminated from the manuscript if space were an issue. 2. It is not clear from the text (page 5), whether only the patients who were pretreated with this drug received additional dose of rimonabant during the manometry session, or it was blind?

We would prefer to keep

On the day of experiment, witnessed intake occurred in the motility unit. This was added to the text.

3. Methods-please include reference for VAS for specific symptoms mentioned in paper.

The reference for the VAS questionnaire was added (Scarpellini et al.,now ref. 31). 1. Despite the pretreatment with rimonabant the baseline manometric measurements did not differ between drug and placebo group. Does it mean that the drug has only short term effect or another explanation exists-needs a comment in discussion. The authors should comment practical clinical applicability based on results of the study.

Based on the pharmacokinetics of the drug, a steady state plasma level is expected a the time of the measurements. The lack of an effect on LES pressure before the meal suggests that the effect of rimonabant on LES pressure is only present in the postprandial state, while preprandial LES pressure is not affected. In contrast to the effect on LES pressure, there is already a preprandial effect of rimonabant on the duration of esophageal body peristalsis, which is in line with the presence of effective plasma levels before drug administration on the day of the experiments.

2. Although rimonabant has some interesting physiological effects on LES which might be helpful for GERD, low safety profile of rimonabant will most probably preclude the clinical use of this medication. This should be mentioned in discussion.

Thus, the future research should be directed to find safer alternative for this medication.

We agree with the reviewer. The poor safety profile for the drugs does not allow new studies in patients to verify the real impact of the present findings and requires the development of new drugs selectively acting on the peripheral CB1 receptor, or centrally acting CB1 receptor antagonists that do not induce depression; although the latter is unlikely to occur. This has been added to the discussion section. Tables&figures.

I think Table 2 is unnecessary. Table 1 is difficult to read. I would suggest to make three columns with physiologic parameters next to each other (LES pressure, Amplitude, and Duration) divided into Placebo/ Rimonabant , while the timing mentioned in first left column only once, instead of three times.

We would prefer to keep Table 2 in the manuscript: inhibition of TLESRs without interference with swallow-induced relaxation is important as it points to involvement of CB1 receptors in neural pathways controlling TLESRs, but not in swallow-induced peristalsis. Inhibition of swallow-induced relaxation is undesirable as this might lead to dysphagia. We have adapted Tables 1 and2 as suggested.

  Rimonabant and esophageal functionScarpellini et al.,[START_REF] Mittal | Characteristics and frequency of transient relaxations of the lower esophageal sphincter in patients with reflux esophagitis[END_REF] 15-min intervals. The influence of drug administration on the basal lower oesophageal sphincter pressure was assessed by comparing the value of the first with the value of the third pre-prandial 30-min interval.Transient lower oesophageal sphincter relaxations were defined according to criteria previously published[START_REF] Holloway | Criteria for objective definition of transient lower esophageal sphincter relaxation[END_REF] as follows: (i) absence of a swallowing signal for 4 s before to 2 s after the onset of lower oesophageal sphincter relaxation; (ii) relaxation rate of ≥ 1 mmHg/s; (iii) time from onset to complete relaxation of ≤ 10 s; and (iv) nadir pressure of ≤ 2 mmHg. Excluding multiple swallows, lower oesophageal sphincter pressure falls that fulfil the last three criteria, but have a duration of > 10 s, can also be classified as transient lower oesophageal sphincter relaxations irrespective of the timing of lower oesophageal sphincter relaxation relative to swallowing.

  Rimonabant and esophageal functionScarpellini et al.,[START_REF] Achem | A prokinetic approach to treatment of gastroesophageal reflux disease[END_REF]

  Rimonabant and esophageal functionScarpellini et al.,[START_REF] Lehmann | Activation of the GABAB receptor inhibits transient lower esophageal sphincter relaxations in the dog[END_REF]

Figure 2 .

 2 Figure 2. LES pressure. (A) Basal LES pressure was comparable for both conditions. Rimonabant pretreatment prevented the meal-induced decrease in LES pressure. * P < 0.05 compared to placebo. (B) Individual results for the LES pressure during the first postprandial hours.

Figure 3 .

 3 Figure 3. TLESRs. (A) Basal LES pressure was comparable for both conditions. Rimonabant pretreatment prevented the meal-induced decrease in LES pressure. * P < 0.05 compared to placebo. (B) Individual results for the total number of postprandial TLESRs.

Figure 4 .

 4 Figure 4. Acid reflux episodes. The number of acid reflux episodes during both the preprandial and post-prandial period was significantly lower after rimonabant. The pretreatment with rimonabant was not associated with a post-prandial increase of acid reflux episodes as after placebo. * p<0.05 compared to placebo. † p<0.05 compared to preprandial.

3 .

 3 Please comment on the difference of results of 3-days pretreatment and acute dose of drug during manometry.

  See above: in terms of duration of peristaltic congtractions, there is already an effect of rimonabant prior to drug administration on the day of the experiments. The effect on LES pressure is only present postprandially, suggesting a selective effect on postprandial LES motor control.

Table 3 )

 3 . Compared to placebo, the number of TLESRs in the rimonabant studies was significantly lower for the whole postprandial period and in the first postprandial hour.

	F o r
	P
	e
	e r
	R
	e
	v i e
	w

and Table

3

. After placebo, ingestion of the meal was associated with a significant increase in the rate of TLESRs during the first and the second post-prandial hour. In contrast, no significant increase in the rate of TLESRs relaxations occurred postprandially in the rimonabant studies (Figure

Table 1 .

 1 Esophageal motility parameters, before and after intake of placebo or rimonabant, and before and after the meal. *p<0.05 compared to placebo; **p<0.001 compared to placebo; † p<0.05 compared to basal.

	o r
	P
	e
	e r
	R
	e
	v i e
	w

Table 2 .

 2 Characteristics of swallow-induced relaxations before and after intake of placebo or rimonabant, and before and after the meal. No significant differences occurred between groups or over time.

		Relaxation (%)		Duration (sec)
		Placebo	Rimonabant 20	Placebo	Rimonabant 20
			mg		mg
	Basal	98.0±0.9	95.6±2.9	11.8±0.5	12.2±0.4
	Post-drug	97.0±1.3	97.1±1.3	12.5±0.7	12.2±0.4
	Postprandial	94.6±1.8	96.5±1.7	11.5±0.7	12.1±0.3
	1 st hour
	Postprandial	95.8±1.8	96.7±1.8	11.8±0.6	12.7±0.4
	2 nd hour

Table 3 .

 3 Number of transient lower oesophageal sphincter relaxation. † P < 0.05 compared to preprandial. *p< 0.05 compared to placebo.

	Pretreatment	Pre-prandial	First post-	Second post-	Total post-
			prandial hour	prandial hour	prandial
	Placebo	0.3±0.1	1.4±0.2 †	1.7±0.4 †	3.1±0.5 †
	Rimonabant 20 mg	0.4±0.2	0.4±0.2*	0.8±0.4	1.2±0.5*

FIGURES

Figure 1. Study outline. Healthy volunteers underwent esophageal sleeve manometry

Table 2

 2

	in the manuscript: inhibition of TLESRs without
	interference with swallow-induced relaxation is important as it points to involvement
	of CB1 receptors in neural pathways controlling TLESRs, but not in swallow-induced
	peristalsis. Inhibition of swallow-induced relaxation is undesirable as this might lead
	to dysphagia.

The rationale for the

[START_REF] Koek | Multivariate analysis of the association of acid and duodeno-gastro-oesophageal reflux exposure with the presence of oesophagitis, the severity of oesophagitis and Barrett's oesophagus[END_REF]

days pretreatment with rimonabant is based on the SmPC information. The bioavailability of rimonabant is not well established. The time to maximum concentration (Tmax) is approximately two hours, and the peak concentration (Cmax) is 196 ± 28.1 ng/ml. The drug's half-life is longer in obese patients (approximately 16 days) than in non-obese patients (nine days). Steady state is reached after 3 days).

Rimonabant and esophageal functionScarpellini et al., 2

ACKNOWLEDGEMENTS

Ricard Farré and Kathleen Blondeau are a postdoctoral research fellows of the FWO Flanders. Emidio Scarpellini is a Rome Foundation research fellow. This work was supported by a Methusalem grant to Jan Tack, M.D., Ph.D.


Louis Daniault
email: louis.daniault@institutoptique.fr

Marc Hanna

Laurent Lombard

Yoann Zaouter

Eric Mottay

Didier Goular

Pierre Bourdon

Frédéric Druon

Patrick Georges

Coherent beam combining of two femtosecond fiber chirped-pulse amplifiers

Ytterbium-doped fiber amplifiers present great advantages in femtosecond regime, including high optical-tooptical efficiency, good thermal handling and high spatial beam quality. However, the tight confinement of the beam inside the fiber and the long interaction length can lead to strong nonlinear effects at high energies, degrading the temporal profile of the output pulses. This constitutes a limit to femtosecond fiber chirped-pulse amplifier (CPA) setups.

A widely investigated way to reduce nonlinearities consists in enlarging the mode field area of the fiber amplifier, so that nonlinear effects appear at higher energies [1,2]. However, the core scaling capability is limited by the transverse multimode character that inevitably appears for large core fibers, reducing the beam quality of the laser source. Moreover, the onset of self-focusing at peak powers of the order of 4 MW is independent of the core size and represents a hard limit for the peak power that can be handled in silica fibers.

Another way to scale the energy is to use several amplifiers and coherently combine their outputs in free space to get the sum of the power while maintaining the beam quality and avoiding detrimental effects at high peak powers in the fiber. This technique is already used in CW [3] and studied in nanosecond regimes [START_REF] Cheung | Advanced Solid-State Photonics, OSA Technical Digest Series[END_REF], and is scalable in terms of achievable energy.

In another context, coherent addition of femtosecond optical pulses has been observed recently with two pulses without any spectral overlap, allowing the generation of single cycle pulses [START_REF] Krauss | [END_REF]. In this case, the motivation is not to scale the power but to obtain the shortest possible pulse.

Here we report on the coherent beam combining of two identical femtosecond Yb-doped fiber CPAs, using active control of the relative phase between the amplifiers by use of an integrated electro-optics phase modulator. Two 4 W femtosecond amplifiers are combined, and the output power is 7:2 W, resulting in a combining efficiency of 90%. The spatial, temporal, and spectral characteristics of a single amplifier are completely recovered in the combined beam, with output pulse duration of 325 fs. This combining experiment opens a novel way for energy/average power scaling of femtosecond amplifiers.

In CW coherent combining experiments, the relative phase between the beams is measured and controlled through a feedback loop. In femtosecond regime, since the optical bandwidth is larger, the optical phase must be matched over the whole bandwidth to efficiently combine the pulses. In other words, the phase, group delay, group-delay dispersion, and so on, should be matched. The tolerance of the system to mismatch of these parameters is inversely related to the pulse width. In what follows, i.e., for pulse durations above 100 fs, we show that a static adjustment of the group delay and groupvelocity dispersion, along with an active correction of the phase, is sufficient to efficiently combine the pulses. Indeed, the group delay and group-delay dispersion fluctuations between the pulses are sufficiently weak so that the final pulse width is not affected.

Our experimental setup, depicted in Fig. 1, begins with a femtosecond air-cooled bulk Yb 3þ : KYW oscillator delivering 260 fs pulses at 1030 nm and 35 MHz repetition rate with 1 W average power. It is followed by a stretcher that broadens the pulses to 150 ps. A half-wave plate with a polarizing beam splitter allows us to separate the seed into two arms with adjustable input power for each. The first one consists of a fiber-coupled LiNbO 3 integrated phase modulator followed by a 1:2 m Yb-doped doubleclad polarization-maintaining 30=150 μm large mode area (LMA) fiber amplifier. The second one contains a 2:40 m single-mode fiber, a free-space delay line, and another 1:2 m Yb-doped double-clad 30=150 μm LMA fiber amplifier. The aim of the single-mode fiber is to match the group-velocity dispersion of both arms, and it also acts as a coarse delay matching element. Both amplifiers are seeded with 100 mW of average power and are pumped to deliver identical powers. At this repetition rate, we operate in completely linear regime, so that self-phase modulation-induced degradation of the combining efficiency does not take place [6]. The outputs of the fiber amplifiers are collimated and overlapped on a 50=50 beam splitter. Provided that the pulses have the same phase, arrival time, temporal shape, polarization, and spatial profile, they can be recombined coherently. The constructive output of the beam splitter goes to the compressor that exhibits a 60% efficiency. The other one is followed by a photodiode with a bandwidth of 1 MHz that detects a signal to be minimized.

The relative phase between the pulses is controlled using the technique described in Ref. [3]. A small voltage modulation at 250 kHz is added to the control signal that drives the phase modulator, yielding a small modulation (0:1 rad) of the relative phase between the two arms. At the output of the photodiode, a lock-in amplifier allows the detection of the small induced amplitude modulation. The resulting error signal is directly proportional to the optical phase difference between the arms. A digital controller then maintains the phase difference close to 0 by adding an offset voltage to the small modulation voltage and manages 2π phase jumps to stay in the phase modulator input voltage range. Since the modulator bandwidth is >1 GHz, the phase jump is instantaneous for pulses at 35 MHz repetition rate. The integration time of the lock-in amplifier is set to 20 μs, meaning that the maximum phase noise frequency that can be corrected is 25 kHz. The controller parameters are set to minimize the RMS fluctuations. The overall feedback loop bandwidth in our case is around 1 kHz, determined by the controller integration. However, this method was shown to be compatible with larger bandwidths, as demonstrated in Ref. [START_REF] Shay | Advanced Solid-State Photonics, OSA Technical Digest Series[END_REF].

The photodiode signal versus time is shown in Fig. 2, both in free-running and locked phase regimes, over 100 s. The free-running signal evolves arbitrarily between the in-phase and out-of-phase states, which determines a contrast ðI max -I min Þ=ðI max þ I min Þ of 84%. This contrast is very sensitive to the spectral intensity of each pulse and depends more loosely on the group delay and group-delay dispersion discrepancies. In phase-locked operation, the signal shows residual intensity fluctua-tions that correspond to an RMS phase noise of λ=40 over the period of 100 s, assuming that the intensity noise is negligible. No readjustment of the delay line is necessary to keep the system locked over hours, indicating that the group-delay drift is small compared to the pulse duration.

Figure 3 shows the phase noise spectral density at the output of the recombination system, computed from the intensity measurement, for free-running and phaselocked operation. Most of the phase noise content is below 1 kHz, and our setup is therefore well adapted to correct this phase fluctuation bandwidth. In particular, the acoustic fluctuations above 100 Hz are compensated by the phase-locking electronics.

Each amplifier delivers 4 W with a 10 W pump power. The recombined power is measured to be 7:2 W before compression, which corresponds to a combining efficiency P out =ðP 1 þ P 2 Þ of 90% and a compressed power of 4:3 W. The output polarization state of each amplifier is measured to be linear to more than 99%. Therefore, polarization has little effect on the efficiency. To achieve this efficiency, the spatial overlap of the beams at the output of the beam splitter must be very carefully optimized. Any difference between the wavefronts, such as defocus and tilt, or in the beam profiles, such as ellipticity and beam position, results in a drastically reduced combining efficiency. The intensity profiles of each arm separately and of the combined beam are shown in Fig. 4. The individual profiles do not exhibit exactly the same shape, and the combined beam is intermediate. This slight spatial discrepancy contributes to the nonperfect combining efficiency.

Finally, we focus on the temporal and spectral properties of the pulses, for each channel independently and for the recombined output. The autocorrelations and spectra are plotted in Fig. 5. The autocorrelations exhibit a very good match with the sech 2 temporal shape, so that a 1.54 deconvolution factor is used to infer the pulse widths. The pulse durations for each isolated arm are 315 fs for the first one and 335 fs for the second one. Thus, for the same compression, the pulses do not have exactly the same duration, indicating a small group-delay dispersion discrepancy. The recombined pulse exhibits 325 fs duration, corresponding to the mean of the initial durations. However, a simple model can be used to show that the contribution of this temporal mismatch to the overall combining efficiency is negligible compared to spatial aspects.

The recombined spectrum is also very similar to the isolated ones. Again, the recombined spectrum is between both individual spectra, with a 4:3 nm FWHM. The time-bandwidth product of the recombined pulse is 0.39, close to the Fourier-transform-limited product of 0.315 for a sech 2 profile. This slightly degraded time-bandwidth product is the same for a single amplifier and is due to a nonperfect alignment of the stretcher and compressor, resulting in high-order spectral phase mismatch. The temporal quality after the combining process is therefore very well preserved.

To conclude, we demonstrate coherent beam combining of two femtosecond fiber CPAs using a single stretcher and compressor, resulting in an output power of 7:2 W before compression, with 90% efficiency. The RMS residual phase error is λ=40, and the final pulse width is 325 fs, with a good spatial and temporal quality. Theoretical considerations to evaluate the tolerance of the system to residual discrepancies between the arms in terms of phase, group delay, dispersion, and B integral indicate that these systems are realistically scalable. To our knowledge, only one similar experiment has been done independently by another research team and reported very recently, using an alternative method based on polarization for the phase control [START_REF] Seise | [END_REF] and for longer pulses. However, the use of the polarization state to generate the error signal implies that only two amplifiers can be combined by one system, and the scalability of this scheme might be challenging. In our case, although the results were obtained with two amplifiers using a beam splitter, the phase feedback loop is compatible with scalable combining geometries such as the tiled aperture technique. Considering that a seed of 10 mW is sufficient to overcome amplified spontaneous emission in our fibers, this setup is thus already applicable to a few tens of amplifiers. Higher scalability can be obtained with preampli-fiers that might be implemented in an all-fiber setup. These experiments represent a promising way towards future high-average and peak power laser systems.

The research described here has been supported by Triangle de la Physique contract 2009-026T.

Fig. 1 .

 1 Fig. 1. (Color online) Experimental setup of the femtosecond coherent combining experiment.

Fig. 2 .

 2 Fig. 2. (Color online) Photodiode signal versus time, over 100 s, in free run and phase-locked regime.

Fig. 3 .

 3 Fig. 3. (Color online) Phase noise density of the recombination, in free run (black) and locked regime (red).

Fig. 4 .

 4 Fig. 4. (Color online) Intensity profiles of (left), (center) the individual channel beams and of (right) the recombined beam at the output of the compressor.

Fig. 5 .

 5 Fig. 5. (Color online) (top) Autocorrelation traces and (bottom) optical spectra for each channel independently (blue and black) and for the recombined beam (red).

March 1, 2011 / Vol. 36, No. 5 / OPTICS LETTERS

0146-9592/11/050621-03$15.00/0 © 2011 Optical Society of America