Bradley Worley
email: bradley.worley@pasteur.fr

Florent Delhommel
email: florent.delhommel@pasteur.fr

Florence Cordier
email: florence.cordier@pasteur.fr

Thérèse E Malliavin
email: therese.malliavin@pasteur.fr

Benjamin Bardiaux

Nicolas Wolff
email: nicolas.wolff@pasteur.fr

Michael Nilges
email: michael.nilges@pasteur.fr

Carlile Lavor
email: clavor@ime.unicamp.br

Leo Liberti
email: liberti@lix.polytechnique.fr

Tuning interval Branch-and-Prune for protein structure determination

Keywords: Distance geometry, Protein structure, Nuclear Magnetic Resonance, Branch-and-Prune

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Within biology, it is well-established that the biochemical function of a molecule is strongly related to its three-dimensional structure. As a direct consequence, substantial effort is devoted within the field of structural biology to the problem of molecular structure determination: given the chemical composition and topology of a molecule, we seek its conformation in R 3 . For the specific problem of protein structure determination, the composition and topology of a protein molecule are completely specified by its amino acid sequence, 1 and we again seek its conformation(s) in three dimensions.

Proteins are highly flexible polymer chains of amino acids, and consequently may have multiple conformations for a given amino acid sequence. Therefore, additional geometric measurements are required in order to obtain a reasonable number of solutions to the protein structure determination problem. Nuclear Magnetic Resonance (NMR) experiments are frequently employed to obtain these measurements. In an NMR experiment, the interaction between the nucleus of each atom in a molecule and a strong magnetic field is measured using precisely timed pulses of radio-frequency radiation [START_REF] Levitt | Spin Dynamics: Basics of Nuclear Magnetic Resonance[END_REF]. The interaction energy of a nucleus as a fraction of the total magnetic field energy is known as its chemical shift. Chemical shifts are highly sensitive reporters of the local electronic environments of their nuclei, making them excellent tools for structure determination. For proteins, whose polymer backbones are formed by repeated N-C α -C units, chemical shifts are dominated by local backbone geometry (Fig. 1), in particular the backbone dihedral angles φ , ψ and ω. The φ and ψ angles are historically referred to as the Ramachandran angles [START_REF] Ramachandran | Stereochemistry of polypeptide chain configurations[END_REF]. For amino acid i of a protein, the φ i angle is defined by the atoms C (i-1) -N (i) -C (i) α -C (i) , the ψ i angle is defined by N (i) -C (i) α -C (i) -N (i+1) , and the ω i angle is defined by

C (i) α -C (i) -N (i+1) - C (i+1) α
. The ω i dihedral angle is usually fixed to 180 • due to the known near-planarity of the peptide bond [START_REF] Berkholz | Nonplanar peptide bonds are common and conserved but not biased toward active sites[END_REF]. When all ω dihedrals are fixed to 180 • , and the local ge- 1 Additional post-translational modifications may alter the chemical composition of a protein molecule, but their effects must be detected by further chemical analysis.

C (i-1)
O (i-1)

N (i) H (i) 1 C (i) α C (i) β H (i) α C (i) O (i) N (i+1) H (i+1) 1 C (i+1) α C (i+1) β H (i+1) α φ i ψ i ω i
Fig. 1 Depiction of a single amino acid unit i of a protein backbone, including any atoms of flanking units i -1 and i + 1 necessary for defining the backbone dihedral angles φ i , ψ i and ω i . While the repetition vertex orders introduced in sec. 2.3 do not include C β atoms, the NMR chemical shifts of these atoms are important indicators of the local dihedral angles (cf. 1).

Fig. 2 NMR spectrum obtained from the 1 H- 15 N heteronuclear single quantum coherence (HSQC) experiment for the HHD2 protein sample described in sec. [START_REF] Brünger | Computational challenges for macromolecular structure determination by X-ray crystallography and solution NMR-spectroscopy[END_REF]. Peaks in the spectrum labeled with a number correspond to pairs of bonded

H (i)
1 and N (i) atoms in the protein, where the number corresponds to the index of the amino acid unit (i). The protein studied here is a short sub-sequence of a much longer amino acid sequence, so the values of i indicate the position of the synthesized amino acids within that longer sequence.

ometry is assumed to be fixed [START_REF] Hinsen | A comparison of reduced coordinate sets for describing protein structure[END_REF], the conformation of a protein's backbone may be completely specified by supplying the Ramachandran angles (φ , ψ) for all of its amino acid units. These angles may be predicted from the chemical shifts of the N, C α , C β , C , H 1 and H α atoms using TALOS-N, an artificial neural network trained on a large database of previously determined NMR protein structures [START_REF] Shen | Protein backbone and sidechain torsion angles predicted from NMR chemical shifts using artificial neural networks[END_REF]. Multiple predictions from the neural network are used to estimate a mean and standard deviation for each backbone dihedral angle, resulting in fairly reliable Ramachandran angle intervals. However, when some chemical shift measurements are missing, the accuracy of the predicted dihedral angles can suffer. A set of multidimensional NMR experiments involving different nuclei [START_REF] Ikura | A novel approach for sequential assignment of 1 H, 13 C and 15 N spectra of proteins: heteronuclear triple-quantum resonance three-dimensional NMR spectroscopy. Application to calmodulin[END_REF] are required to determine the chemical shift of each atom, using a process known as sequential assignment (cf. §4). Fig. 2 illustrates one such experimental result in protein NMR, a two-dimensional spectrum that is used to correlate the chemical shifts of backbone

H (i)
1 and N (i) atoms. In addition to the local geometic information provided by chemical shifts, NMR may also be used to obtain global geometric information in the form of distances between pairs of atoms. In such experiments, the intensity of the NMR signal at the chemical shifts of two atoms, e.g. H (i) and H ( j) , is approximately proportional to the inverse sixth power of their distance, d(H (i) , H ( j) ) -6 . Because the signal intensity decays rapidly with increasing distance, this information is usually only available for atoms that are within 5-6 Å in the structure. Additional effects, including inevitable molecular motion during the measurement, can also result in a change of signal between nearby atoms. Finally, these signal intensity measurements are often perturbed by strong systematic measurement errors, so they are generally converted into interval distances using "rules of thumb" developed by NMR spectroscopists [START_REF] Wüthrich | NMR of Proteins and Nucleic Acids[END_REF].

Conventional NMR protein structure determination protocols combine known distances and angles between bonded atoms with predicted Ramachandran angles and interval distance measurements into a non-convex objective function that smoothly penalizes any deviations of the structure from the target geometry [START_REF] Brünger | Computational challenges for macromolecular structure determination by X-ray crystallography and solution NMR-spectroscopy[END_REF][START_REF] Mucherino | Finding low-energy homopolymer conformations by a discrete approach[END_REF]. As this objective contains many local optima, metaheuristics such as simulated annealing with multiple random initializations are employed. Nevertheless, there is no guarantee that the resulting structures are within the feasible set specified by the distance, angle and dihedral constraints.

The protein structure determination problem may be recast into a DISTANCE GE-OMETRY PROBLEM (DGP) by converting all aforementioned geometric constraints into distance constraints between pairs of atoms. Formally, the protein is represented as a graph G = (V, E, d), where V represents the atoms and E holds the atom pairs which have a known distance. The distances are either exact values or intervals, so there exists a partitioning of the edge set E = E D ∪ E I , where E D and E I hold exact and interval distances, respectively. We denote the set of positive real valued intervals as IR + . Exact distances in E D are given by known bond lengths, angles, and dihedrals in protein structures. Three-atom angles and four-atom dihedral angles that are known from protein chemistry shall be collected into the sets Θ and Ω , respectively. Given this information, we may formally re-introduce the INTERVAL DISCRETIZ-ABLE MOLECULAR DISTANCE GEOMETRY PROBLEM [START_REF] Liberti | Euclidean Distance Geometry and Applications[END_REF] as the following,

E D → R + E I → IR +
and a vertex ordering R = (v 1 , . . . , v n ) on V satisfying the following requirements:

-The subgraph of G induced by

V 0 = {v 1 , v 2 , v 3 } is a clique with all edges in E D -For all j ∈ R \V 0 we have 1. d j-1, j , d j-2, j ∈ E D 2. d j-3, j ∈ E -For all j ∈ R \ {v 1 , v 2 }
the strict triangular inequality is satisfied by d j-1, j , d j-2, j and d j-2, j-1 is there an embedding x :

V → R 3 such that x u -x v = d uv ∀{u, v} ∈ E D and x u -x v ∈ [d uv , d uv ] ∀{u, v} ∈ E I ?
In general, the vertex ordering R may visit each vertex in the graph one or more times. We define such a REPETITION ORDER [START_REF] Cassioli | Discretization vertex orders in distance geometry[END_REF] as a sequence R : N → V ∪ 0 with length |R| ∈ N (such that R i = 0 for all i > |R|) when it satisfies the requirements of the iDMDGP. A repetition order entry r i is referred to as repeated when there exist one or more elements v j with j < i -2, such that v j = v i . As shown in [6, §4], there exist iDMDGP instances whose vertex orders must contain repetitions.

The interval Branch-and-Prune (iBP) algorithm is a combinatorial method for finding solutions to the iDMDGP. The iDMDGP solution space is a search tree, within which each path from the root node to a leaf node represents a distinct protein conformation. The iBP algorithm, explained in detail previously [START_REF] Lavor | The interval Branch-and-Prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF][START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF][START_REF] Lavor | Clifford Algebra and the Discretizable Molecular Distance Geometry Problem[END_REF], recursively traverses the search tree in order to enumerate solutions to a given problem instance. If at any point in the search, the partial solution becomes infeasible with respect to the constraints, then the leaves of the current search tree node are "pruned" and iBP backtracks until a solution is identified. By searching the entire tree, iBP is capable of systematically exploring the feasible set of any given problem instance.

This work introduces a variant of iBP that uses an iterative tree traversal algorithm based on embedding equations derived from Clifford algebra [START_REF] Lavor | Clifford Algebra and the Discretizable Molecular Distance Geometry Problem[END_REF]. This iterative variant is ideal for iDMDGP instances that have highly repetitive vertex orders [START_REF] Cassioli | Discretization vertex orders in distance geometry[END_REF], as it requires no extra matrix computations for repeated atoms in the order. We then pair our iterative iBP algorithm with such a highly repetitive discretization vertex order that directly uses information on dihedral angles during tree traversal. When many backbone dihedral angles are known to high precision in a given problem instance, this method and vertex order produce significantly smaller search trees than prior methods. For the purposes of this work, the proposed algorithm and vertex order enable iBP to obtain solutions to problem instances containing strong backbone dihedral information. For such problem instances, it is challenging for iBP to obtain any solutions using previously introduced vertex orders that do not branch on backbone dihedrals. We show that our new implementation solves this family of instances efficiently, even when previous iBP implementations fail to obtain any solutions.

The remainder of this paper is organized as follows. In Section 2, we introduce the iterative iBP algorithm with its embedding equations, along with new repetition vertex orders for protein structures. In Section 3, we describe features of the new iterative iBP implementation that simplify its use in routine applications while retaining the remarkable flexibility and generality of iBP. In Section 4 we present some computational results, specifically focusing on benchmarking the ability of iBP to find a single solution. Finally, we conclude the paper in Section 5.

Algorithmic considerations

Iterative embedding relations

iBP is most naturally defined as a recursive tree traversal algorithm [START_REF] Lavor | The interval Branch-and-Prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF][START_REF] Costa | Discretization orders for protein side chains[END_REF] that uses recursively defined affine transformation matrices [START_REF] Thompson | Calculation of Cartesian Coordinates and Their Derivatives from Internal Molecular Coordinates[END_REF] for computing embedded coordinates. However, recent use of Clifford algebra has yielded embedding equations that are non-recursive [START_REF] Lavor | Clifford Algebra and the Discretizable Molecular Distance Geometry Problem[END_REF], the key results of which are recalled here.

For any vertex j in a repetition order R, we seek the embedded coordinate x j ∈ R 3 , given distances to the three preceding vertices j -1, j -2 and j -3. From the Fig. 3 Embedding relations for computing the position of a vertex x 4 from its predecessors (x 1 , x 2 , x 3 ), via the parameters (δ 4 , θ 4 , τ 4 , σ 4 ).

properties of the iDMDGP, the distances d j-1, j , d j-2, j and d j-3, j are known, where d j, j-3 is potentially an interval. The remaining distances in the clique formed by the four vertices ( j -3, j -2, j -1, j) are calculable from the coordinates x j-1 , x j-2 and x j-3 . We shall introduce the quantities δ j , θ j , τ j and σ j for embedding vertex j, where δ j d j-1, j (cf. Fig. 3). The angle θ j is uniquely defined, and is obtained from the cosine law using the relevant distances,

θ j cos -1 d 2 j-1, j + d 2 j-2, j-1 -d 2 j-2, j 2d j-1, j d j-2, j-1 (1)
From the requirement that the vertex order satisfies the strict triangular inequality, θ j is strictly within (0, π), ensuring the vertices ( j -3, j -2, j -1) are non-collinear. The variable τ j cos ω j is related to the dihedral angle ω j formed by the vertices ( j -3, j -2, j -1, j). When ω j is known from either protein chemistry or measurement, such that ω j ∈ Ω , we may directly compute τ j , as well as the sign variable σ j ∈ {-1, +1}:

σ j sin ω j | sin ω j | if sin ω j = 0 1 if sin ω j = 0
In cases where ω j / ∈ Ω , the iDMDGP nevertheless guarantees that the distance d j-3, j is available, and we may compute τ j from the cosine law for a trihedron [START_REF] Lavor | Clifford Algebra and the Discretizable Molecular Distance Geometry Problem[END_REF]:

τ j = 2d 2 j-2, j-1 d 2 j-3, j-2 + d 2 j-2, j -d 2 j-3, j -d j-3, j-2, j-1 d j-2, j-1, j 4d 2 j-3, j-2 d 2 j-2, j-1 -d 2 j-3, j-2, j-1 4d 2 j-2, j-1 d 2 j-2, j -d 2 j-2, j-1, j (2)
where

d j-3, j-2, j-1 d 2 j-3, j-2 + d 2 j-2, j-1 -d 2 j-3, j-1 d j-2, j-1, j d 2 j-2, j-1 + d 2 j-2, j -d 2 j, j-1
Given δ j , θ j , τ j , and σ j , the embedded coordinates of vertex j are given by the following equation:

x j = p 1 + τ j p 2 + σ j 1 -τ 2
j p 3 where p 1 , p 2 , p 3 ∈ R 3 depend only on x j-1 , x j-2 , x j-3 , δ j and θ j ,

p 1 = - δ j r 12 cos(θ j ) - r 12 δ j x j-1 -cos(θ j ) x j-2 p 2 = - δ j r 12 sin(θ j ) r 12 × r 23 r 12 2 r 23 -(r 12 • r 23 ) r 12 p 3 = - δ j r 12 sin(θ j ) r 12 × r 23 r 12 (r 12 × r 23 )
and we have introduced r 12 , r 23 ∈ R 3 for notational simplicity,

r 12 = x j-1 -x j-2 r 23 = x j-2 -x j-3
By explicitly parameterizing each embedding with the sign σ j , we obtain an unambiguous solution for the coordinate x j . An iBP search tree that uses these equations therefore stores a set of values (δ , θ , τ, σ , x) at each of its nodes. Thus, whenever ω j -and consequently σ j and τ j -is known for vertex j, this yields a single branch at that level in the tree.

In the present iBP implementation [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF] that uses recursively defined affine transformation matrices, iBP must compute and store a matrix for each repeated vertex in a repetition order. To compute the matrix of a repeated vertex, the sign parameter σ j must be determined by systematic search during branching: the value of σ j that yields an embedded coordinate nearest to the originally determined coordinate is used to build the transformation matrix. For vertex orders containing a large degree of repetition, this introduces unnecessary computations. Using the above equations, embedding a new vertex x j requires only (δ j , θ j , ω j , σ j ) and the embedded coordinates of the last three vertices in the order, so repeating a vertex in a discretization order introduces a copy of the original vertex coordinate, but requires no additional computations.

Like the iBP implementation using affine transformation matrices, our implementation requires the calculation of τ j from distances by eq. ( 2) when ω j / ∈ Ω . While this can lead to numerical instability in large instances with imperfect distance data, the fact that our implementation directly computes τ j when ω j ∈ Ω effectively mitigates this instability when dihedral constraints are available.

Thus far, the discussion of the Clifford algebraic embedding equations has focused on situations when either ω j or d j, j-3 are known exactly, resulting in a single value of τ j and either one or two values of σ j , respectively. When d j-3, j is an interval distance, such that { j -3, j} ∈ E I , we shall denote its lower and upper bounds as d j-3, j and d j-3, j , respectively. Similarly, when ω j is an interval dihedral, i.e. ω j ∈ Ω I , we denote its bounds as ω j and ω j . In a process that is described in more detail below, these intervals are then discretized through linear interpolation in order to obtain values that may be used in the embedding equations.

Algorithm 1 The advance method 1: Input: Index I and tree size N in N |R| , incremented element position j 2: Output: Least modified element position j * 3: I k ← 1 ∀k ∈ { j + 1, . . . , |R|} 4: for all k ∈ { j, j -1, . . . , 1} do 5:

I k ← I k + 1 6: j * ← k 7: if I k > N k then 8: I k ← 1 9:
else 10:

return j * 11:

end if 12: end for 13: return j *

Iterative tree traversal

In order to obtain solutions to the iDMDGP, iBP discretizes the search space by breaking each interval into a finite set of values, resulting in a search tree. We introduce the multi-index N ∈ N |R| , referred to as the tree size, which holds the number of branches at each level of this tree. Herein R is the repetition order associated with the iDMDGP instance. The number of branches at level j of the tree is given by N j , and depends on the geometric information available for vertex j. When j is a repetition of a previously embedded vertex, we have N j = 1. When d j-3, j is an exact distance, there are two possible embeddings of vertex j, so N j = 2. Finally, when d j-3, j is an interval distance, N j = 2B where B is a user-specified discretization factor. In cases where ω j is known-as opposed to knowing only cos ω j via d j-3, j -the branch count N j is reduced by a factor of two, as σ j takes a single value.

By default, interval discretization in iBP is done uniformly, with each of the B discrete points equally spaced upon its arc formed by the three-sphere intersection sub-problem. Due to discretization, iBP is a heuristic method. For any iDMDGP instance having a non-empty feasible set, there is no guarantee that iBP will select discretization points which satisfy the original iDMDGP [START_REF] Goncalves | Recent advances on the interval distance geometry problem[END_REF]. However, when sufficiently large values of B are used in concert with a small tolerance during distance feasibility pruning, the chances of obtaining solutions is increased. Thus, in the absence of alternative methods for discretizing the interval d j and ω j values, we defer to uniform subdivision.

In the following, we shall employ multi-index notation to describe the iterative tree traversal routine, with multi-indices denoted by capitalized boldface letters. Given a tree size N ∈ N |R| , we introduce the partial order operator for all |R|-indices I ∈ N |R| as follows:

I ≤ N → true if I j ≤ N j ∀ j ∈ {1, . . . , |R|} false otherwise
An index I describes a valid-though potentially infeasible-path through the tree if it satisfies 1 ≤ I ≤ N, where 1 is a vector of ones.

In the iterative formulation of iBP described in this paper, paths through the search tree are enumerated by advancing an index I (Algorithm 1). Advancing an index I at its last element I |R| yields a sequence of indices-and thus paths through a tree-that implicitly effects a depth-first search in that tree. On the other hand, incrementing an index I at some element I j , with j < |R|, is equivalent to pruning the tree at that level, as it skips all indices that would have been produced by incrementing at |R|. The term index is used to describe I due to the similarity it shares to the indices of multidimensional arrays of shape N. In this analogy, the leaves of the search tree are each given one element in a multidimensional array, which has |R| dimensions and N j elements along dimension j. Alternatively, the act of advancing an index may be considered equivalent to a set of combined operations on a stack S.

Proposition 1 (Stack equivalence) An index I ∈ N |R| , combined with a level j ≤ |R|, represents a stack S, such that S = (I 1 , I 2 , I 3 , . . . , I j ). By introducing the following equivalence between index operations and stack operations,

(i.) I k ← z, k > j: no operation. (ii.) I k ← z, k ≤ j: (a.) pop S until |S| = k -1, (b.) push z onto S. (iii.) j ← k, k > j: push 1 onto S until |S| = k.
we see that the advance method maintains a stack within the first j elements of its index I.

Proof (Stack equivalence) (i.) As S only contains the first j elements of I, the modification of an index element I k for k > j does not modify S. (ii.) This applies to lines 5 and 8 of Algorithm 1, and ensures that S k = I k for k ≤ j.

Note that this changes the size of the stack, which is reflected in the algorithm by the least modified index j * (iii.) This ensures that the initialization ( j ← 1) and the act of moving to the next tree level ( j ← j + 1) are equivalent to pushing the next available tree node onto the stack S.

Using these indices, we may construct an iterative variant of iBP that implicitly traverses and prunes a search tree by advancing an index (Algorithm 2). The cal-cAngle and calcTorsion methods in Algorithm 2 calculate the bond angle θ i and the cosine of the dihedral angle, τ i , using equations 1 and 2 above, using the provided vertex coordinates to compute the required distances.

Ramachandran-defined vertex orders

In order to leverage known backbone dihedral angles when enumerating solutions to the iDMDGP, a new set of vertex orders was introduced. We define an initial order Algorithm 2 The iterative iBP algorithm 1: Input: iDMDGP instance with graph G, repetition order R, and tree size N. 2: Output: Set of solutions X . 3: Initialize I ← 1, j = 1, and X ← / 0 4: while I ≤ N do 5:

while j ≤ |R| do 6:

if r j is a repetition then 7:

Copy x j from the original embedding x i {where i < j s. t. r i = r j } 8: else 9:

δ ← d j-1, j 10:

θ ← calcAngle(x j-1 , x j-2 , d j-1, j , d j-2, j ) 11:

if ∃ω j ∈ Ω then

12: ω ← I j N j (ω j -ω j ) + ω j ) 13: σ ← sin ω /| sin ω | {Or 1, if sin ω = 0} 14: τ ← cos ω 15: else 16: σ ← 2(I j mod 2) -1 17: d ← I j /2 N j /2 (d j-3, j -d j-3, j ) + d j-3, j 18: τ ← calcTorsion(x j-1 , x j-2 , x j-3 , d j-1, j , d j-2, j , d ) 19:
end if 20:

Embed x j given x j-1 , x j-2 , x j-3 , δ , θ , τ, σ 21:

if vertex x j is infeasible then 22:

(I, j) ← advance(I, N, j) 23:

break {Jump to outermost while (line 4).} 24:

end if 25:

end if 26:

if j = |R| then 27:

X ← X ∪ x {The current embedding x is a solution} 28:

end if 29: j ← j + 1 30:

end while 31:

(I, j) ← advance(I, N, |R|) 32: end while R 1 , an inner order R i and a final order R n for the first, inner, and last amino acid units of protein graphs, respectively. The vertex orders are as follows:

R 1 = N (1) , H (1) 1 , H (1) 2 , C (1)
α , N (1) , H

α , C (1) α , C (1) R i = N (i) , O (i-1) , C (i-1) α , C (i-1) , N (i) , C (i) α , C (i) , N (i+1) , C (i-1) , N (i) , C (i) α , H (i) 1 , N (i) , C (i) α , C (i) , H (i) α , C (i) , C (i) α R n = N (n) , O (n-1) , C (n-1) α , C (n-1) , N (n) , C (n) α , C (n) , C (n-1) , N (n) , C (n) α , H (n) 1 , N (n) , C (n) α , C (n) , H (n) α , C (n) , C (n) α , O (n) 1 , C (n) , O (n) 2 (1)
where i denotes the amino acid unit index and n denotes the total number of amino acid units within this context. These orders make extensive use of repeated vertices in order to achieve direct branching on the backbone dihedrals φ , ψ and ω.

Fig. 4 Reorder-dependent differences in sampling of Ramachandran space, using identical branching factors (B = 64). Black crosses and red dots represent conformations sampled with and without application of dihedral restraints to (φ , ψ), respectively. a. Old repetition orders. b. New Ramachandran-defined repetition orders. Fig. 4 illustrates how the new orders more directly and uniformly sample Ramachandran space than previously published vertex orders [START_REF] Costa | Discretization orders for protein side chains[END_REF]. An iBP calculation was performed, using both vertex orders, to find solutions to the iDMDGP instance of a backbone-only tetrapeptide, with the direct distance feasibility (DDF, [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF]) and van der Waals (VDW) pruning devices enabled. Because the previous orders do not branch on the φ or ω backbone dihedral angles of proteins, they do not enable iBP to cover this space in a regular way, and often sample nearby points in (φ , ψ)-space more than once during tree traversal (Fig. 4a). An even stronger difference in behaviour is observed when constraints are placed on an instance's φ and ψ angles (Fig. 4, red dots). Using previously published orders, iBP treats the dihedral constraints as extra pruning edges, thus maintaining the same irregular sampling pattern and pruning any partial conformations that fall outside the constraints. Using the new orders that directly branch on protein backbone dihedral angles, iBP adaptively changes its sampling pattern (Fig. 4b, red dots) to better cover the region specified by the constraints. The red rectangle in the second panel of Fig. 4 is indeed a regular grid of sampled conformations. The central regions of (φ , ψ)-space were not sampled by either vertex order, as they were pruned by the VDW device.

3 Implementation-specific technical innovations

CHARMM-syntax force fields

The basic input to protein structure determination procedures is the amino acid sequence of the protein of interest. From the sequence, software packages draw information from "force field" libraries in order to construct the topology (graph structure) and parameters (graph edge weighting) of the target molecule. Routinely employed software packages [START_REF] Weiner | AMBER: Assisted model building with energy refinement. A general program for modeling molecules and their interactions[END_REF][START_REF] Brooks | CHARMM: A program for macromolecular energy, minimization, and dynamics calculations[END_REF][START_REF] Brünger | Crystallography and NMR System: A new software suite for macromolecular structure determination[END_REF][START_REF] Schwieters | The Xplor-NIH NMR Molecular Structure Determination Package[END_REF] employ a common force field syntax derived from the Chemistry at HARvard Molecular Mechanics (CHARMM) software package [START_REF] Brooks | CHARMM: A program for macromolecular energy, minimization, and dynamics calculations[END_REF].

To ensure extensibility, flexibility and interoperability with these packages, our implementation of iBP also uses CHARMM-syntax force fields to construct iDMDGP instances. In addition, iBP accepts a superset of the CHARMM parameter file syntax that allows users to specify interval distances, angles and dihedrals within force fields. The force field used by iBP was derived from the protein based force field of the Crystallography and NMR System (CNS, [START_REF] Brünger | Crystallography and NMR System: A new software suite for macromolecular structure determination[END_REF]) with minor additions. For example, the topology information of a protein backbone is defined in iBP using the following notation: Vertex orders in iBP are also specified using a custom syntax that retains the flexibility of CHARMM-style force fields. Each amino acid unit defined in the force field topology files is given a corresponding order. As an example, the Ramachandrandefined vertex order of a protein backbone corresponding to the above topology entry is specified as follows:

! BBI: interior (1<i<n)
reorder BBI N, -O, -CA, -C, N, CA, C, +N, -C, N, CA, H1, N, CA, C, HA, C, CA end where -O denotes atom O of the previous amino acid unit, and +N denotes atom N of the next. At runtime, iBP constructs the iDMDGP instance of a protein from its sequence using this combined topology, parameter, and vertex order information. The validity of the resulting problem instance is then checked, and tree traversal is initiated via Algorithm 2. Therefore, it is now straightforward to use iBP within routine structure determination efforts. In addition, this flexibility enables iBP to handle structure determination involving post-translational modifications and non-natural amino acids, provided a repetition order can be crafted.

Pruning and timing statistics

Our implementation of iBP admits the inclusion of new pruning devices via a system of callback functions. On tree initialization, each pruning device registers its callbacks at each level of the tree. During tree traversal, iBP executes each registered Fig. 5 The first solution returned by iBP for the HHD2 iDMDGP instance, using only dihedral constraints predicted from NMR chemical shifts. callback function to determine the feasibility of newly embedded vertices. Finally, when traversal is completed or terminated prematurely by the user, each pruning device outputs vertex-specific pruning results. This is useful for identifying any distance, angle or torsion constraints that are geometrically inconsistent, for example.

Furthermore, the index-based traversal algorithm enables iBP to roughly estimate its runtime. Given an index I in a tree of size N, we define the "width" of the sub-tree from 1 to I as the number of leaves in that sub-tree,

w(I, N) = |R| ∑ j=1 (I j -1) |R| ∏ k= j+1 N k
Therefore, the width of the entire tree is w(N, N)+1, and the width of the yet-untraversed portion of the tree is given by w(N -I + 1, N). Assuming the rate of tree traversaldefined as the number of leaf nodes traversed per unit time-is relatively constant on average, we estimate the remaining runtime of an iBP tree traversal as follows:

t r = w(N -I + 1, N) w(I, N) t e
where t r and t e are the remaining and elapsed runtimes, respectively. Using these equations, iBP may output tree size and periodic timing information to the user. This assumption of a constant traversal-and thus pruning-rate of iBP is ensured by measuring w(N -I + 1, N) at sufficiently spaced time intervals.

Fig. 6 Comparison of the first solution returned by iBP for the HHD2 iDMDGP instance (Fig. 5, shown here in dark grey) with the ten lowest-energy structures from ARIA/CNS (colored, thin traces), using the larger initial structure set and slower annealing rate.

4 Computational results

Experimental data

The vertex orders in this work are ideally suited to pairing with backbone dihedral angles, such as those predicted from measured NMR chemical shifts (e.g. by TALOS-N, [START_REF] Shen | Protein backbone and sidechain torsion angles predicted from NMR chemical shifts using artificial neural networks[END_REF]). To illustrate this, the structure of an HHD2 protein domain (residues G420-R499 of Whirlin isoform-4) was solved by constraining all backbone dihedral angles to the intervals predicted by TALOS-N (Fig. 5). NMR chemical shifts were measured from a uniformly 15 N, 13 C-labeled protein sample that was prepared at a concentration of 300 µM in 250 µL of a buffer solution (150 mM NaCl, 50 mM Tris-HCl, 5% D2O, pH 7.5). All data were collected at 25 • C on a Bruker Avance III 900 MHz spectrometer with a three-channel cryogenically cooled probe. The following experiments were performed, yielding chemical shifts for the following types of atomic nuclei:

-BT-HSQC [START_REF] Lescop | A set of BEST triple-resonance experiments for time-optimized protein resonance assignment[END_REF]: H 1 , N. (cf. Fig. 2) -HNHA [START_REF] Vuister | Quantitative J correlation: a new approach for measuring homonuclear threebond J HNHα coupling constants in 15 N-enriched proteins[END_REF]:

H 1 , N, H α .
-BT-HNCO [START_REF] Lescop | A set of BEST triple-resonance experiments for time-optimized protein resonance assignment[END_REF], BT-HNCO+ [START_REF] Gil-Caballero | HNCA+, HNCO+, and HNCACB+ experiments: improved performance by simultaneous detection of orthogonal coherence transfer pathways[END_REF]:

H 1 , N, C .
-BT-HNCACB, BT-HNCOCACB [START_REF] Lescop | A set of BEST triple-resonance experiments for time-optimized protein resonance assignment[END_REF]:

H 1 , N, C α , C β .
From these experiments, about 99% of the expected backbone chemical shifts were assigned (99% of H 1 , N, C α , C β and C ; 91% of H α ). The assigned chemical shifts, along with the amino acid sequence of HHD2, were then used to obtain predicted intervals for φ and ψ backbone dihedral angles using TALOS-N. The ω backbone dihedral angles were fixed to 180 • , following standard practice in the structural bioinformatics field. The resulting iDMDGP instance contained 464 vertices, and had a vertex order length |R| = 1378.

Structure calculations

To obtain the structure illustrated in Fig. 5, the iterative iBP variant described in Algorithm 2 was run using the Ramachandran-defined vertex orders described in Section 2.3. A branching factor of B = 16 and a branch epsilon of ε = 0.01 Å were used, resulting in an effective branching factor B eff, j that varies at each level of the tree according to the following equation:

B eff, j min B, d j, j-3 -d j, j-3 ε
The direct distance feasibility (DDF, [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF]) and van der Waals (VDW) pruning devices were enabled during iBP tree traversal, as well as a new pruning device that tested for direct dihedral angle feasibility. In short, dihedral feasibility pruning ensures that all quartets of atoms with corresponding dihedral angles in Ω are consistent with their respective dihedral constraints. This dihedral feasibility pruning device effectively generalizes previously developed iBP pruning devices related to proteins, including the α-helix and chirality pruning devices introduced in [START_REF] Cassioli | An algorithm to enumerate all possible protein conformations verifying a set of distance constraints[END_REF].

To provide a basis for comparison, conventional NMR structure calculations were performed by molecular dynamics simulated annealing (MDSA) using ARIA/CNS [START_REF] Brünger | Crystallography and NMR System: A new software suite for macromolecular structure determination[END_REF][START_REF] Rieping | ARIA2: automated NOE assignment and data integration in NMR structure calculation[END_REF], which was recently evaluated as one of the most effective software tools for NMR structure determination [START_REF] Rosato | The second round of Critical Assessment of Automated Structure Determination of Proteins by NMR: CASD-NMR-2013[END_REF]. In MDSA, the motion of each protein structure is simulated by numerically solving Newton's equations of motion from various initial velocities, during which the temperature of the system is reduced from ∼ 10 4 K to ∼ 50 K. The structures from simulated annealing are then subjected to a round of local descent. For both the dynamics and descent stages, the objective function is an approximation of the molecule's potential energy, which includes terms for the known local geometry and the NMR-derived geometric constraints [START_REF] Mucherino | Finding low-energy homopolymer conformations by a discrete approach[END_REF].

Within ARIA/CNS, default parameters were used to randomly generate two sets of 20 and 100 initial structures, based solely on the predicted (φ , ψ) angles from TALOS-N. The smaller set of structures was subjected to local descent using a short MDSA calculation, whereas the larger set was subjected to a longer MDSA calculation with a ten-fold slower annealing rate.

While complete traversal of the search tree, which contained 10 147 leaves, was estimated to require 10 140 minutes, the solution illustrated in Fig. 5 was obtained by iBP within a few seconds, thanks to the strong constraints supplied from TALOS-N. In contrast, average times of the short and long MDSA computations were 29 seconds and 186 seconds, respectively. However, none of the structures produced by short MDSA were in the dihedral angle feasible set, and only 71% of structures from long MDSA were feasible. A comparison of the overall folds produced by iBP and the long MDSA run of ARIA/CNS is given in Fig. 6; while the secondary structures and general fold are the same among all structures, only the iBP structure is guaranteed to be feasible in the iDMDGP.

Sensitivity to uncertainty

In general graphs, iBP exhibits exponential worst-case complexity [START_REF] Liberti | The Discretizable Molecular Distance Geometry Problem seems Easier on Proteins[END_REF]. Fortunately, the fact that pruning edges of similar length are distributed sufficiently uniformly over the vertex order ensures that branches do not occur frequently, resulting in tractable instances. However, when intervals are employed for branching in iBP, discretization can potentially spoil the favorable properties of the method [START_REF] Goncalves | Recent advances on the interval distance geometry problem[END_REF]. Lower discretization factors B increase the probability that no path through the search tree remains feasible after applying the discretization. Conversely, increasing the discretization factor grows the search tree exponentially, resulting in dramatically longer runtimes.

To analyze the sensitivity of the proposed iBP algorithm and vertex order to interval uncertainty, a set of experiments was performed using the first solution from the TALOS instance (the target structure) as a basis. First, dihedral restraints were obtained by computing the backbone φ and ψ dihedral angles from the target structure, and varying degrees of uncertainty (∆ ) were added to the dihedrals of residues 11-13, 30-34, 46-49, and 61-65. These residues were found in "loop" regions of the target structure, which are generally more flexible in protein structures. In addition, long-range distance restraints were obtained by collecting the following set of distances from the target structure:

D = { x i -x j < d max ; |i -j| ≥ 5}
where d max specifies the maximum admissible distance in the set. The final distance restraints were obtained by forming intervals 0.5 Å wide around each distance in D. We shall refer to the size of the set D as m dist in all further discussions. For the problem instances resulting from each pair of values (∆ , d max ), iBP was run using identical parameters to those used to solve the original TALOS instance, with the exception that the DDF tolerance was expanded to 0.1 Å. The time required for iBP to obtain a single solution, subject to a time limit of 5 hours, was recorded as t 1 for each instance, and the root mean square deviation of each obtained solution to the target structure was also recorded (Fig. 7). The complete set of results of this analysis is given in Tables 1, 2, and 3. This analysis illustrates several general behaviors of iBP for solving protein instances. When dihedrals are specified with high precision (Table 1), the conformational space is small, and iBP rapidly obtains a solution. As ∆ is increased, iBP still rapidly obtains a solution when only a few long-range distances are supplied. However, their deviations from the target structure increase with ∆ , and supplying more distances decreases the deviations at the expense of increased computation time. Finally, the inclusion of a large number of distance restraints tends to increase runtime past the 5 hour limit, especially when the conformational space is large due to large ∆ . While one may expect runtimes to decrease as more distances are added, due to the graph becoming more complete [START_REF] Liberti | The Discretizable Molecular Distance Geometry Problem seems Easier on Proteins[END_REF], the opposite effect is observed. This is a result of the negative interplay between interval discretization and distance pruning in iBP, and might be alleviated using error-tolerant pruning devices such as mean distance error (MDE) [START_REF] Goncalves | Recent advances on the interval distance geometry problem[END_REF].

Finally, it is important to note that none of the artificial instances produced from the original TALOS structure generated solutions when running iBP with previously defined vertex orders [START_REF] Costa | Discretization orders for protein side chains[END_REF], due to the irregular sampling patterns of Ramachandran space that those orders produce.

Conclusions

This paper introduces a new implementation of the interval Branch-and-Prune algorithm for molecular structure determination [START_REF] Lavor | The interval Branch-and-Prune algorithm for the discretizable molecular distance geometry problem with inexact distances[END_REF], and describes several modifications that provide enhanced performance for the specific task of protein structure elucidation. Comparison against an existing state of the art method in the field demonstrates the advantages of iBP in practical structure determination problems. Our new iBP implementation is completely open source, and is available under the MIT license at http://github.com/geekysuavo/ibp-ng.

  INTERVAL DMDGP: given a simple weighted undirected graph G = (V, E, d) where E = E D ∪ E I and, d :

Fig. 7

 7 Fig. 7 Summary of root mean square deviations (rmsd) between the first solution produced by iBP from the TALOS instance and the target structure, as a function of dihedral uncertainty (∆ ) and distance threshold d min . Full results are listed in Tables 1, 2, and 3.

Table 1

 1 Results of solving the TALOS instance under low dihedral angle uncertainties (∆= 0 • , 1 • , 2 • ).

	∆ ( • ) d max ( Å) m dist	t 1 (s) rmsd ( Å)
	0	5	0	1.7	0.014
	0	6	1	1.7	0.014
	0	7	6	1.7	0.014
	0	8	29	1.7	0.014
	0	9	81	1.7	0.014
	0	10	163	1.7	0.014
	0	12	314	1.7	0.014
	0	15	625	1.7	0.014
	0	20	1150 1.7	0.014
	1	5	0	1.7	0.014
	1	6	1	1.7	0.014
	1	7	6	1.7	0.014
	1	8	29	1.6	0.014
	1	9	81	1.7	0.014
	1	10	163	1.7	0.014
	1	12	314	1.7	0.014
	1	15	625	1.7	0.014
	1	20	1150 1.7	0.014
	2	5	0	1.8	0.014
	2	6	1	1.7	0.014
	2	7	6	1.7	0.014
	2	8	29	1.7	0.014
	2	9	81	1.7	0.014
	2	10	163	1.7	0.014
	2	12	314	1.7	0.014
	2	15	625	1.7	0.014
	2	20	1150 1.7	0.014

Table 2

 2 Results of solving the TALOS instance under moderate dihedral angle uncertainties (∆ = 5 • , 10 • , 15 • ).

	∆ ( • ) d max ( Å) m dist	t 1 (s)	rmsd ( Å)
	5	5	0	1.7	0.522
	5	6	1	1.7	0.522
	5	7	6	1.7	0.522
	5	8	29	1.6	0.522
	5	9	81	1.7	0.555
	5	10	163	2.1	0.360
	5	12	314	96.4	0.317
	5	15	625	-	-
	5	20	1150 -	-
	10	5	0	1.3	1.130
	10	6	1	1.2	1.130
	10	7	6	1.2	1.130
	10	8	29	1.7	0.720
	10	9	81	79.1	0.611
	10	10	163	638.0	0.529
	10	12	314	-	-
	10	15	625	-	-
	10	20	1150 1378.0	0.289
	15	5	0	0.8	1.732
	15	6	1	0.8	1.732
	15	7	6	0.8	1.732
	15	8	29	1.2	1.557
	15	9	81	1.3	1.381
	15	10	163	18.9	1.398
	15	12	314	1400.8	1.169
	15	15	625	-	-
	15	20	1150 -	-

Table 3

 3 Results of solving the TALOS instance under high dihedral angle uncertainties (∆ = 20 • , 25 • , 30 • ).

	∆ ( • )	d max ( Å) m dist	t 1 (s)	rmsd ( Å)
	20	5	0	1.4	2.028
	20	6	1	1.3	2.028
	20	7	6	1.5	2.028
	20	8	29	2.8	1.334
	20	9	81	4.4	2.130
	20	10	163	8.2	1.632
	20	12	314	1004.2	1.389
	20	15	625	-	-
	20	20	1150 -	-
	25	5	0	1.4	3.222
	25	6	1	1.3	3.566
	25	7	6	0.8	3.566
	25	8	29	3.2	2.742
	25	9	81	3.3	2.317
	25	10	163	-	-
	25	12	314	-	-
	25	15	625	-	-
	25	20	1150 -	-
	30	5	0	1.6	3.589
	30	6	1	2.1	3.589
	30	7	6	1.9	3.589
	30	8	29	17.1	2.472
	30	9	81	20.7	2.266
	30	10	163	106.5	1.882
	30	12	314	-	-
	30	15	625	-	-
	30	20	1150 -	-

Acknowledgements Bradley Worley thanks the Pasteur Foundation for postdoctoral support. The authors wish to thank Institut Pasteur, Ecole Polytechnique, CNRS (TGIR-RMN-THC FR3050), and the French research agency (ANR-10-BINF-03-08 "Bip:Bip") for financial support. Carlile Lavor wishes to acknowledge the Brazilian research agencies CNPq and FAPESP for support.