<figure><img src="image_1.png" /></figure>

HAL
open science

# L'extension de la géométrie des distances

Simon B Hengeveld

► To cite this version:

Simon B Hengeveld. L'extension de la géométrie des distances : Etude de trois applications, en particulier en biologie structurale. Bio-informatique [q-bio.QM]. Université de Rennes, 2024. Français.

NNT: tel-04607576

HAL Id: tel-04607576

https://theses.hal.science/tel-04607576

Submitted on 10 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

Public Domain
THÈSE DE DOCTORAT DE

L'UNIVERSITÉ DE RENNES

ÉCOLE DOCTORALE N° 601
Mathématiques, Télécommunications, Informatique, Signal, Systèmes, Électronique
Spécialité : Informatique

Par
**Simon HENGEVELD**

Extending Distance Geometry methods : not only distances

A study of three applications, with a focus on structural biology

Thèse présentée et soutenue à Rennes, le 8 février 2024

Unité de recherche : IRISA

Rapporteurs avant soutenance :

Adrien GOËFFON Professeur, LERIA/Université d'Angers
Juan CORTÉS Directeur de recherches, LAAS-CNRS à Toulouse

Composition du Jury :

Président : Thérèse MALLIAVIN Directeur de Recherches (HDR), LORIA-CNRS à Nancy

Examinateurs : Adrien GOËFFON Professeur des universités (HDR), LERIA/Université d'Angers

Juan CORTÉS Directeur de Recherches (HDR), LAAS-CNRS à Toulouse

Douglas GONÇALVES Professeur, Universidade Federal de Santa Catarina

Silvia BOTTINI Chaire de professeur junior (CPJ), INRAE à Nice

Thérèse MALLIAVIN Directeur de Recherches (HDR), LORIA-CNRS à Nancy

Dir. de thèse : Antonio MUCHERINO Maître de conférence (HDR), IRISA/Université de Rennes
<figure><img src="image_2.png" /></figure>
ACKNOWLEDGEMENT

First of all, I want to give my thanks to my supervisor, Antonio Mucherino, for his
guidance and support throughout the PhD thesis. I would also like to thank the review-
ers (Adrien Goëffon and Juan Cortés) and the committee members (Douglas Gonçalves,
Silvia Bottini and Thérèse Malliavin) for their time, involvement and fruitful feedback.
Furthermore, I want to acknowledge the different collaborators that I have had the plea-
sure of working with on the various projects carried out as part of my PhD: Thérèse
Malliavin, Jung-Hsin Lin, Leo Liberti, Wagner da Rocha, Frank Plastria, David Pelta,
Nara Rubiano da Silva, Douglas Gonçalves and Paulo H. Souto Ribeiro.

I am extremely grateful to my girlfriend, Sofia, who always supported me through
both good and difficult moments these last three years. Her encouragement and love was
immensely important and helped me finish this PhD. Next, I want to express my deepest
appreciation to my parents, who helped me through the tough times and gave me very
good advice. Apart from offering me work-related guidance, they were always available
when I called on them and were of great emotional comfort. I would also like to mention
my grandmother and my siblings, Paul, Hans and Naomi-Joan who also provided support,
each in their own way.

My friends and (some ex) housemates, Chris, Imogen and Sjoerd are also owed my
gratitude. The many discussions with Chris helped me solve various problems that I en-
countered throughout the PhD. The evening gaming sessions with Sjoerd (and sometimes
my brother Paul) helped me get rid of any frustrations that I might have built up during
the day.
# RESUMÉ EN FRANÇAIS

## Introduction

Distance Geometry (DG) est l'étude des concepts géométriques basés sur la distance. Un des premiers résultats historiques dans le contexte de DG est la formule de Heron, dont l'origine remonte au 1er siècle de notre ère. Cette formule peut être utilisée pour calculer l'aire d'un triangle uniquement en fonction de ses côtés (les distances entre ses sommets). De nombreux outils et techniques pratiques permettent de mesurer les distances entre objets, comme les capteurs des réseaux de capteurs [1] ou les expériences Nuclear Magnetic Resonance pour les distances inter-atomiques dans les protéines [2]. DG est la branche des mathématiques qui s'occupe de la façon dont nous pouvons obtenir plus d'informations sur l'ensemble des objets en utilisant uniquement des distances par paires. Dans cette thèse, nous étendons la discussion afin de considérer non seulement les distances euclidiennes, mais également d'autres concepts géométriques tels que les angles et les orientations locales ou différents types de distances. Le problème central dans le domaine de DG est le Distance Geometry Problem (DGP) [3], dans lequel nous sommes chargés de trouver les positions (relatives) des objets, de telle sorte que les distances mesurées entre eux soient satisfaites. Dans cette thèse, nous étudions différents aspects du problème de la géométrie des distances. Le DGP est un problème avec de nombreuses applications. Notre objectif principal réside dans l'application de la biologie structurale, où l'objectif est de déterminer la structure d'une protéine sur la base des distances interatomiques. Outre la détermination de la structure des protéines, nous discuterons également du problème du reciblage de mouvement, qui est lié à la dynamique DG [4, 5]. La dernière application que nous étudions est celle des cartes adaptatives. Notre travail pour cette application comprend le développement et la mise en œuvre d'un modèle linéaire qui montre des résultats prometteurs.

Il existe diverses méthodes et algorithmes qui peuvent trouver des solutions au DGP. Elles vont des méthodes basées sur la décomposition matricielle, l'optimisation globale, les heuristiques et les méta-heuristiques. Dans cette thèse, nous nous concentrons sur un algorithme qui exploite la discrétisabilité de certaines instances. Cette classe d'instance
est connue sous le nom de Discretizable Distance Geometry Problem (DDGP) [6]. Lorsque
nous discrétisons l'espace de recherche, nous permettons d'énumérer toutes les solutions
possibles au problème posé. L'algorithme capable d'effectuer cette tâche est branch-and-
prune (BP), qui fonctionne en phase de branchement et d'élagage. Notez que cet algo-
rithme ne peut énumérer complètement toutes les solutions possibles que lorsque nous
utilisons principalement des données de distance précises. Dans de nombreuses applica-
tions, nous sommes confrontés à un degré élevé de bruit, ce qui conduit à des intervalles de
distance. Il existe diverses implémentations de l'algorithme BP conçues pour gérer cette
incertitude. Dans cette thèse, nous nous concentrons sur un en particulier: MDJEEP [7].
Cette implémentation, initialement écrite en C, a été adaptée à Java pour cette thèse. De
plus, dans le cadre de cette thèse, cet algorithme a été étendu pour utiliser non seulement
les distances, mais également les angles de torsion. En particulier, elle a été étendue pour
exploiter le signe des angles dièdres dans la phase de branchement. Cela signifie que la
taille de l'espace de recherche est considérablement réduite, ce qui augmente les chances
de l'algorithme d'identifier des solutions de bonne qualité.

Biologie structurale

L'algorithme BP est particulièrement utile dans le contexte de la biologie structurale. Certaines instances qui surviennent dans le contexte d'expériences Nuclear Magnetic Resonance (NMR) peuvent être discrétisées [8]. De plus, dans ces cas, l'inclusion des informations sur l'angle de torsion devrait améliorer l'efficacité de la méthode BP. Pour étudier cela, nous avons présenté plusieurs expériences avec l'algorithme étendu BP, implémenté en Java. La première série d'expériences concerne des instances issues de données réelles NMR. Dans la littérature, il est courant d'étudier BP pour la biologie structurale en utilisant des données simulées au lieu de données réelles [7–14]. En effet, le bruit dans les données réelles est si élevé que les méthodes BP ont tendance à échouer. Les expériences présentées dans cette thèse montrent que, au moins pour les petites protéines, les structures générées à l'aide de la méthode BP étendue sont de bonne qualité. Cependant, à mesure que la taille des protéines augmente, l'incertitude capturée dans les instances augmente également. Cette incertitude a un impact important sur la qualité des résultats. Cela nous amène à la deuxième série d'expériences. Une source d'inexactitude dans les cas provient des variations dans la géométrie covalente des protéines. Nous présentons une analyse qui montre que les angles de liaison et l'angle de torsion ω peuvent présenter de
grands écarts d'un acide aminé à l'autre. Cependant, dans les cas que nous utilisons, ces variations ne sont pas capturées et ces angles (et les distances associées) sont souvent considérés comme constants. Dans les expériences que nous présentons, nous montrons l'effet de ces variations sur la méthode BP et les solutions générées. Dans les résultats, nous constatons une corrélation claire entre la taille des protéines et la qualité des solutions, montrant que pour des protéines plus grosses, l'algorithme produit des structures moins bonnes. Les travaux futurs sur ce sujet visent principalement à résoudre le problème de ces variations, ainsi qu'à améliorer la gestion de l'incertitude dans les instances.

## Géométrie de distance dynamique

Nous examinons ensuite la géométrie dynamique des distances en nous concentrant principalement sur les mouvements humains. Les mouvements humains sont des animations généralement enregistrées par capture de mouvement. Ces mouvements décrivent un personnage squelettique effectuant un mouvement à travers différentes images. La raison pour laquelle nous nous intéressons particulièrement aux mouvements humains est que, comme chez les humains, les protéines présentent également une structure squelettique. Par conséquent, avec cette étude, nous posons les bases de futurs projets qui pourraient se concentrer sur la modélisation de la dynamique des protéines. Dans cette section de la thèse, nous nous concentrons sur le problème du retargeting de mouvement [15–17]. Ici, l'enjeu est de copier le mouvement d'un acteur vers un deuxième personnage, de morphologie différente. Ce problème peut être décrit en utilisant Distance Geometry dynamique et des travaux antérieurs montrent que cette approche basée sur la distance peut être fructueuse [4, 5]. Nous étendons davantage l'approche basée sur la distance, en travaillant directement dans l'espace des angles d'Euler. De plus, nous proposons une nouvelle représentation des mouvements humains, inspirée de la biologie structurale, qui pourrait être particulièrement utile dans le contexte du reciblage des mouvements humains [18, 19]. En utilisant cette nouvelle représentation, nous présentons une analyse statistique, mettant en valeur la pertinence de la représentation. Les prochaines étapes sur ce sujet consistent notamment à améliorer encore la stratégie basée sur la distance pour le reciblage de mouvement. Deuxièmement, les travaux devraient être étendus à la modélisation de la dynamique des protéines.
# Cartes adaptatives

La dernière application sur laquelle nous nous concentrons est celle des cartes adaptatives [20, 21]. Ici, nous générons des cartes où les distances entre les Point of Interest (POIs) clés sur la carte ne sont pas euclidiennes, mais plutôt basées sur des mesures subjectives de proximité, qui dépendent de l'utilisateur. L'un des défis lors de la génération de cartes adaptatives est de déplacer les POI de manière à ce que ces distances subjectives soient satisfaites, tout en s'assurant que la nouvelle carte ressemble toujours à l'ancienne carte à usage général. Ce problème de relocalisation des POI peut être formulé comme un type de DGP, où nous incluons des contraintes basées sur les orientations locales entre les POI. Ces mêmes contraintes permettent de linéariser le problème. Pour cette thèse, un nouveau modèle linéaire a été proposé et implémenté. De plus, nous présentons diverses expériences avec le nouveau modèle, montrant des résultats prometteurs. Enfin, nous examinons une possibilité d'amélioration de la méthode, qui consiste à essayer différentes rotations des axes des orientations locales. Les travaux futurs se concentreront sur l'amélioration de cette méthode de rotation.

# Représentations binaires et DGP en dimension 1

Enfin, l'auteur a apporté des contributions mineures à deux autres sujets. Le premier travail se concentre sur les représentations binaires de problèmes combinatoires et sur la manière dont nous pouvons utiliser les méta-heuristiques pour trouver des solutions. Cela concerne DG car les solutions aux instances de la sous-classe DDGP peuvent être représentées à l'aide d'une chaîne binaire. Le deuxième travail considère le DGP en dimension 1. En particulier, nous nous sommes concentrés sur un type spécifique d'instances «paradoxales», qui semblent triviales mais sont en fait difficiles à résoudre par les méthodes branch-and-prune. Nous proposons des méthodes alternatives, basées sur une reformulation matrice par vecteur, qui permettent de résoudre plus efficacement ces instances paradoxales.
<figure><img src="image_3.png" /></figure>
# TABLE OF CONTENTS

<table><tr><td>List of Figures</td><td>18</td></tr><tr><td>List of Tables</td><td>19</td></tr><tr><td>List of Acronyms and Abbreviations</td><td>20</td></tr><tr><td>1 Introduction</td><td>23</td></tr><tr><td>1.1 Defining the Distance Geometry Problem</td><td>25</td></tr><tr><td>1.2 Different DGPs and their applications</td><td>28</td></tr><tr><td>1.3 Not only distances</td><td>31</td></tr><tr><td>1.4 Methods for solving the DGP</td><td>35</td></tr><tr><td>1.4.1 Operating on Euclidean Distance Matrices</td><td>35</td></tr><tr><td>1.4.2 Stochastic Proximity Embedding</td><td>38</td></tr><tr><td>1.4.3 Nonlinear Programming</td><td>40</td></tr><tr><td>1.4.4 Meta-heuristics</td><td>41</td></tr><tr><td>1.4.5 Spectral Projected Gradient</td><td>43</td></tr><tr><td>1.4.6 Distance norms and Linear Programming</td><td>46</td></tr><tr><td>1.4.7 Build-Up and Branch-and-prune</td><td>47</td></tr><tr><td>1.5 The Java Package</td><td>48</td></tr><tr><td>2 Discretizable Distance Geometry</td><td>52</td></tr><tr><td>2.1 Introduction</td><td>54</td></tr><tr><td>2.2 Branch-and-prune in dimension 3</td><td>57</td></tr><tr><td>2.3 Coordinate generation</td><td>60</td></tr><tr><td>2.4 Handling interval distances</td><td>63</td></tr><tr><td>2.4.1 Sampling the arcs</td><td>63</td></tr><tr><td>2.4.2 A coarse-grained representation and local optimization</td><td>64</td></tr><tr><td>2.5 Torsion angles and the search space</td><td>65</td></tr></table>
<table><tbody><tr><td><b>3</b></td><td><b>Molecular Structure Determination</b></td><td><b>68</b></td></tr><tr><td>3.1</td><td>The problem of Protein Structure Determination</td><td>70</td></tr><tr><td>3.2</td><td>Experimental techniques and machine learning</td><td>72</td></tr><tr><td>3.2.1</td><td>X-ray crystallography</td><td>72</td></tr><tr><td>3.2.2</td><td>Cryo-Electron Microscopy</td><td>73</td></tr><tr><td>3.2.3</td><td>Nuclear Magnetic Resonance Spectroscopy</td><td>74</td></tr><tr><td>3.2.4</td><td>Machine learning and ALPHAFOLD</td><td>76</td></tr><tr><td>3.3</td><td>Inter-atomic distances and torsion angles</td><td>77</td></tr><tr><td>3.4</td><td>Distance-based approaches</td><td>80</td></tr><tr><td>3.4.1</td><td>Simulated Annealing</td><td>81</td></tr><tr><td>3.4.2</td><td>General DGP methods</td><td>82</td></tr><tr><td>3.4.3</td><td>Exploiting the discretizability of the instances</td><td>84</td></tr><tr><td>3.5</td><td>Experiments and results</td><td>85</td></tr><tr><td>3.5.1</td><td>Branch-and-prune with real NMR instances</td><td>85</td></tr><tr><td>3.5.2</td><td>Covalent geometry and branch-and-prune</td><td>90</td></tr><tr><td>3.6</td><td>Discussion</td><td>97</td></tr><tr><td><b>4</b></td><td><b>Dynamical Distance Geometry & Motions</b></td><td><b>99</b></td></tr><tr><td>4.1</td><td>Introduction</td><td>101</td></tr><tr><td>4.2</td><td>Human motion retargeting</td><td>102</td></tr><tr><td>4.2.1</td><td>Representing human motions</td><td>103</td></tr><tr><td>4.2.2</td><td>Motion retargeting in the literature</td><td>105</td></tr><tr><td>4.3</td><td>Distance-based motion retargeting in Euler space</td><td>109</td></tr><tr><td>4.4</td><td>A representation for human motions inspired by structural biology</td><td>113</td></tr><tr><td>4.4.1</td><td>The vector-torsion representation</td><td>113</td></tr><tr><td>4.4.2</td><td>Motion analysis</td><td>117</td></tr><tr><td>4.5</td><td>Discussion & Future work</td><td>125</td></tr><tr><td><b>5</b></td><td><b>Adaptive Maps & Linear Programming</b></td><td><b>126</b></td></tr><tr><td>5.1</td><td>Introduction</td><td>128</td></tr><tr><td>5.2</td><td>A linear program for POI relocation</td><td>130</td></tr><tr><td>5.2.1</td><td>Implementing the model</td><td>131</td></tr><tr><td>5.3</td><td>Experiments: the city of Granada</td><td>132</td></tr><tr><td>5.4</td><td>Rotating the axes</td><td>135</td></tr><tr><td>5.4.1</td><td>Experiments</td><td>136</td></tr></tbody></table>
<table><tr><td>5.5</td><td>Discussion</td><td>139</td></tr><tr><td>6</td><td>Discussion, ongoing and future work</td><td>141</td></tr><tr><td>6.1</td><td>Discussion and future work</td><td>143</td></tr><tr><td>6.2</td><td>Ongoing work</td><td>145</td></tr><tr><td>6.2.1</td><td>Binary representations</td><td>146</td></tr><tr><td>6.2.2</td><td>The DGP in dimension 1</td><td>147</td></tr><tr><td colspan="2">Author's Contributions</td><td>149</td></tr><tr><td colspan="2">List of Publications</td><td>149</td></tr><tr><td colspan="2">Bibliography</td><td>151</td></tr></table>
2.4 In dashed lines, the circles obtained by intersecting the two spheres centered in $u_1$ and $u_2$. (a) The strategy of selecting a predefined number $B$ samples. In this example, $B = 3$, so six branches would be created. (b) The arc is smaller than the resolution parameter $\rho$. We fully cover the arc with a box and initialize the candidate position $v_+$ in the center of the arc. We create 2 branches. (c) The arc is larger than the resolution parameter $\rho$. We cover the arc with equally sized boxes and initialize the four candidates positions $v_+$ in the center of their sub-arcs. We create four branches. . . . . . . . . . . 65

3.1 Three different levels of protein organization. Left: primary structure (linear sequence of amino acids). Middle: secondary structure (structural patterns, alpha helix, and beta sheets). Right: tertiary structure (three-dimensional confirmation of proteins). Image modified of original provided by National Human Genome Research Institute . . . . . . . . . . . . . . . . . . . . . . . . . 70

3.2 An example of a two-dimensional NOESY spectrum for a small molecule, where the chemical shifts on the axes correspond to $^1$H protons (labeled in ppm). At the top of the spectrum, the different frequencies are labeled with the hydrogen that emits it. In the bottom right there is a two-dimensional diagram of the molecule that is being analysed. The red peaks in the diagonal are the self-peaks for each of the hydrogen atoms. The cross-peaks (blue) away from the diagonal are the ones we are interested in. For example, the indicated cross-peaks in the figure tell us that both $H_8$ atoms are in close proximity to the $H_4$ and $H_7$ atoms in the molecule. We see that there is a peak for $H_{8b}$ and the $H_3,H_9$ pair, but no peak between the pair and $H_{8b}$, possibly indicating that $H_{8b}$ is not close to $H_3$ and $H_9$. Figure taken from [168] and modified. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

3.3 The structures computed by Java MDjeep (blue) with lowest RMSD superimposed with the structures from the PDB. The structures were not subject to energy minimization or any other post-processing step. . . . . . . . . . . . 89
3.4 The average values of the bond angles between the heavy backbone atoms (N, CA, C, O) as well as on the torsion angle $\omega$ of the peptide plane. Np denotes the nitrogen atom N of the next residue, Cm denotes the carboxyl carbon C of the previous residue, and CA denotes the CA atom. The distributions of these angle values are uni-modal, which means that we can use their averages as a suitable descriptor. On the left, we see the averages of the angles according to the amino-acid types, while on the right we see them according to the Ramachandran region of the residue. The horizontal dashed lines correspond to the values of the angles from the force field parameters [180]. 92

3.5 The Ramachandran plot illustrating the 12 different regions used in the data analysis. Note that this figure is only an approximation of the regions based on a 5-degree grid. 93

3.6 A histogram with the resulting RMSD scores for the two instance classes, with bins of size 2Å. 96

3.7 The two structures computed by Java MDJEEP for instances A (pink) and B (blue) overlapped with the structure measured with X-ray crystallography (gold). Examples for four proteins: 1EZG, 2ZW2, 1MK0 and 5YIU. Underneath, we see the PDB code of the protein as well as the RMSD of the two produced structures. 97

4.1 An example of a skeletal structure ($H, \chi$). The initial T-pose along with the labels are commonly used in BVH files [199]. Joints with $|\chi| = 0$ are shown next to their parent joint, marked in italics. For the arm, only the left-side joints are labeled, for the legs only the right-side labels are shown. The root of the tree $H$ is the “Hips” joint. Note that this is the original root of $H$, and not the fictive root node $v_0$, which is instead placed at the origin of the coordinate system. This image and other human motion images that we will present in this thesis are generated using the free software Blender. 104
4.7 An example of $\zeta_v$ and $\tau_v$. Visible are the legs of a human skeleton, taken from a BVH file describing a running motion. The left side illustrates the vector angle $\zeta_v$, which is the angle between the RightFoot, the RightLeg, and the RightUpLeg joints. The center shows how we arrive at the torsion angle $\tau_v$ using two planes defined by RightFoot and its three ancestors RightLeg, RightUpLeg, and Hips (or RightHips). On the right, we see what happens when $\tau_v = 270^\circ$, while $\zeta_v$ remains $90^\circ$. 116

4.8 The RightLeg and the LeftUpLeg are two examples of joints that do not have enough ancestors to be represented by the vector-torsion angle pair. We can still use other consecutive reference vertices in the anatomy *H* to compute these angles and perform the same analysis for these joints. In the figure, we show the two sets of joints used to define the angles for the RightLeg and the LeftUpLeg. For the former, we used the joints RightUpLeg, Hips, and Spine (the dark blue path). For the latter, we used the joints Hips, Spine, and Spine1 (the orange path). 118

4.9 The plot for the Head. 119

4.10 The plot for the Neck. The encircled points correspond to the unnatural position seen on the right in Figure 4.27. 119

4.11 The plot for the Spine1. 119

4.12 The plot for the Spine. 119

4.13 The plot for the LeftArm. 120

4.14 The plot for the RightArm. 120

4.15 The plot for the LeftForeArm. 120

4.16 The plot for the RightForeArm. The encircled points correspond to the unnatural position seen on the right in Figure 4.27. 120

4.17 The plot for the LeftHand. 121

4.18 The plot for the RightHand. 121

4.19 The plot for the LeftUpLeg. 122

4.20 The plot for the RightUpLeg. 122

4.21 The plot for the LeftLeg joint. The points in the green circle correspond to the uncommon poses in Figure 4.28. 122

4.22 The plot for the RightLeg joint. The points in the green circle correspond to the uncommon poses in Figure 4.28. 122

4.23 The plot for the LeftFoot. 122
<table>
  <tbody>
    <tr>
      <td>4.24</td>
      <td>The plot for the RightFoot.</td>
      <td>122</td>
    </tr>
    <tr>
      <td>4.25</td>
      <td>The plot for the LeftToeBase.</td>
      <td>123</td>
    </tr>
    <tr>
      <td>4.26</td>
      <td>The plot for the RightToeBase.</td>
      <td>123</td>
    </tr>
    <tr>
      <td>4.27</td>
      <td>The two motion capture errors. On the left, we see a bug for the Neck (motion: #138_31, frames: 70-688). The corresponding points in the scatter plot are encircled in red in Figure 4.10. On the right, we see an error for the RightForeArm (motion: #143_38, frames: 625-660). The corresponding points in the scatter plot are encircled in red in Figure 4.16.</td>
      <td>124</td>
    </tr>
    <tr>
      <td>4.28</td>
      <td>Examples of outlying data points that correspond to valid human postures. Left: a frame of the character landing after a wide-legged jump (motion: #121_19, frame: 847). Right: a frame of the character in a sitting position (motion: #111_09, frame: 357).</td>
      <td>125</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Information related to the 10-point instance in [21]. Top left: a visual representation of the local orientations of the POIs in the original realization, which our linear model is supposed to preserve in the found solution. Bottom right: a heatmap, where the colors encode the difference in values between the Euclidean distances (in the original realization x<sub>0</sub>) and the given proximity distances (walking distances, see Table 5.1). The lower the difference, the smaller the value attributed to the corresponding square in the heatmap; all difference values are normalized between 0 and 1.</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>The solution found by our linear model to the POI relocation problem related to data in Table 5.1. The points P<sub>*</sub> indicate the original locations of the POI in the initial map. In black, we see the new positions for the POIs.</td>
      <td>134</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>An analysis of the data set with 10 POIs (see Table 5.1). The x-axis: the virtual axes rotation in degrees. On the y-axis: the objective value of the corresponding LP.</td>
      <td>136</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Overlaps the original points (labeled) with the two sets of resulting points. In black, we see the solution with an objective value of 10273, and in gray we see the solution with an objective value of 7896. Gridlines are added for the points in black, showing the relative constraints relating to the pair (P4,P8) and (P1,P10) are swapped.</td>
      <td>137</td>
    </tr>
  </tbody>
</table>
LIST OF TABLES

<table>
  <tbody>
    <tr>
      <td>2.1</td>
      <td>An overview of the similarities and differences between the three implementations.</td>
      <td>59</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Different properties of the proteins and peptides considered for these experiments. Includes the number of distances, per type, that can form our DDGP instances.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>The resulting RMSD scores for the proteins (in Å).</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Details for the two instance classes that were experimented with. Instances A are completely artificial while instance B mixes artificial with information from force-field parameters.</td>
      <td>94</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>The resulting RMSD statistics for the 435 proteins for instances A and B.</td>
      <td>95</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>The three reference joints for each joint that was involved in the analysis. The reference joint with the smallest numerical label is the closest; the one with the largest numerical label is instead the farthest. In some cases, the reference joints are not the joint ancestors given by the graph structure.</td>
      <td>118</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>The walking distances between 10 pairs of POIs, that we use in our computational experiment. Reproduced from [21] (Table A1).</td>
      <td>133</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Results of using the exhaustive rotation method on various point sets.</td>
      <td>138</td>
    </tr>
  </tbody>
</table>
LIST OF ACRONYMS AND
ABBREVIATIONS

<table><tr><td>APA</td><td>Alternating Projections Algorithm.</td></tr><tr><td>BH</td><td>Basin Hopping.</td></tr><tr><td>BP</td><td>branch-and-prune.</td></tr><tr><td>CNNs</td><td>convolutional neural networks.</td></tr><tr><td>DDGP</td><td>Discretizable Distance Geometry Problem.</td></tr><tr><td>DG</td><td>Distance Geometry.</td></tr><tr><td>DGP</td><td>Distance Geometry Problem.</td></tr><tr><td>DVOP</td><td>Discretization Vertex Order Problem.</td></tr><tr><td>dynDGP</td><td>Dynamical Distance Geometry Problem.</td></tr><tr><td>EDM</td><td>Euclidean Distance Matrix.</td></tr><tr><td>GO</td><td>Global Optimization.</td></tr><tr><td>LP</td><td>Linear Programming.<br>linear program.</td></tr><tr><td>MBH</td><td>Monotonic Basin Hopping.</td></tr><tr><td>MDGP</td><td>Molecular Distance Geometry Problem.</td></tr><tr><td>MDS</td><td>Multi-dimensional scaling.</td></tr><tr><td>MIP</td><td>Mixed-integer programming.</td></tr><tr><td>MSAs</td><td>Multiple Sequence Alignments.</td></tr></table>
<table><tr><td>NeRF</td><td>Natural Extension Reference Frame.</td></tr><tr><td>NLP</td><td>Nonlinear Programming.</td></tr><tr><td>NMR</td><td>Nuclear Magnetic Resonance.</td></tr><tr><td>NOE</td><td>Nuclear Overhauser effect.</td></tr></table>

<table><tr><td>PBH</td><td>Population Basin Hopping.</td></tr><tr><td>PCA</td><td>Principal Component Analysis.</td></tr><tr><td>PDB</td><td>Protein Data Bank.</td></tr><tr><td>POI</td><td>Point of Interest.</td></tr><tr><td>POIs</td><td>Points of Interest.</td></tr></table>

<table><tr><td>Revorder</td><td>Referenced Vertex Ordering.</td></tr><tr><td>RMSD</td><td>Root mean square deviation.</td></tr></table>

<table><tr><td>SA</td><td>Simulated Annealing.</td></tr><tr><td>SBB</td><td>Spatial Branch-and-Bound.</td></tr><tr><td>SDP</td><td>Semidefinite Programming.</td></tr><tr><td>SoME</td><td>Summer of Math Exposition.</td></tr><tr><td>SPE</td><td>Stochastic Proximity Embedding.</td></tr><tr><td>SPG</td><td>Spectral Project Gradient.</td></tr><tr><td>SPM</td><td>Successive Projection Methodology.</td></tr><tr><td>VNS</td><td>Variable Neighbourhood Search.</td></tr></table>
<figure><img src="image_4.png" /></figure>
INTRODUCTION

Distance Geometry (DG) is the study of geometric concepts based on distance. An early historical result in the context of DG is Heron’s formula, which was proven in the 1st century CE. This formula can be used to compute the area of a triangle solely based on its sides (the distances between its vertices). Many practical tools and techniques allow for the measurement of distances between objects, such as the sensors in sensor networks [1] or Nuclear Magnetic Resonance experiments for inter-atomic distances in proteins [2]. DG is the branch of mathematics occupied with how we can obtain more information about the set of objects by using only pairwise distances. In this thesis, we extend the discussion so that we do not only consider Euclidean distances, but also other geometric concepts such as angles and local orientations or different distance types. The core problem within the field of DG is the Distance Geometry Problem (DGP), in which we are tasked with finding the (relative) positions of objects, such that the measured distances between them are satisfied. This thesis is focused on different variants of the DGP, its discretization, and some of its major practical applications, with an emphasis on structural biology.

In this introductory chapter, we start by giving formal definitions for the DGP (Sec-
tion 1.1). Next, we will touch on the applications of the DGP (Section 1.2). Afterwards,
in Section 1.3 we show that sometimes it is desirable to include other geometric informa-
tion in the instances besides distances. We will discuss different methods of solving the
general DGP found in the literature in Section 1.4. Finally, in Section 1.5, we will describe
the Java package that was created and used for the various experiments conducted as part
of this thesis. The second chapter of the thesis will focus on one of the most interesting vari-
ants of the DGP: the Discretizable Distance Geometry Problem (DDGP), whose instances
arise in the context of molecular biology. The rest of the thesis is structured according to
three applications of the DGP, to each of which a chapter is dedicated (Chapters 3 to 5).
Chapter 3 focuses on structural biology, Chapter 4 discusses motion retargeting and hu-
man motions and Chapter 5 explores the application of adaptive maps. Finally, Chapter 6
closes off with a short summary and a brief discussion of some ongoing work to which the
author has made minor contributions as part of the thesis.
<table><tr><td>1.1</td><td>Defining the Distance Geometry Problem</td><td>25</td></tr><tr><td>1.2</td><td>Different DGPs and their applications</td><td>28</td></tr><tr><td>1.3</td><td>Not only distances</td><td>31</td></tr><tr><td>1.4</td><td>Methods for solving the DGP</td><td>35</td></tr><tr><td>1.4.1</td><td>Operating on Euclidean Distance Matrices</td><td>35</td></tr><tr><td>1.4.2</td><td>Stochastic Proximity Embedding</td><td>38</td></tr><tr><td>1.4.3</td><td>Nonlinear Programming</td><td>40</td></tr><tr><td>1.4.4</td><td>Meta-heuristics</td><td>41</td></tr><tr><td>1.4.5</td><td>Spectral Projected Gradient</td><td>43</td></tr><tr><td>1.4.6</td><td>Distance norms and Linear Programming</td><td>46</td></tr><tr><td>1.4.7</td><td>Build-Up and Branch-and-prune</td><td>47</td></tr><tr><td>1.5</td><td>The Java Package</td><td>48</td></tr></table>
## 1.1 Defining the Distance Geometry Problem

The Distance Geometry Problem (DGP) [3] is a problem that involves determining the positions of a set $V$ of $n$ points in a $K$ dimensional space, based on a (possibly incomplete) set $\mathcal{D}$ of pairwise proximity measures between these points. In general, $\mathcal{D}$ is a set of distances $\delta(u, v)$ for different pairs of points $(u, v) \in V$. Historically, this problem has had different names [22], such as the *position-location problem* [23] or the *K-embeddability problem* [24]. In the context of sensor network localization, the problem is often referred to as the *Graph Realization problem* [1]. The DGP was proven to be NP-complete for $K = 1$, while for $K > 1$, it was shown to be NP-hard [24]. The most studied variant of the DGP is the *assigned version*, where the assignment of every distance $\delta \in \mathcal{D}$ to a corresponding pair of points $\{u, v\} \in V \times V$ is known. There also exists the less researched *unassigned DGP* [25]. In this thesis, we will focus on the *assigned variant*.

An instance of the assigned DGP may be represented in different ways. One of the most widely used representations in the literature [3, 18, 19, 22, 26–30], and the representation that we will use throughout this thesis, utilizes a graph structure to describe the problem. For this representation, we use a simple, undirected graph $G = (V, E, d)$. The vertices of $G$ correspond to the set of points that we want to find a realization for. The existence of an edge $\{u, v\} \in E$ in $G$ indicates that we have access to the distance from $u$ to $v$. Finally, $d$ is the weight function in $G$ which pairs every edge $\{u, v\} \in E$ to its corresponding distance $\delta \in \mathcal{D}$:

$$d: \{u, v\} \in E \rightarrow \delta(u, v) \in \mathcal{D}. \qquad (1.1)$$

Using the graph representation, the DGP can be formalized as follows [3, 22]:

**Definition 1** Given a simple weighted undirected graph $G = (V, E, d)$ and an integer $K \in \mathbb{Z}_+$ the (assigned) DGP asks whether a realization

$$x: v \in V \rightarrow x_v \in \mathbb{R}^K \qquad (1.2)$$

exists, such that

$$\forall \{u, v\} \in E, \quad ||x_u - x_v|| \in \delta(u, v) \qquad (1.3)$$

There are two important details to note here. Firstly, it is natural to think of the distances $\delta$ and the norm $||\cdot||$ in the constraints (1.3) to be Euclidean. However, there are
different methods in which the proximity between objects is estimated via non-Euclidean measures. Despite this, we will still attempt to find a realization in the Euclidean space. For example, in Section 1.4.6 and Chapter 5, we will see that in order to linearize the DGP, we need to use a linear distance measure, such as the $L_1$ and $L_\infty$ norms [31, 32]. Secondly, various applications will have distances that may include noise (for example, distances obtained through physical experiments). For this reason, the constraints (1.3) leave space for the distances $\delta(u, v)$ to be non-exact (an interval) $[\underline{\delta}(u, v), \bar{\delta}(u, v)]$:

$$ \underline{\delta}(u, v) \leq \|x_u - x_v\| \leq \bar{\delta}(u, v), $$

as well as the exact case (degenerate interval):

$$ \|x_u - x_v\| = \delta(u, v). $$

In some cases, where we have a lot of noise on the input distances, it may not be possible to find a realization such that all distances are satisfied. In these cases, it can instead be useful to consider the DGP as a Global Optimization (GO) problem, where we aim to minimize the error on the distances [26]:

**Definition 2** Given a simple weighted undirected graph $G = (V, E, d)$ and an integer $K \in \mathbb{Z}_+$ the unconstrained optimization-based DGP asks to determine a realization

$$ x: v \in V \rightarrow x_v \in \mathbb{R}^K $$

such that a penalty function $\sigma$ is minimized.

The penalty function $\sigma$ should measure the violation of the distance constraints. For example, we may use [4, 33]:

$$ \sigma(x) = \frac{1}{2} \sum_{\{u,v\} \in E} (\|x_u - x_v\| - \delta(u, v))^2 \quad (1.4) $$

Note that when all distances are exact, the problems given in Definitions 1 and 2 are equivalent. An important property of $\sigma(x)$ is that it is differentiable at $x$ when no two points are in the same position. Thus we can compute the gradient $\nabla\sigma(x)$ and use methods such as gradient descent to find local minima. In Section 1.4, we will see several optimization techniques that find solutions for the DGP by minimizing such penalty functions.
When the distance graph $G$ is complete, or at least dense, it may make sense to represent the DGP using an Euclidean Distance Matrix (EDM) [34–36], which is a hollow square matrix containing the pairwise squared distances between the points in $V$. Consider some order on the points in $V = \{v_1, \dots, v_n\}$, where $n = |V|$. The entry of the EDM at indices $(i, j)$ contains the squared distance $\delta(v_i, v_j)^2$. This is a matrix of real values, which means an EDM cannot be used to represent interval distances, but only exact ones. This may be too limiting for some of the applications of the DGP. Figure 1.1 showcases the graph representation of the DGP and compares it to the EDM representation.

<figure><img src="image_5.png" /><figcaption>Figure 1.1 – On the left we see the graph representation of the DGP, where the presence of an edge represents the fact that we know a distance. The colours of the edges represent the value of the distance. On the right, we see the EDM representation of the same two instances. On the top, we have a fully complete instance where we have access to all pairwise distances. On the bottom, we see a cycle graph, where we miss half of the distances, leading to a sparser matrix. Inspiration was taken from Fig. 1 in [36].</figcaption></figure>

Note that in the graph representation of the DGP, we do not have a specific vertex order on $V$, while the EDM does include an implicit vertex ordering. The advantage of the EDM representation is that we can use several techniques to find an embedding for $V$ that
operates directly on the EDM (see Section 1.4.1). However, there are several complications with these techniques. In many real-world applications, we will have missing distances, and thus the EDM will be incomplete. One can only use the matrix-based techniques to find a realization of the distances when the EDM is complete. This leads to a problem known as *EDM completion*, where efforts are made to complete the missing distance information in the matrix, for example through the use of Semidefinite Programming. More information on these techniques is given in Section 1.4.1.

## 1.2 Different DGPs and their applications

Distance Geometry Problems can come in many forms. Instances of these different types of problems can be found in different practical applications. Common themes are the notions of discretization, graph rigidity, non-Euclidean distance norms, and the fact that in practice we often have a large degree of noise on the input distances. As part of this thesis, practical work was done on four of these applications (molecular structure determination, motion retargeting, adaptive maps, and clock synchronization).

A very important application, and the main focus of this thesis, is *molecular structure determination* [27, 33, 37–45]. In this application, the vertices *V* of the graph *G* correspond to atoms of a molecule, and the distance information *d* associated with the edges *E* of the graph are the inter-atomic distances. The goal is to determine the three-dimensional structure of the molecule such that the inter-atomic distances are satisfied. The instances of the DGP in the context of molecular structures often belong to a special subclass of the problem. The distance graphs that we obtain satisfy special conditions that let us discretize the search space [6, 8, 26]. This means that we may be able to enumerate all three-dimensional structures that a protein can fold into. These types of problems are referred to as instances of the Discretizable Distance Geometry Problem (DDGP). More detail on this instance type can be found in Chapter 2. The problem of molecular structure determination is further examined, where several experiments will be presented. In the literature, most experiments with DDGP methods for protein structure determination are done with artificially generated instances from known structures. In Chapter 3, experiments will be presented with instances created from real protein data. An example of some of the results of these experiments can be seen in Figure 1.2, which shows the structures of two proteins with complex secondary structures, that were computed by exploiting several available distances between the atoms of the molecules.
<figure><img src="image_6.png" /><figcaption>Figure 1.2 – Examples of conformations of two proteins (PDB codes from left to right: 1YD0A, 4WZXA). These structures were computed using the inter-atomic distances and other geometric information, using software created as part of this thesis. The image was generated using molecule viewing software Chimera [46].</figcaption></figure>

It is known that, because of the inter-atomic forces at play, the atoms of a protein move over time, forcing the protein to take on different conformations. Instances of the DGP where we have such a dynamic aspect, where the points that we want to embed may move over time, are referred to as instances of the Dynamical Distance Geometry Problem (dynDGP). An example of an application of this dynamical variant is *motion retargeting*, where we want to copy the animation of one character to another while avoiding unwanted contact between the joints of the character. To do this, we may exploit the distances between the joints in the original character, and attempt to map these to the second animation. To solve such a dynDGP, essentially we have to solve a series of DGP instances. In the case of motion retargeting, we are solving one DGP instance per frame. Aside from motion retargeting, the dynDGP has other applications such as simulating pedestrians in a crowd [47, 48] or aircrafts in the air [49]. More information about the dynDGP will be given in Chapter 4, in which several experimental results will be presented on the problem of motion retargeting.

The concept of discretizability, which was mentioned in relation to the application of molecular structure determination, is closely linked to the notion of *graph rigidity* [23, 50–54]. In fact, if a graph is rigid, we know the number of possible incongruent realizations of the graph is finite [22, 55]. A graph is said to be rigid if its edges will not *bend* or
*flex* under an applied outside force. This leads us to a historical application of Distance Geometry (DG): statics [56, 57], which is the study of forces that act on non-moving bodies in physical systems that are in static equilibrium. More details of how this application relates to DG can be found in [22].

Another interesting application of the DGP is *wireless sensor network localization* [1, 58–63], which is a problem where rigid graph theory is very relevant [64]. Here, we have a set of wireless sensors in a network. To properly route communication signals through the network, the router must determine the positions of the sensors. Every sensor, given a certain radius, is aware of its neighbours. This awareness is expressed in the form of a distance, which means that we have an instance of the DGP. In our graph representation for this problem, the vertices $V$ are the sensors and edges $E$ reflect the distances between them. These sensors are prone to making errors, so similarly to when we deal with inter-atomic distances for molecular structure determination, the measured distances between the sensors may contain a lot of uncertainty. Rigid graphs have the property of sustaining various kinds of deformations due to translation, rotation, and reflection. Hence, it is fruitful to use the concepts of rigid graphs, to find the locations of the sensors, even when we have noise on the distances between them.

A classical application of the Distance Geometry Problem is *dimensionality reduction* [65–70]. Dimensionality reduction is the transformation of data points in a high-dimensional space to a lower-dimensional space. In the context of Distance Geometry, the main idea is to start by computing the distances between the data points in the input. This is an example of a case in which we use distance norms other than the Euclidean one. This is because distances computed at higher dimensions are likely to be non-realizable in lower dimensional spaces. Therefore, other distance norms may be used, such as the $L_1$ norm (see Section 1.4.6 for a definition of this norm). Once we have the distances between the data points, we can solve an instance of the DGP with a different target dimension $K$, lower than the original dimension of the data points. The DG-based techniques for solving dimensionality reduction problems are known as Multi-dimensional scaling (MDS) [34, 36, 65], and work directly on EDMs (see Section 1.4.1).

Another application of Distance Geometry in which non-Euclidean distance norms are important is the more recently studied problem of *Adaptive maps*. More often than not, geographical maps are meant for general use. Because of this, the distances between Points of Interest (POIs) on modern maps are based on the Euclidean distance between them. Adaptive maps on the other hand are maps that are modified based on characteristics
specific to the end-user, such as the slope of the route or the accessibility of the user. To create such an adaptive map, the original geographical map needs to be modified to match these new user-specific distances. This corresponds to solving a DGP in dimension 2. To make sure the adaptive map still bears a resemblance to the original, some extra linear constraints are added to the problem. This makes Linear Programming one of the best tools for finding solutions to the problem. More information on Linear Programming as a tool to solve DGPs can be found in Section 1.4.6 while more details on adaptive maps and several related computational experiments can be found in Chapter 4.

In dimension 1, the $L_1$, $L_2$, and $L_\infty$ norms are the same. In this dimension, we have the *clock synchronization* problem, where we are given a set of clocks, a subset of offset measurements between pairs of clocks as well as some offsets to a central clock, for which we know the time. The goal is to compute the precise time on each of the clocks [61, 71, 72]. This is the main practical application for the DGP in dimension $K = 1$. The clocks represent the vertices $V$ of our graph $G$ and the edges $E$ and their weight $d$ corresponds to the temporal offsets between the clocks. Section 6.2.2 focuses on this problem and the DGP in dimension 1. Specifically, we focus on one special difficult subclass of the DGP in dimension 1, in which our distance graph $G$ is a cycle graph.

One of the reasons that the DGP is an interesting and multifaceted problem is that applications seemingly appear everywhere. For example, a niche application, closely related to sensor networks, is that of controlling and locating autonomous underwater vehicles [73, 74]. An application of the DGP can even be found in music [75], where Distance Geometry is used to provide a quantitative analysis of rhythm and music. Aside from these, DG is also relevant for visualization [76, 77] and graph coloring [78].

## 1.3 Not only distances

In some of the applications mentioned in the previous section, we may have more geometric data available that is not captured by the input distances. Therefore, it is desirable to include this extra information in the instances so that we may exploit it.

### Vector and torsion angles

In the context of molecular structure determination, we have access to two types of angles: *vector* and *torsion* angles. Vector angles $\zeta(v_1, v_2, v_3) \in [0, \pi/2]$ are induced by three vertices $v_1, v_2$ and $v_3$, where we know $\delta(v_1, v_2)$ and $\delta(v_2, v_3)$. These vector angles can
be easily converted to a third distance $\delta(v_1, v_3)$ by using the cosine rule, without losing
any information. In the case of proteins, these angles come in the form of bond angles,
which are angles created by every pair of atomic bonds that include a common atom.

In contrast to these types of angles, torsion angles cannot be so easily translated to
distance information. Given four vertices $v_1, v_2, v_3$ and $v_4$. The torsion angle $\tau \in [-\pi, \pi]$ is the dihedral angle created by the two planes induced by $v_1, v_2, v_3$ and $v_2, v_3, v_4$. From this point onward, we will confuse the two terms dihedral and torsion angle.

Examples of one torsion angle and two vector angles induced by four vertices are shown
in Figure 1.3. In the case we know the three distances $\delta(v_1, v_2)$, $\delta(v_2, v_3)$ and $\delta(v_3, v_4)$, as
well as the vector angles $\zeta(v_1, v_2, v_3)$ and $\zeta(v_2, v_3, v_4)$ we may convert the torsion angle to
a distance $\delta(v_1, v_4)$.

One way to do this is by computing the coordinates of the points $v_1$ and $v_4$ with the
origin set at the midpoint of the vector $\overrightarrow{v_2v_3}$. We imagine the particular case where the
torsion angle $\tau = 0$. Then, from this specific case, we generalize for all other cases by
rotating the vector $\overrightarrow{v_1v_2}$ around the $x$ axis, with an angle of $\tau/2$. We rotate $\overrightarrow{v_3v_4}$ in the
opposite direction, with angle $-\tau/2$. This gives us the below formulas for $v1$ and $v4$:

$$v_1 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \tau/2 & -\sin \tau/2 \\ 0 & \sin \tau/2 & \cos \tau/2 \end{pmatrix} \begin{pmatrix} -\delta(v_1, v_2) \cos \zeta_1 \\ +\delta(v_1, v_2) \sin \zeta_1 \\ 0 \end{pmatrix} - \begin{pmatrix} \delta(v_2, v_3)/2 \\ 0 \\ 0 \end{pmatrix},$$

$$v_4 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \tau/2 & \sin \tau/2 \\ 0 & -\sin \tau/2 & \cos \tau/2 \end{pmatrix} \begin{pmatrix} -\delta(v_3, v_4) \cos \zeta_2 \\ +\delta(v_3, v_4) \sin \zeta_2 \\ 0 \end{pmatrix} + \begin{pmatrix} \delta(v_2, v_3)/2 \\ 0 \\ 0 \end{pmatrix}.$$

All that remains is to compute the distance from $v_1$ to $v_4$ using the Euclidean distance
formula.

However, converting such a torsion angle to a distance is not a lossless process. In
fact, $+\tau$ and $-\tau$ lead to the same distance $\delta(v_1, v_4)$. It is therefore desirable to include
the torsion angles directly in the instances and exploit them in other ways, such that the
information of the sign is not lost.
<figure><img src="image_7.png" /><figcaption>Figure 1.3 – An illustration of two vector angles $\zeta_1 = \zeta(v_1, v_2, v_3)$ and $\zeta_2 = \zeta(v_2, v_3, v_4)$ as well as the dihedral angle $\tau$</figcaption></figure>

To achieve this, we represent the problem using a directed weighted hypergraph $G(V, E, d)$.
In such a graph, an edge (called hyperedge) can join any number of vertices. Such a hyperedge $e \in E$ consists of a pair of subsets of $V$. In our hypergraph we will only have two types of edges, those that relate to distances and those that relate to torsion angles. We will refer to the edges that correspond to distances as $E_\delta \subset E$. An edge that relates to the distance between the vertices $u$ and $v$ is simply represented using a pair of singletons: $(\{u\}, \{v\}) \in E_\delta$. For every hyperedge $(\{u\}, \{v\}) \in E_\delta$, we should have the corresponding hyperedge $(\{v\}, \{u\}) \in E_\delta$. The edges that relate to torsion angles concern quadruplets of vertices which together (in a specific order) form a torsion angle $\tau$. These edges $E_\tau \subset E$ will consist of a pair of a triplet of vertices with a singleton. For example, if a hyperedge relates to a torsion angle formed by the quadruplet $v_1, v_2, v_3$ and $v_4$, we will have an edge $(\{v_1, v_2, v_3\}, \{v_4\}) \in E_\tau$. Note that these (directed) edges imply a specific vertex order. Vertex orders and their importance are discussed in detail in Chapter 2. In order to complete our description, the hypergraph is paired with the weighting function $d$. The function $d$ maps each hyperedge $(\{u\}, \{v\}) \in E_\delta$ to an interval distance $\delta(u, v)$, where $\delta(u, v) = \delta(v, u)$. Every edge $(\{v_1, v_2, v_3\}, \{v_4\}) \in E_\tau$ is mapped to an interval $[\underline{\tau}(v_1, v_2, v_3, v_4), \bar{\tau}(v_1, v_2, v_3, v_4)]$ described by $\tau(v_1, v_2, v_3, v_4)$. Note that just like before, the intervals $\delta(u, v)$ and $\tau(v_1, v_2, v_3, v_4)$ may be degenerate, i.e., they may correspond to an exact value.
This hypergraph lets us define a new variant of the DGP, in which torsion angles given in the input should also be observed:

**Definition 3** Given a weighted directed hypergraph graph $G = (V, E, d)$, where $d$ maps edges $E$ to distances $\delta$ and torsion angles $\tau$, the (torsion angle enriched) DGP asks whether a realization

$$x: v \in V \rightarrow x_v \in \mathbb{R}^3 \qquad (1.5)$$

exists, such that

$$\forall (\{u\}, \{v\}) \in E_{\delta}, \quad \|x_u - x_v\| \in \delta(u, v)$$

and

$$\forall (\{v_1, v_2, v_3\}, \{v_4\}) \in E_{\tau}, \quad \text{dihed}(x_{v_1}, x_{v_2}, x_{v_3}, x_{v_4}) \in \tau(v_1, v_2, v_3, v_4)$$

where dihed($x_{v_1}, x_{v_2}, x_{v_3}, x_{v_4}$) corresponds to the dihedral angle that is formed by the positions of $v_1, v_2, v_3$ and $v_4$.

In Section 2.5, we will see that including this extra information into our instances will let us decrease the size of the search space in the context of the branch-and-prune (BP) method. Furthermore, experiments in Section 3.5.2 will illustrate the usefulness of torsion angles in the context of molecular structure determination. Molecular structure determination is not the only application for which we are interested in these angles. In Chapter 4, we will see that they are also of interest for motion retargeting, where they allow us to define a new representation for human motions.

## Local relative orientations

Instances of the DGP relating to adaptive maps showcase another interesting case where we have additional geometric information. Because we have access to the original geometric map, we have access to the relative orientations between each pair of points that we want to embed. Specifically, we have access to the relative orientations for every pair of vertices $(u, v) \in E$. In the case of adaptive maps where the dimension $K = 2$, these orientations are included in the instances in terms of linear constraints by using four quadrants (NW, NE, SW, and SE). The edge set $E$ is partitioned into four subsets that correspond to these quadrants. For example, when we have $(u, v) \in \text{NW}$ we know that the positions for $v$ should be to the *North-West* of the positions for $u$. Including these constraints let us define a new variant of the DGP, where the local orientations
are imposed on the desired realization $x$. This variant of the problem can be solved by constrained optimization methods. The linear nature of these constraints makes Linear Programming an important contender for solving these instances. More information can be found in Section 1.4.6 and Chapter 5.

# 1.4 Methods for solving the DGP

Many algorithms have been proposed for solving the Distance Geometry Problem. There are methods using matrix operations directly on the EDM. We also have the option to use fast-running non-linear optimization methods that work mostly well on small instances. For applications where it makes sense to linearize the problem, we can opt for Linear Programming to find good-quality solutions effectively. Finally, if our instances meet specific conditions we can use geometric build-up methods that have the capability of providing us with all possible realizations to the problems at hand. In this section, we will describe the most important methods and frameworks. We also touch on the strengths and weaknesses of each of these methods.

## 1.4.1 Operating on Euclidean Distance Matrices

We are given a collection of $n$ points in dimension $K$, described by a coordinate matrix $X$ with $K$ rows and $n$ columns, where every column represents the set of coordinates of one of the points. We can now exploit the Euclidean distance formula to devise a function that computes the Euclidean Distance Matrix given such a (centered) coordinate matrix $X$ [36]:

$$ \mathrm{edm}(X) = \mathbf{1}\mathrm{diag}(X^{\top}X)^{\top} - 2X^{\top}X + \mathrm{diag}(X^{\top}X)\mathbf{1}^{\top} \quad (1.6) $$

where $\mathbf{1}$ is a column vector of ones and $\mathrm{diag}(A)$ is a column vector with the diagonal entries of the matrix $A$. Solving the DGP, given an EDM, corresponds to inversing the equation shown in (1.6), as long as all the distances are exact and we are considering the Euclidean norm. This is precisely the approach that Multi-dimensional scaling (MDS) [34, 36, 65] methods take to find solutions for the DGP. Recall from Section 1.1 that when we have instances with interval distances, we cannot use an EDM to represent the problem. Therefore we cannot use classical MDS techniques for these instances. However, as we will see towards the end of this section, EDM completion methods can be used to deal with
instances with interval distances. Classical MDS was implemented as part of the Java package for this thesis and included as a benchmark for other methods (see Section 1.5).

Note that in (1.6), we see that the EDM can actually be defined as a function of the Gram matrix G :

$$
G = X^{\top} X \tag{1.7}
$$

Given a complete Euclidean Distance Matrix D , the first step of classical MDS is to compute G , starting from D . We can do this by exploiting the Geometric centering matrix J :

$$
J = I - \frac{1}{n} \mathbf{11}^{\top} \qquad (1.8)
$$

$$
G = -\frac{1}{2} J D J \qquad (1.9)
$$

where I is the identity matrix. Schoenberg’s theorem [79] shows that there is a one-to-one relation between a matrix being an Euclidean Distance Matrix and its Gram matrix G being positive semi-definite (all of its eigenvalues are non-negative) [34]. This means that if D is a valid EDM we can use the eigenvalue or singular value decomposition of G to compute X , which is the last step of classical MDS. This method was already proposed in 1901 by the name Principal Component Analysis (PCA) [80], and has been widely used for finding realizations for EDMs [42, 81]. This method works reasonably well, even when the distances include noise. This is because it uses only the top- K eigenvalues to compute the realization matrix X , where K is the dimension we are trying to embed in. However, remember classical MDS only works if we have access to a complete EDM, with all pairwise squared distances. The application that is most relevant to these methods is that of of dimensionality reduction (see Section 1.2), because in this application we often know all pairwise distances between the input data points and can thus construct a complete EDM.

In the case where we have an incomplete EDM, it may still be possible to use EDM- based methods to solve these incomplete instances. The general approach is by first using an algorithm to complete the missing distance information in EDM. Then, as a final step, we can use classical MDS to compute the coordinate matrix X corresponding to a realization of the EDM. One can directly exploit the rank of the EDM. Because the coordinate matrix X has K rows, the ranks of X and X $^{⊤}$X are at most K . The other two terms in (1.6) have rank 1, which means that any valid EDM has a maximal rank of K + 2 [36]. Because we know the embedding dimension K and thus the rank of the
matrix, one approach is to try to iteratively complete the EDM, while enforcing the (noisy) distances that we know, until the desired rank is reached. Alternatively, because we know that the EDM will have a low rank, we may use several general algorithms for low-rank matrix completion that have been proposed in the literature. One example of this is OPTSPACE [82], which is based on spectral methods followed by local manifold optimization. Apart from gradient descent methods, we can also use an alternating least squares minimization approach, such as ALTERN complete proposed in [83].

Aside from exploiting the rank of the EDM, we can try to complete it by propagating triangle and tetrangular inequalities from the known distances. The well-known EMBED algorithm works using this approach [37]. These inequalities give us lower and upper bounds on the missing distances. To be able to use classical MDS, we need an EDM with only real, exact values, and no interval distances. The EMBED algorithm overcomes this issue by iteratively selecting a "trial" matrix by randomly choosing values within the distance intervals. These random values can be chosen uniformly or one can use metrization, to compute a random metric space whose distances all lie within their respective limits [37]. The obtained matrix will most likely not be a valid EDM, which is why the second step includes running local optimization to attempt to convert the invalid EDM to a valid one. A similar approach is the Alternating Projections Algorithm (APA) [84]. This is based on Successive Projection Methodology (SPM) [85], which projects a starting alternatively on two convex sets. Instead of first exploiting triangle and tetrangular inequalities, here we simply generate a random starting point EDM $D$, such that all entries of $D$ either satisfy the known distance constraints or are equal to 0 otherwise. Then, we alternatively project $D$ on two convex sets: the set of matrices with a zero-diagonal and the set of matrices that correspond to a positive semi-definite Gram matrix $G$ (1.9). At the intersection of these two sets lie the completed, valid EDMs.

Semidefinite Programming (SDP) [1, 35, 38, 45, 58, 59] has been widely used for EDM completion as well. Semi-definite programming can be used to solve convex relaxations of a global optimization problem, with the requirement that the variables in the optimization problem correspond to entries in a positive semi-definite matrix. Because we know that the gram matrix $G$ is positive semi-definite, we can cast EDM completion as a semidefinite program, and then we solve this (convex) relaxation. In fact, the EDM completion problem corresponds to a special case of the general Semidefinite Programming feasibility problem [86]. These semidefinite programs may be solved using primal-dual interior point-solvers. However, this is an expensive process and may lead to inexact solutions.
Instead, it is possible to avoid SDP solvers altogether when a facial reduction technique
is used [87, 88]. Interestingly, this method finds solutions by regarding the problem in
three equivalent forms: as a graph realization problem, as the EDM completions problem
as well as a rank restricted SDP. The face of the SDP cone corresponding to a given
clique in the graph is characterized. Next, the intersection of two faces corresponding to
overlapping cliques is characterized. Such an intersection may be obtained by computing
the singular value decomposition. This intersection step is equivalent to completing blocks
of the EDM. This process is repeated until the EDM is fully completed, which can be
done in a finite number of intersection steps. Note that this method cannot be used for
all instances of the problem. In some cases, the entire problem may not be reduced to a
single clique. In this case, it is still necessary to use a SDP solver to complete the EDM.

In general, EDM-based methods can be effective and produce good-quality results. When the EDM is complete, classical MDS is a fast-running method that provides good-quality results even if we have some noise present in the input distances. However, having a complete EDM is a very strict requirement for most important applications. For example, in molecular structure determination, we are missing many distances making the EDM incomplete (see Figure 1.1, bottom-right). When we are missing only a few distances, the EDM completion methods can still lead to good-quality results. Another advantage of some of these EDM completion methods (EMBED and SDP approaches) is that they can handle interval distances in the input as well. However, these algorithms do have some pitfalls. A very important shortcoming is that when the EDM is sparse even the methods based on EDM completion tend to fall short [45]. This means that for many of the harder applications of Distance Geometry, these methods will not be very effective.

Furthermore, when we have instances with a large number of points, these large ma-
trices use a lot of memory space. Furthermore, the MDS based methods do not scale so
well: generally, O(n²) or worse. In the next two sections, we describe (meta) heuristics
and non-linear optimization techniques which aim to address this problem of scalability,
and will generally run much faster on larger applications.

1.4.2 Stochastic Proximity Embedding

One straightforward heuristic has been proposed for the DGP, known as Stochastic Proximity Embedding (SPE) [89]. For a description of the algorithm using pseudocode, see Algorithm 1. We start by initializing with a random initial realization *x* and choosing a learning rate *λ*. From here, we try to improve *x* in iterations, until a maximum number
of iterations is reached. To do this, we have a double loop: the inner loop runs for several steps S while the outer loop runs for an amount of cycles C. In the inner loop, at every step, we select one random vertex pair (u, v) for which we know the true distance $\delta(u, v)$. Then, we adjust the coordinates $x_u$ and $x_v$ such that the distance $||x_u - x_v||$ better approximates $\delta(u, v)$. The amount by which we "correct" the coordinates is calculated by using the learning rate $\lambda$. In the outer loop, after each set of steps, we decrease our $\lambda$ by a predetermined amount. This means that in the earlier cycles, we will update the coordinates more drastically. When we start to converge towards a solution, only small modifications are made.

**Algorithm 1 Stochastic Proximity Embedding (SPE)**

**Input:** Cycles C, steps S, learning rate $\lambda$ and decrease factor $\delta_\lambda$, small tolerance value $\epsilon$

**Output:** Realization $x$

$x \leftarrow$ initial random realization

while $c \le C$ do

while $s \le S$ do

Select a random distance $\delta(u, v)$

$d_{uv} \leftarrow ||x_u - x_v||$

if $d_{uv} \in \delta(u, v)$ then

$r_{uv} \leftarrow \delta(u, v)$ if it is exact, otherwise the nearest bound to $d_{uv}$

$x_u \leftarrow x_u + \frac{\lambda r_{uv} - d_{uv}}{2(d_{uv} + \epsilon)} (x_u - x_v)$

$x_v \leftarrow x_v + \frac{\lambda r_{uv} - d_{uv}}{2(d_{uv} + \epsilon)} (x_v - x_u)$

end if

$s \leftarrow s + 1$

end while

$c \leftarrow c + 1$

$\lambda \leftarrow \lambda - \delta_\lambda$

end while

return $x$

SPE was implemented as part of this thesis to serve as a benchmark for the other solvers (see Section 1.5). One of the main advantages of SPE is that it is very straightforward to implement. SPE is a great alternative to EDM based methods, as the computation time scales in order of $O(n)$ and because it can trivially handle interval distances. Furthermore, the inner loop of the SPE procedure lends itself to parallelization, which will make it run even faster. Another advantage of SPE over the MDS techniques is that it does not require the distance graph to be dense at all. However, these advantages come at a price. As it is such a simple heuristic, we cannot always expect high-quality solutions.
Because of its stochastic nature, we have no guarantees on the accuracy of the solution.
In the next section, we will see some (global) optimization techniques that are a little
more robust than SPE.

1.4.3 Nonlinear Programming

In a heuristic like SPE, at every iteration, we only modify up to two positions of the
realization *x*. However, we can also update the positions *x*<sub>*v*</sub> for all *v* ∈ V, at the same
time, by regarding the problem as a Global Optimization (GO) problem (see Definition 2).
We may use different Nonlinear Programming (NLP) methods to attempt to solve this
global optimization problem. These methods exploit the fact that penalty functions such
as (1.4) are differentiable on realizations *x* when there are no pairs of vertices that have
the same position [26].

When we use the graph representation of the DGP (see Definition 1) we may use some
mixed combinatorial methods based on graph decomposition. The idea is to divide our
graph *G* into a set of smaller graphs. Then, in the first phase, we use local optimization
to solve these subgraphs. In the second phase, we attempt to combine the solutions to
the subgraphs to find a solution for the global optimization problem. The first algorithm
to use this graph decomposition approach was ABBIE [90], where local optimization
techniques were used for both phases. Alternatively, Semidefinite Programming techniques
may be used for both phases as well [91, 92]. This second, global optimization phase of this
approach is the most challenging step and is what can complicate matters. Often, heuristic
approaches must be taken to do this effectively. An example of this is the algorithm 3D-
ASAP [93]. A drawback of these graph decomposition methods is that they only work
on exact distances, and extending them to interval data is not trivial [22].

Another option is Spatial Branch-and-Bound (SBB), which is a divide-and-conquer
technique used to find the solution to global optimization problems. SBB uses a deter-
ministic search using three steps: branching, bounding and, pruning. We start with the
branching phase, in which we divide the global problem into regions. Each of these regions
lets us compute lower and upper bounds on the objective function by utilizing local op-
timization. This is called the bounding step. Finally, in the pruning step, we discard any
regions that have lower bounds larger than the best local optimum found so far. We repeat
these three steps, on increasingly smaller regions, until the computed bounds converge,
which means we have found a local optimum. This general purpose approach was tried
for the DGP in [94]. The advantage is that it is deterministic, while most of the other
approaches discussed in this section are stochastic. However, this method is extremely
slow, even for small-sized instances.

There are also several algorithms based on smoothing, in which local optimization is used to find minimizers for increasingly finer smoothings of the original global optimization problem. At each iteration, we end up with fewer local minima, until we finally end up with a convex optimization problem that we can solve with local optimization. This is known as a *global continuation* method. A well-known example of this approach is the DGSOL algorithm [43, 44]. A similar algorithm (DCA) was introduced in [95, 96], which differs mostly in their choice of local optimizer. DGSOL is a fast-running algorithm that provides good-quality solutions for small instances. To extend to (slightly) larger instances, DGSOL was combined with the meta-heuristic VNS [97] (discussed as a stand-alone method in the next section). This lead to better quality solutions for larger instances, but increased the computation time of the approach.

An example of a global constrained optimization method that was introduced specifi-
cally to solve instances related to molecular structure determination is GNOMAD [98].
This is a multi-level iterative algorithm that groups the vertices at the highest level. Next,
it attempts to identify an order within each of the groups based on the contribution of
each vertex to the total error on the distance constraints. Then, in this order, the vertices
are moved using a series of local optimization steps, while making sure all constraints
remain satisfied. This algorithm can handle exact distance constraints as well as intervals.

These non-linear optimization methods can be a good compromise between the ac-
curacy of MDS methods and the fast computation times of SPE. Not all methods work
well with interval data, such as the graph decomposition methods. On the other hand,
methods based on smoothing can be adapted such that they use *hyperbolic smoothing*,
which works for interval data [99]. In general, we cannot guarantee the level of quality
of the solutions. However, SBB is an example that does provide us with a bound on the
expected error. Despite this, the SBB method is not the best option, as it is extremely
slow on larger instances [94].

1.4.4 Meta-heuristics

A common way to tackle the optimization problem given in Definition 2 is by means
of *meta-heuristic searches*. The general idea of these methods is to repeatedly find local
minima for the penalty function *σ*, until something close to the global minimum has been
found. Different meta-heuristics have been tried in the context of Distance Geometry. Two
that were tested in [94] were variations of Multi-start [100] and Variable Neighbourhood Search (VNS) [101]. The most simple of these is Multi-start, which works by repeated local searches (gradient descent) from randomly selected starting points. Multi-start was tested on the DGP in [102]. VNS is a more sophisticated meta-heuristic. Instead of running gradient descent from random starting points, it performs the searches in neighbourhoods of already identified local optima, in an attempt to find the global optimum. These two methods were effective for finding solutions for the DGP, but only for instances on the small side, with *n* < 50 [94]. For larger instances, the method is too slow.

Two meta-heuristic global optimization methods for the DGP were presented in [102]. These methods are variations of the Basin Hopping (BH) algorithm, which was inspired by Monte Carlo minimization and used for finding minimum energy structure for molecules [103]. For solving instances of the DGP, two versions of this algorithm were tried: Monotonic Basin Hopping (MBH) and Population Basin Hopping (PBH). MBH looks like VNS in the way that it repeatedly finds local minima and performs neighbourhood changes. However, in contrast to VNS, a change of neighbourhoods is done by performing a *local move*, which is a procedure that allows us to move from one local minimum to another. The nature of these local moves is very important to the efficacy of the algorithm. Because these algorithms are stochastic, a better choice for the local move procedure will lead to a higher chance of finding a global optimum. Two options for local moves are discussed in [102], showing that a problem-specific local move is the best option. PBH is an extension to MBH that includes performing several MBH runs in parallel, in order to increase the chance of success. Aside from these algorithms, there are a group of methods that are based on the Simulated Annealing (SA) meta-heuristic. These are specialized for protein structure determination and are discussed in Section 3.4.1.

The stochastic meta-heuristics find good-quality solutions for small instance sizes [22, 94, 102]. However, the larger the instances get, the more difficult it becomes for some of these methods to perform well. Just like the optimization techniques mentioned in the previous section, most of these methods do not give any guarantee about the level of the quality of the solutions. They are perhaps best used as a refinement tool, in combination with another DG method. In general, the meta-heuristic methods are quite easily adapted such that they work with interval distances as well as exact distances. These methods differ also in another facet: some of them are general purpose while others need are more specialized. The advantage of the general purpose methods is that generally there is a lot of information about them in the literature, and often their implementation is
straightforward. The specialized methods may be harder to implement, but in the end,
they may provide us with better results, especially if we further fine-tune them to the
specific application that corresponds to the instance of the DGP that we are attempting
to solve.

1.4.5 Spectral Projected Gradient

We focus our attention on a specific local optimization method: Spectral Project Gradient (SPG). SPG is often used in the context of Distance Geometry [4, 39]. We highlight this method because it was used in various experiments in this thesis, across different applications. In particular, the algorithm can be used to refine a solution obtained via another approach. For example, it has shown to be an effective refining tool for instances relating to molecular structure determination, where it was used in combination with SDP and classical MDS [104]. SPG can be used for constrained as well as unconstrained optimization. In the Java package (see Section 1.5) for this thesis, we used the constrained variant in combination with branch-and-prune (see Section 1.4.7 and Chapter 2), which has been shown to be a fruitful approach in [9]. Aside from this, we will show experiments in the context of Dynamical Distance Geometry Problem, where the unconstrained variant of SPG was used as the main method for solving the motion retargeting problem. SPG was used as a tool for the motion retargeting problem also in [4, 19, 105, 106].

Given an initial realization $x_0$, the spectral projected gradient is a method that will
try to improve the quality of the solution iteratively by minimizing the penalty function
$\sigma$. This initial realization $x_0$ can be chosen arbitrarily, but when we use SPG as a refine-
ment method, $x_0$ will generally already be of a decent quality. The main component of
SPG is gradient descent, which is a well-known local optimization procedure and one of
the principal ingredients of deep learning methods such as neural networks [107–110]. For
gradient descent, the step size is an important parameter. If the step size is too large,
we might step over a minimum. When it is too small we could take forever to converge.
We attempt to address this issue by using a spectral projected gradient method, which
means that the step size used throughout the gradient descent is spectral [111], based on
the underlying local Hessian, rather than using a standard decrease in the objective func-
tion. The spectral projected method is combined with a non-monotone line search [112] to
speed up the convergence. For this line search method, global convergence was shown for
smooth and non-convex functions, and R-linear convergence was proven for strongly con-
vex functions [113]. It is important to note that these convergence properties were shown
for unconstrained optimization problems. In practice, the method is still very effective for constrained optimization.

The pseudocode for the full procedure is given in Algorithm 2 and will be referenced throughout the explanation. This pseudocode is based on [106, 111] and the implementations of SPG in the context of Distance Geometry. SPG operates through a main loop, which runs for a maximal number of iterations *I*, which is given as an input parameter.

<caption>**Algorithm 2** Spectral Project Gradient (SPG)</caption>

**Input:** Initial solution $x_0$, penalty function $\sigma$, maximal iterations *I*, thresholds $\epsilon$ and $\alpha_\epsilon$, stepsize bounds $\mu_{\min}, \mu_{\max}$, decrease parameter $\gamma$, nonmonotonicity degree $\eta$
**Output:** Final solution $x$
1: $i \leftarrow 0, q_0 \leftarrow 0$
2: $m_0 \leftarrow \sigma(x_0)$
3: **while** $i < I$ and $\sigma(x_i) > \epsilon$ **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;**if** $i = 0$ **then**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mu_0 \leftarrow 1$
6: &nbsp;&nbsp;&nbsp;&nbsp;**else**
7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$y_{i-1} \leftarrow \sigma(x_i) - \nabla\sigma(x_{i-1})$
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$s_{i-1} \leftarrow x_i - x_{i-1}$
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\mu_i \leftarrow (s_{i-1}^\top s_{i-1})/(s_{i-1}^\top y_{i-1})$
10: &nbsp;&nbsp;&nbsp;&nbsp;$\mu_i \leftarrow \min(\mu_{\max}, \max(\mu_{\min}, \mu_i))$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Using safeguards
11: &nbsp;&nbsp;&nbsp;&nbsp;**end if**
12: $d_i \leftarrow -\mu_i \nabla\sigma(x_0)$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Using $\mu_i$ to compute $d_i$
13: $d_i = P_\Omega(x_i - d_i) - x_i$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Projection constraints (only for constrained version)
14: $\alpha_i \leftarrow 1$
15: **while** $\sigma(x_i + \alpha_i d_i) > m_i + \gamma \alpha_i \nabla\sigma(x_i)^\top d_i$ **do** &nbsp;&nbsp;&nbsp;&nbsp; ▷ Nonmonotone line search
16: &nbsp;&nbsp;&nbsp;&nbsp;$\alpha_i \leftarrow \alpha_i / 2$
17: **end while**
18: **if** $\alpha_i < \epsilon_\alpha$ **then**
19: &nbsp;&nbsp;&nbsp;&nbsp;**return** $x_i$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Terminate early (stopping criteria 3)
20: **end if**
21: $x_{i+1} \leftarrow x_i + \alpha_i d_i$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Computing next solution and preparing next iteration
22: $q_{i+1} \leftarrow \eta q_i + 1$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ $q_i$ is a helper variable for $m_i$
23: $m_{i+1} \leftarrow (\eta q_i m_i + \sigma(x_{i+1})) / q_{i+1}$
24: $i \leftarrow i + 1$
25: **end while**
26: **return** $x_i$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Output the final solution

Recall the penalty function $\sigma(x)$ which can be used to measure the violation of the
distance constraints for a realization $x$:

$$ \sigma(x) = \frac{1}{2} \sum_{\{u,v\} \in E} (\|x_u - x_v\| - \delta(u, v))^2 $$

Because this function is differentiable (when no two points have the same position), we can use its gradient to improve a solution $x$. At every iteration $i \in 0, \dots, I$, a new, improved solution $x_{i+1}$ is computed by performing a step in the opposite direction $d_i$ of the gradient $\nabla\sigma(x_i)$. The direction $d_i$ can be computed from the gradient using the spectral stepsize $\mu_i$ (line 12) and the projection $P_\Omega$ [114, 115] (line 13):

$$ d_i = P_\Omega(x_i - \mu_i \nabla\sigma(x_i)) - x_i, $$

$\Omega$ is the feasible convex set onto which we know how to project solutions $x$. Incorporating the projections $P_\Omega$ is what allows us to perform constrained optimization by using SPG. The spectral stepsize $\mu_i$ is computed based on the previous step from iteration $i-1$ to $i$ (line 9):

$$ \mu_i = \frac{s_{i-1}^\top s_{i-1}}{s_{i-1}^\top y_{i-1}}. $$

where $s_{i-1} = x_i - x_{i-1}$ and $y_{i-1} = \nabla\sigma(x_i) - \nabla\sigma(x_{i-1})$. This stepsize should be bounded by the parameters $\mu_{\min}$ and $\mu_{\max}$ as a safeguard (line 10). The last ingredient we need is the multiplier $\alpha$, obtained from the non-monotone line search (lines 14-17). Initially, $\alpha_i$ is set to 1.0. Next, $\alpha_i$ is reduced by half until the following stopping criteria are met (line 15):

$$ \sigma(x_i + \alpha_i d_i) \leq m_i + \gamma \alpha_i \nabla\sigma(x_i)^\top d_i, $$

where $\gamma \in (0, 1)$ is a sufficient decrease parameter and $m_i$ is a variable based on the penalty values $\sigma$ of the previous iterations. Often, a parameter $\eta$ is included in the line search as well, which controls the degree of nonmonotonicity and is used to update $m_i$ throughout the iterations. More information about these parameters and variables as well as an analysis of this step can be found in [113]. Once we have obtained the three ingredients $\mu_i$, $d_i$ and $\alpha_i$, we proceed to compute the next solution $x_{i+1}$ by taking a step in the right direction (line 21):

$$ x_{i+1} = x_i + \alpha_i d_i. $$

The above three steps (computation of spectral stepsize, non-monotone line search, and stepping in the opposite direction of the gradient) are repeated until one of several
stopping criteria is reached:

1. The maximal number of iterations *I* is reached (line 3).

2. The objective value $\sigma(x_i)$ reaches below a certain threshold $\epsilon$, where $\epsilon$ is given as a parameter (line 3).

3. The stepsize multiplier $\alpha$ reaches below a threshold $\alpha_{\epsilon}$, where $\alpha_{\epsilon}$ is given as a parameter (line 18).

### 1.4.6 Distance norms and Linear Programming

The optimization variant of the DGP in Definition 2 is non-linear. Firstly, the penalty function $\sigma$ has quadratic terms. Next, the Euclidean distance norm ($L_2$), which we generally have in the applications of the DGP, is non-linear:

$$L_2(x_u, x_v) = \sqrt{\sum_{i}^{K} (x_u^i - x_v^i)^2} \qquad (1.10)$$

where $x_u^i$ is the $i$-th coordinate of vertex $u$ in the realization $x$.

However, it is possible to linearize this problem when we replace the Euclidean norm ($L_2$) with a linearizable distance measure. In [31], it was proposed to replace the Euclidean distances by the $L_{\infty}$ norm:

$$L_{\infty}(x_u, x_v) = \max_{i} |x_u^i - x_v^i| \qquad (1.11)$$

This norm has some useful properties when the DGP instance contains interval distances. When we use the $L_{\infty}$ norm, these bounds are now linear. In [32] this idea was expanded upon. Here, not only the $L_{\infty}$ norm was considered, but also the $L_1$ norm, which is also known as taxi-cab or Manhattan distance:

$$L_1(x_u, x_v) = \sum_{i}^{K} |x_u^i - x_v^i| \qquad (1.12)$$

These alternative distance norms cannot be linearized without the introduction of binary variables. Using these binary variables, it is possible to replace nonlinear functions with piecewise linear forms. This allows us to transform the non-linear optimization problem described in Definition 2 to a linear one. We can specify Mixed-integer programming (MIP) formulations of the now linear problem, which is a special type of linear program
(LP) which allows for a mix of binary and integral variables. The MIP formulations for
the $L_1$ and $L_\infty$ norms can be found in [32].

In practice, using Linear Programming for general instances of the DGP is not so
effective as it becomes very slow with larger input sizes, due to the combinatorial explosion
of constraints [32]. This approach is most promising in the context of the additional
local orientation constraints discussed in Section 1.3. For the application of adaptive
maps, this approach provides us with good-quality solutions with a low computation
time. Experiments with this approach are presented in Chapter 5.

1.4.7 Build-Up and Branch-and-prune

In this section, we describe two unique approaches. The methods discussed here do not
rely on solving the DGP as an optimization problem or matrix decomposition problem.
Instead, these methods iteratively embed the vertices by exploiting geometric properties
of the combinations of distances that we have access to.

The BUILD-UP algorithm was proposed in the context of the DGP for molecular
structure determination [116–118]. As the name of the algorithm suggests, we attempt
to build up an embedding for the vertices V in a sequential manner. Because the focus
was on determining protein structures, the dimension K is set to 3. We start by fixing
the coordinates for the first four vertices. Next, we find a vertex v for which we know a
distance to at least four vertices that we have already embedded. These four vertices are
called reference vertices for v. When we have a vertex v and four such reference vertices, we
can use trilateration to compute the coordinates for v. We repeat this process until there
are no vertices left to embed or we cannot embed any more vertices (in which case the
algorithm fails). The original version of the algorithm had an issue with an accumulation
of round-up errors, which meant that there were large inaccuracies in the coordinates
of the vertices embedded towards the end. An update to the algorithm was made [118]
such that these numerical errors were put under control. A second modification to the
algorithm made it so only three reference vertices are needed when embedding a vertex
v [119]. It is important to note that this approach only works if we have the right number
of distances available at the right time. Essentially, we require a specific order on the
vertices. Because we do not know beforehand which vertex order will lead to a successful
run of the algorithm, this approach is heuristic in nature.

This lead to the development of the branch-and-prune (BP) framework [10]. Similar
to the BUILD-UP algorithm, we iteratively embed the vertices. This framework works for
all dimensions *K*, as long as when embedding a vertex *v* we have at least *K* reference
distances for *v*. An important difference with the BUILD-UP algorithm is that in the BP
framework, much research has been done on how to identify vertex orders which guarantee
that every vertex *v* will have a sufficient number of reference distances. This means that,
given the correct vertex order, the BP algorithms are not heuristics, but exact methods
(as long as we are dealing with exact distances). Furthermore, instead of identifying one
position at every step, we may identify two or more possible positions for a vertex *v*.
Effectively, this allows us to discretize the search space, making it possible for us to
enumerate every possible solution to an instance of the DGP that allows for a vertex
order that satisfies our requirements. An in-depth description of the branch-and-prune
algorithm and its implementations is given in Chapter 2.

The several papers that researched the BUILD-UP algorithm showed that can be very
effective for determining protein structures given somewhat sparse distance graphs [117,
118]. The resulting solutions had low Root mean square deviation (RMSD) scores when
compared to the structures obtained from x-ray crystallography. The biggest drawback
of the BUILD-UP algorithm is that we are not guaranteed to find a solution. When we
embed the vertices using a greedy method, we do not know if we will always have the right
reference distances available. A second important drawback, when comparing BUILD-UP
to BP, is that the requirements for a vertex order for the BUILD-UP algorithm are stricter
than for BP.

Another advantage of the BP framework is that it has several extensions that allow us
to solve instances that contain interval distance data [11, 12]. The BP framework is one of
the most promising avenues of research for solving the DGP, especially in the context of
structural biology. This framework was implemented in the Java package and intensively
used for experiments in this thesis.

1.5 The Java Package

Several different methods for the DGP were implemented as part of this thesis. Most
of these methods were grouped into one Java package. While this Java package is not
published in its entirety (as it is still in development), a subset of the classes used in
the context of protein structure determination was made public on GitHub<fnref n="1" />. More in-
formation about this repository is found in Section 3.5.1. The class diagram of the full

<footnote n="1">1. https://github.com/simonheng/BP_ProteinFileReader</footnote>
Java package is shown in Figure 1.4. This diagram shows the relations between the several classes in the package based on dependency and inheritance. To avoid the diagram from getting cluttered, duplicate dependencies are not visualized. For example, the utility classes MATHS and GEOMETRY are used by almost all other classes, but these arrows are not included if the dependency is already passed through a different class.

<figure><img src="image_8.png" /><figcaption>Figure 1.4 – A class diagram representing the java package that was created for this thesis.</figcaption></figure>

In the top left of the diagram we have four classes with subroutines that are used
by almost all other classes in the project. PAIR is an object class which represents an
immutable tuple, while the MATHS and GEOMETRY are utility classes containing static
methods corresponding to subroutines such as complex number division and sphere in-
tersections (Section 2.3). MATRIX is an object class for matrices, based on the JAMA
package<fnref n="2" />.

We use the graph representation of the DGP, which is why we include the object
classes GRAPH and WEIGHT to describe weighted directed and undirected graphs. The
DGP object class is an extension of the GRAPH class where the weights are DISTANCE

<footnote n="2">2. https://math.nist.gov/javanumerics/jama/</footnote>
objects. The vertices of the graph described by the DGP class are extensions of the ab-
stract COORDINATES class. For each of the different applications of the DGP, the package
contains an object class that extends the DGP class which represents the instances for
the application in question. For protein structure determination (Chapter 3) we have the
class MDGP, which stands for Molecular DGP. The vertices of the graph that MDGP
describes are of the ATOM type, which is a special implementation of the CARTESIAN
coordinate subclass. MDGP also relies on three other classes which deal with molecu-
lar input data: FORCEFIELD, NMR and PDB. More information about this is given in
(Chapter 3). For motion retargeting (Chapter 4), we have the BVH class, which uses the
EULER coordinate extension class for its vertices. The VECTOR TORSION object class
represents the new representation that we propose for human motions, described in Sec-
tion 4.4. Finally, for adaptive maps (Chapter 5) we have the ADAPTIVEMAP class which
uses CARTESIAN vertices, where the dimension K is 2.

The different methods of solving the DGP are included in the Java package as imple-
mentations of the SOLVER interface. The solutions to the problems may be represented
using the REALIZATION class. Some of these solvers work for any extension of the DGP
class while other work for only specific ones. The MDJEEP class contains an implemen-
tation of the branch-and-prune framework (Chapter 2) which only works if the DGP is
discretizable. As we will see in Chapter 3, this is the case for the instances that arise
for protein structure determination, which is why in the Java package MDJEEP is used
specifically to solve instances represented by the MDGP class. The class REVORDER is
used to make sure that we have an order of vertices which allows us to discretize the
search space (see Section 2.1 and Definition 5). As we will see in Chapter 5, we use Linear
Programming in the context of Adaptive Maps, which is why the LINEARSOLVER class
works specifically with the ADAPTIVE MAP class. The LINEARSOLVER class relies on
the Java interface provided by CPLEX [120] using the Concert Technology. The CLASSI-
CAL MDS class works on DGP objects as long as the coordinates are CARTESIAN and
the graph is complete. The classical MDS (Section 1.4.1) method was implemented as a
benchmark for the linear programming technique used for adaptive maps. As a benchmark
for the more general methods, SPE (Section 1.4.2) was implemented in the SPE class,
which is compatible with any of the DGP (sub)classes that use CARTESIAN coordinates.
The implementation of SPG (Section 1.4.5) works with any extension of the DGP object.
When the coordinates used are CARTESIAN, SPG uses an exact method to compute the
gradient. This is the case when SPG is used as a refinement step for MDJEEP (more
on this in Section 2.4.2). In the context of motion retargeting (Chapter 4), where the coordinates are *EULER* (or perhaps in the future *VECTOR TORSION*), SPG relies on a numerical differentiation method (finite differences) [121] to approximate the gradient.
# DISCRETIZABLE DISTANCE GEOMETRY

Given certain assumptions, the search space relating to an instance of the Distance Geometry Problem can be *discretized*. In fact, we can construct a tree in which each node corresponds to a candidate position for a vertex *v*. Every path from the root node to a leaf node corresponds to a possible solution to the DGP instance at hand. Combinatorial methods may be used to explore the tree and enumerate every possible solution for the given DGP instance. In particular, we will look at the branch-and-prune framework and its implementations for the dimension *K* = 3. As part of this thesis, one of the two open-source branch-and-prune implementations was implemented into Java with some extra features.

In this chapter, we begin with the main theory of discretization and vertex orders
in Section 2.1. Next, in Section 2.2, we describe the algorithmic framework of branch-
and-prune, which is the method that exploits the discretization. Section 2.3 discusses
several methods that can be used for generating coordinates during the branching step of
BP. Afterwards, in Section 2.4 we investigate how the different implementations handle
interval distances in the input. Finally, we discuss how including torsion angles in the
DDGP instances lets us decrease the size of the search space (Section 2.5).
<table><tr><td>2.1</td><td>Introduction</td><td>54</td></tr><tr><td>2.2</td><td>Branch-and-prune in dimension 3</td><td>57</td></tr><tr><td>2.3</td><td>Coordinate generation</td><td>60</td></tr><tr><td>2.4</td><td>Handling interval distances</td><td>63</td></tr><tr><td>2.4.1</td><td>Sampling the arcs</td><td>63</td></tr><tr><td>2.4.2</td><td>A coarse-grained representation and local optimization</td><td>64</td></tr><tr><td>2.5</td><td>Torsion angles and the search space</td><td>65</td></tr></table>
## 2.1 Introduction

It has been shown that, when several assumptions are met, a DGP instance $G = (V, E, d)$ with $K > 0$ can be discretized [6]. The discretization of the search space of such an instance makes it possible for combinatorial methods to enumerate all feasible solutions for the problem at hand. This is of great interest, especially in the applications of molecular structure determination (see Chapter 3), because it allows us to identify distinct protein structures given inter-atomic distances. In fact, in the context of protein structure determination, the instances often satisfy the required assumptions for discretization. For more information about which inter-atomic distances are required to discretize, see Sections 3.3 and 3.4.3. The subclass of the DGP which captures such instances is referred to as the Discretizable Distance Geometry Problem (DDGP). The DDGP was shown to be NP-hard [6, 8].

As part of this thesis, an explanatory video concerning this subclass of instances was created and entered to the Summer of Math Exposition (SoME) 22 contest. It can be found on Youtube<fnref n="1" />. In order to discuss this subclass and the necessary discretization assumptions, we need to partition the set of edges $E$. Let $E'$ be the subset of edges in $G$ that are related to exact distances. As a consequence, $E \setminus E'$ is the subset of edges in $G$ that contain all the interval distances. We can now give a formal definition of the DDGP [6, 26]:

**Definition 4** A simple weighted undirected graph $G = (V, E, d)$ represents a DDGP instance in dimension $K$ if and only if there exists a vertex ordering on $V$ such that the following two assumptions are satisfied:

(a) $G[\{1, 2, \dots, K\}]$ is a clique whose edges are in $E'$;

(b) $\forall v \in \{K+1, \dots, |V|\}$, there exist $K$ vertices $u_1, u_2, \dots, u_K \in V$ s.t.

(b.1) $u_1 < v, u_2 < v, \dots, u_K < v$;

(b.2) $\{\{u_1, v\}, \{u_2, v\}, \dots, \{u_{K-1}, v\}\} \subset E'$ and $\{u_K, v\} \in E$;

(b.3) $\mathcal{V}_S(u_1, u_2, \dots, u_K) > 0$ (if $K > 1$),

where $G[\cdot]$ is the subgraph induced by the subset of vertices of $V$ given in argument and $\mathcal{V}_S(\cdot)$ is the volume of the simplex $S$ generated by a valid realization of the vertices $u_1, u_2, \dots, u_K$.

<footnote n="1">1. https://www.youtube.com/watch?v=ZTR8txn2wBU</footnote>
The first assumption (**a**) lets us fix the coordinate space of the instance at hand. Because of assumptions (**b.1**) and (**b.2**), at least *K* other distinct vertices, preceding the current vertex *v* in the given vertex ordering, can be used as a reference for positioning the vertex *v*. Because of this notion, these vertices are referred to as *reference vertices* and the distances $\delta(u_1, v), \delta(u_2, v), \dots, \delta(u_K, v)$ are referred to as *reference distances*. Moreover, assumption (**b.2**) makes sure that at most one of the *K* available distances has a large enough degree of uncertainty to be considered as an interval distance, while at least *K* − 1 other distances can be considered as “exact”. These two previous conditions ensure that for every vertex $v \in V$ to be placed, there is only a “small” subset (under some conditions, a discrete subset) of feasible positions for that vertex. This is because the *K* reference distances can be used to define *K* Euclidean objects, which we can then intersect to obtain a discrete number of positions for *v*. For each reference vertex $u_i$ for which $\{u_i, v\} \in E'$, this object is a *K*-sphere, which we construct with radius $\delta(u_i, v)$. If one of the reference vertices ($u_K$) relates to an interval distance, we will instead construct a *K*-spherical shell, where the thickness of the shell is linked to the size of the interval. The assumption (**b.2**) is necessary because if more than one of the reference vertices would relate to an interval distance, the intersection between the Euclidean objects would not correspond to a disjoint set of arcs or points which would mean that we cannot use the algorithms that will be discussed in this chapter. The final assumption (**b.3**) makes sure that these *K* Euclidean objects will have at most two disjoint sets as intersection. In the case that all *K* reference distances are exact, these *K*-spheres will intersect in at most two positions where we know that all reference distances are satisfied. These are the candidate positions for *v*. When we instead intersect *K* − 1 spheres and 1 spherical shell, we obtain at most two disjoint arcs. The vertex *v* may be anywhere inside these two arcs, which means the problem is still continuous. These arcs may then be sampled so that we can still obtain a discrete number of positions for *v*. More details on this will be given in Section 2.4. This discretization lets us model the domain of the search space as a tree, where we have a layer for every vertex *v*. Each possible position of *v* will correspond to a node in this layer. An example of such a tree is given in Figure 2.1. Note that the first *K* vertices (that form a clique in *G*) are fixed, so we have only one candidate position for each of them (modulo isometries).

To be able to satisfy the assumptions in Definition 4 we need a specific order on vertices. Given an unordered instance of the DGP and reordering it such that the assumptions are met is not a trivial task. For instances related to molecular structure determination, different handcrafted orders have been proposed in the literature [11, 122, 123]. Reorder-
<figure><img src="image_9.png" /><figcaption>Figure 2.1 – An example of a tree representing the search space of an instance of the DDGP with $|V| = 6$ and $K = 3$. A path in the tree, corresponding to a possible solution to the instance is highlighted.</figcaption></figure>

ing the vertices of $V$ such that the discretization assumptions are met corresponds to a problem in the literature known as the Referenced Vertex Ordering (Revorder) [124] problem. Let $\sigma$ be a vertex order with ranks $1, \dots, n$, where $n = |V|$ and let $R_{\sigma}(i)$ be the reference set of the vertex $v$ with rank $i$ in the order $\sigma$. The Revorder problem can be defined as a decision problem as follows [124]:

**Definition 5** Given a simple directed graph $G = (V, E)$, a positive integer $K$ is there a vertex order $\sigma$ such that:

$$ \begin{cases} \forall i \in \{1, \dots, K\}, & |R_{\sigma}(i)| = i - 1 \\ \forall i \in \{K + 1, \dots, n\}, & |R_{\sigma}(i)| \ge K \end{cases} $$

A vertex order that satisfies the above definition is called a referenced order. Note that this definition does not distinguish between reference vertices relating to exact and interval distances. However, the definitions and the algorithms discussed in this section can be trivially extended so that the assumptions from the previous section are precisely met.

In the literature, several works have been published which study a problem that is similiar to the Revorder problem [8, 125]. This problem is known as the Discretization Vertex Order Problem (DVOP) [125] and is in fact a specific case of the Revorder problem where an extra assumption is included. This assumption is known as the consecutivity assumption, which requires that the reference set of a vertex $v$, together with $v$ itself, are consecutive in the vertex order. However, this is not strictly required for discretization, and thus in later works this additional assumption was dropped.
For the DVOP, a greedy algorithm was proposed in [8]. This algorithm was greatly improved upon in [124], which extended it to the more general Revorder problem. The greedy search starts from an initial clique (of size $K$) and attempts to build a referenced order from there. The running time of this greedy algorithm is $O(n \log(n) + m)$, where $n = |V|$ and $m = |E|$. The same paper also proposed other algorithms and branch-and-bound enumeration schemes that pertain to the Revorder problem. The greedy algorithm was implemented into Java during this thesis, and used for various experiments (the REVORDER class in Figure 1.4).

Once we have identified a suitable order for the vertices, either by using an algorithm or by following a handcrafted order, we may proceed with exploring the full search space of the tree. This can be done by utilizing the branch-and-prune algorithm, which is discussed in the next section.

## 2.2 Branch-and-prune in dimension 3

The branch-and-prune (BP) [10–12] algorithm is the method that exploits the discretization of the DDGP instances. We will focus on the case where $K = 3$, which matches the application of structural biology. In this section, we will describe the general branch-and-prune framework. In the next sections, we will describe how each phase of the algorithm may be handled or implemented differently.

The algorithm operates in three main phases, the *initialization phase*, the *branching phase*, and the *pruning phase*. In the initialization phase, we embed the first three vertices $u, v$, and $w$ in the order, which is possible because these three vertices are a clique in the graph (see assumption (a) in Definition 4). For example, we may fix the first vertex directly at the origin: $x_u = (0, 0, 0)$. We can then use the distance $\delta(u, v)$ to place the second vertex, e.g. $x_v = (-\delta(u, v), 0, 0)$. To place the third vertex $w$, we draw two circles. One around $u$, using as radius $\delta(u, w)$ and one around $v$ using radius $\delta(v, w)$. These circles intersect in two points, which are both candidate positions for $w$. We can simply pick one of the two and use it to place the third vertex $w$ (the other intersection point leads to a symmetric solution). Essentially, in this phase, we are fixing the coordinate space.

After the initialization phase, the next step is to iteratively or recursively embed the rest of the vertices, exploiting the vertex order. We do this by switching between the branching phase and the pruning phase until we have identified all valid candidate positions for every vertex $v \in V$.
In the branching phase, we perform the sphere intersections to compute the possible
embeddings of a vertex $v$. This step is illustrated in Figure 2.2. Different methods can be
used to compute these candidate positions, and this differs among the implementations.
Section 2.3 will discuss these various methods. In the case that one of the reference
distances is an interval, several approaches can be taken so that we may still discretize
the search space. This will be further investigated in Section 2.4. It is important to note
that when all distances are exact, BP is an exact method that guarantees to enumerate
all possible solutions to the DGP instances at hand. However, when we have interval
distances, the sampling methods that the implementations use are heuristic, and thus the
algorithm loses this property of exactness.

<figure><img src="image_10.png" /><figcaption>Figure 2.2 – An illustration how a vertex $v$ is embedded in dimension $K = 3$, using three reference distances. The distance $\delta(u_3, v)$ may be an interval distance.</figcaption></figure>

Once we have obtained a candidate position for a vertex $v$ in the branching phase, we validate it during the pruning phase. In this phase, we verify if there exists any other vertex $u < v$ for which we have a known distance $\delta(u, v)$ and that is not one of the reference vertices. For each of such vertices $u$, we check whether the computed position for $v$ satisfies the distance $\delta(u, v)$. If not, we reject this candidate position and prune the corresponding branch. The distances $\delta(u, v)$ are referred to as *pruning distances*. When a pruning distance $\delta(u, v)$ is checked, the implementations generally use a tolerance parameter $\epsilon$, allowing for a small error. The reason for this is to take into account noise that may be included in the input distances. Other information, aside from distances, may be verified as well in this phase. In [126], experiments were conducted in the context of molecular structure determination, where different pruning devices were used. One of the pruning devices
that was tested was based on verifying torsion angles $\tau$, that were part of the instances, as we described in Section 1.3. However, as we will see in Section 2.5, it may be more advantageous to directly use the torsion angle information during the branching phase so that we can avoid creating infeasible branches.

Two open-source BP implementations for $K = 3$, which are both written in C, are publicly available: IBP-NG<fnref n="2" /> [127] and MDJEEP<fnref n="3" /> [7]. A multi-threaded implementation of IBP-NG is also available, known as TAIBP<fnref n="4" /> [128]. MDJEEP was incorporated in the Java package for this thesis (Figure 1.4), so that several new features could be tested. While, compared to C, Java is slower at execution time, it is easier to permute which means different ideas can be implemented and tested at a higher rate. Table 2.1 shows concisely how the various implementations overlap and how they differ. One aim of the Java implementation of MDJEEP was to combine the best parts of both MDJEEP and IBP-NG.

<table><caption>Table 2.1 – An overview of the similarities and differences between the three implementations.</caption><thead><tr><th></th><th>IBP-NG (C)</th><th>MDJEEP (C)</th><th>MDJEEP (Java)</th></tr></thead><tbody><tr><td>Error tolerance</td><td>Tolerance parameter</td><td>Tolerance parameter $\epsilon$</td><td>Tolerance parameter $\epsilon$</td></tr><tr><td>Tree traversal</td><td>Iterative (depth-/breadth-first)</td><td>Recursive (depth-first)</td><td>Recursive (depth-first)</td></tr><tr><td>Coordinate generation (Section 2.3)</td><td>Clifford algebra [129]</td><td>Matrix approach [130]</td><td>Matrix approach [130]</td></tr><tr><td>Intervals (Section 2.4)</td><td>Sampling the arcs</td><td>Refinement by SPG</td><td>Refinement by SPG</td></tr><tr><td>Torsion angles (Section 2.5)</td><td>Used for branching</td><td>Purely distance based (sign is lost)</td><td>Used for branching</td></tr></tbody></table>

To further illustrate the BP algorithm and its framework, the pseudocode for the main recursive step in MDJEEP, based on the Java implementation, is shown in Algorithm 3. In the Java implementation, and in this pseudocode, at every recursive step we build a candidate list $L$ (line 6) for a vertex $v$ (at position $i$ in the order). During the coordinate generation phase (lines 7-13), we populate the candidate list with either directly com-

<footnote n="2">2. https://github.com/geekysuavo/ibp-ng</footnote>

<footnote n="3">3. https://github.com/mucherino/mdjeep</footnote>

<footnote n="4">4. https://github.com/tmalliavin/TAIBP</footnote>
puted candidate positions (three exact distances) or samples from the arcs (one interval
distance). More information about this step can be found in Section 2.3. Next, for each of
these candidate positions, we verify its feasibility (line 15). If the candidate $x_v$ is feasible,
we move on to the next layer $i + 1$. If it is not feasible, MDJEEP carries out an extra
refinement step using local optimization, in order to attempt to satisfy the pruning dis-
tances, before moving on to the next candidate $x_v$. More about this refinement step can
be found in Section 2.4.2.

<caption>**Algorithm 3** The main recursive step of MDJEEP: BP(i, x)</caption>

Input: Index *i* of the vertex *v* in a discretization order *σ*, incumbent solution *x*
1: **if** *i* > |*V*| **then**
2:     **save** current solution *x* ▷ We are done
3:     **else**
4:         *v* ← vertex at position *i* in *σ*
5:         *u*₁, *u*₂, *u*₃ ← reference vertices for *v*
6:         *L* ← empty list of candidate positions
7:         **if** *u*₃ ∈ *E* \ *E*′ **then**
8:             **compute** at most two candidate arcs ▷ See Section 2.3
9:             **add** samples from both arcs to *L* ▷ See Section 2.4
10:         **else**
11:         **compute** at most two candidate positions ▷ See Section 2.3
12:         **add** candidate positions to *L*
13:         **end if**
14:         **for** *x*ₐ ∈ *L* **do**
15:             **if** *x*ₐ is feasible **then**
16:                 BP(*i* + 1, *x*) ▷ Move on to the next candidate vertex in the order
17:                 **else** ▷ Local optimization refinement (see Section 2.4.2)
18:                 **refine** current solution *x* and verify feasibility again
19:                 **if** still infeasible move on to next candidate position *x*ₐ
20:                 **else** BP(*i* + 1, *x*), with refined solution *x*
21:             **end if**
22:         **end for**
23:     **end if**

2.3 Coordinate generation

In the branching phase of BP, when *K* = 3, we compute the intersection of either three
spheres or two spheres and one spherical shell in order to identify the candidate positions
of a vertex $v$. When we have three reference vertices $u_1, u_2$ and $u_3$, this corresponds to solving the following system of quadratic equations:

$$ \left\{ \begin{array}{l} ||x_v - x_{u_3}||^2 = \delta(u_3, v)^2 \\ ||x_v - x_{u_2}||^2 = \delta(u_2, v)^2 \\ ||x_v - x_{u_1}||^2 = \delta(u_1, v)^2 \end{array} \right. \qquad (2.1) $$

In one of the earliest implementations of BP [6], where only exact distances were considered, solutions to the above quadratic system were identified by solving two linear systems [131]. However, this approach lead to a build-up of errors, which was especially noticeable for larger instances. That is why in later works, a different method was used altogether. Instead of solving the quadratic system given in (2.1), we can instead exploit the torsion angle $\tau$, defined by the two planes defined by the triplet of vertices $(u_3, u_2, u_1)$ and $(u_2, u_1, v)$. Figure 2.3 shows the relation between this dihedral angle, and the candidate positions for $v$. The torsion angle $\tau$ can be computed by exploiting the cosine law (see [8]). The necessary distances are available either because of assumption (b), or we may calculate them from the positions $x_{u_1}, x_{u_2}$ and $x_{u_3}$. If $\delta(u_3, v)$ is an interval distance, we can instead compute the corresponding interval $(\underline{\tau}, \bar{\tau})$ and use the methods discussed below to compute the endpoints of the arc. When we compute $\tau$ from the distances we do not have any information about the sign and are forced to consider both candidate positions $v_+$ and $v_-$. Alternatively, $\tau$ can be part of the DGP instance (see Section 1.3). When we have the signed $\tau$ included in the input, we may a priori disregard the branches corresponding to the incorrect sign of the angle. This is further discussed in Section 2.5.

Using $\tau$, the positions $x_{u_1}, x_{u_2}$ and $x_{u_3}$ as well as the three distances $\delta(u_1, v), \delta(u_2, v)$ and $\delta(u_3, v)$, we can compute the Cartesian coordinates $x_v$. Different methods have been used for this step. In [8, 11], a matrix multiplication method was proposed which was very efficient in practice. However, it relied on the consecutivity assumption (mentioned in Section 2.1), which means that $u_1, u_2$ and $u_3$ have to be consecutive in the order. Instead, in [130], a new method was put forward that does not require the consecutivity. Furthermore, practical tests showed that the method performs well and that there is little to no build-up of errors. This approach is used in the BP solver MDJEEP and was also implemented as part of this thesis and used for the experiments discussed in later chapters. This coordinate generation method prevents error build-up by avoiding the accumulation of matrices when we step from one vertex to the next one. Instead, only one matrix $U$ is
<figure><img src="image_11.png" /><figcaption>Figure 2.3 – The candidate positions $v_+$ and $v_-$ can be computed by exploiting the dihedral angle $\tau$. In case we know the sign of $\tau$ (given in the input), we can consider only one of the two candidate positions or arcs.</figcaption></figure>

defined, which can be used to convert a vertex position in the coordinate system defined by $u_1$ to the canonical coordinate system. The three columns of $U$ itself are the unitary vectors $\hat{x}, \hat{y}, \hat{z}$ which correspond to the three axes of the coordinate system centered in $u_1$. The first column, $\hat{x}$ is defined by the vector $\overrightarrow{u_2u_1}$. The third column, $\hat{z}$ can be obtained by the vectorial product $\overrightarrow{u_2u_1} \times \overrightarrow{u_2u_3}$. Finally, the second column $\hat{y}$ is the result of the vectorial product of the other two column vectors: $\hat{y} = \hat{x} \times \hat{z}$. Once we have $U$, we can use it to obtain the coordinates of a candidate position $x_v$ as follows:

$$x_v = x_{u_1} + U \begin{bmatrix} -\delta(u_1, v) \cos \tau \\ \delta(u_1, v) \sin \zeta \cos \tau \\ \delta(u_1, v) \sin \zeta \sin \tau \end{bmatrix} \quad (2.2)$$

where $\zeta$ is the vector angle defined by the triplet $u_2, u_3, v$ which can be obtained using the cosine law. A very similar method for the conversion of torsion space to Cartesian space is presented in [132]. Here, a matrix called $\hat{M}$ (corresponding to the above $U$) is used. The method is referred to as the Natural Extension Reference Frame (NeRF). NeRF is used in the software Rosetta [133], which can be used for analyzing protein structures, which is an application where torsion angles are often encountered.

Aside from the above methods, there are many different ways to compute the candidate positions $x_v$ from the geometric information that we have available during the
branching step. For example, the BP implementation IBP-NG relies on a specific (usually
handcrafted) vertex order and Clifford algebra when computing the vertex positions [127,
129].

## 2.4 Handling interval distances

In the previous section, we discussed the tools that we can use to compute the candidate positions for a vertex $v$. However, in the case that we have a vertex $v$ for which one of the reference distances is an interval, we can only use the coordinate generation methods to obtain the endpoints of the arc as well as the corresponding interval on the dihedral angle $\tau$. In this case, we need to take some extra measures to obtain a discrete number of candidate positions for each vertex $v$. There are two main methods that we can use, which we will discuss below. As mentioned earlier, both these methods are heuristic which means that for interval distances, BP is no longer an exact method, as we cannot guarantee that all possible solutions will be identified.

### 2.4.1 Sampling the arcs

The first method we will look at is straightforward and is how the implementation IBP-NG handles the interval distances (see [127]). The algorithm is parameterized by a parameter, which in IBP-NG is called the *branching parameter* $B$. In case we encounter a vertex $v$ for which we have a reference interval distance, we simply sample $B$ points from the computed arc. In practice, this is done by cutting the interval distance into $B$ pieces. For each sub-interval, we select the middle point as a candidate position for $v$ and use it to create a new branch. For a visual example, see Figure 2.4 (a). This means that in practice, IBP-NG will generally create 2 branches for every $v$ with exact reference distances, and $2B$ for every $v$ with an interval reference distance. In the case that IBP-NG is supplied with torsion angle information, only 1 or $B$ branches are created at every branching step (we can disregard half the branches). The choice of the parameter $B$ is vital for this approach. In case it is picked too small, the chance is high that the correct sample position is missed and the tree is pruned entirely. When it is picked too large, the size of the tree may explode and the algorithm will take very long to terminate. Even when $B$ is a large number, there is a significant chance that this approach will not be able to find a solution for the DDGP instance at hand. This is especially the case in which we
have a large number of interval distances or a high degree of noise. To avoid this problem,
MDJEEP relies on a different method for handling intervals.

2.4.2 A coarse-grained representation and local optimization

In MDJEEP, a combination of a coarse-grained representation of the search space [9]
and a refinement step is used. For this, a resolution parameter ρ is included for the input.
The key point is that we define boxes around each candidate position of a vertex. In a
later refinement step, the points will be allowed to move within these boxes. Consider the
case in which we have a vertex v, and all its reference distances are exact. We are quite
certain about the position of v, so we define a very small box, where its dimensions are
equal to the tolerance parameter ϵ. When we have a vertex v with two exact reference
distances and one interval, we first compute the length of the arc induced by the interval
distance. If the length of the arc is smaller than ρ, we fit the entire arc into one box
(see Figure 2.4 (b)). When it is larger than ρ, we divide the arc into equally sized sub-arcs
smaller than ρ, and cover each such piece with a box (see Figure 2.4 (c)). Within these
boxes, the vertex is initially positioned on the arc or sub-arc, near the center of the box.
For each such box, we create a branch in our tree. When we use this method, the tree
that represents the search space is different from the classical branch-and-prune tree. Each
node no longer only describes a candidate position of a vertex v, but also its surrounding
box.

These boxes come into play in the pruning phase of the algorithm. When we encounter
a pruning distance that is not satisfied we do not immediately prune. First, we run a
constrained local optimization method, minimizing the error on the distances and moving
all the vertices within their boxes. For this refinement step, we use a Spectral Project
Gradient (SPG) method as described in Section 1.4.5. The SPG method is included in the
Java package (see Figure 1.4). In case local optimization manages to satisfy the pruning
distances, we can continue with the next vertex. Otherwise, the branch is pruned. In [13],
it was noted that when we use box sizes that are too small, this can lead to issues for
SPG. This is because bounds that are too strict may harm its ability to converge to a
local optimum. Furthermore, in each refinement step, the reference vertices of a vertex v
may be moved. Therefore, the box corresponding to a candidate position of v should be
expanded. In practice, this issue is largely solved by re-centering the box around the new
candidate position obtained after each refinement step. These expanding and re-centering
routines are also included in the Java version of MDJEEP. Apart from its use for the size
<figure><img src="image_12.png" /><figcaption>Figure 2.4 – In dashed lines, the circles obtained by intersecting the two spheres centered in $u_1$ and $u_2$.</figcaption></figure>

(a) The strategy of selecting a predefined number *B* samples. In this example, *B* = 3, so six branches would be created.

(b) The arc is smaller than the resolution parameter $\rho$. We fully cover the arc with a box and initialize the candidate position $v_+$ in the center of the arc. We create 2 branches.

(c) The arc is larger than the resolution parameter $\rho$. We cover the arc with equally sized boxes and initialize the four candidates positions $v_+$ in the center of their sub-arcs. We create four branches.

of the boxes, the resolution parameter $\rho$ has a second application as well [9]. If at least one solution has already been found we can use $\rho$ to verify if any new solution is too close to this previously found solution. For example, when the BP algorithm is currently exploring alternative candidates for a given vertex $v$ to only consider positions $x_v$ that are at a distance larger than $\rho$ away from the position of $v$ in the previously found solution. Alternatively, after generating a new solution, one can measure the Root mean square deviation (RMSD) between the new solution and the previous one. If the RMSD is lower than some threshold (one could use $\rho$), two solutions should be considered too close, and the current branch should be discarded.

Compared to the sampling method, this refinement approach improves the chances of identifying solutions. This was confirmed in [134], where an intermediate version of Java MDJEEP and IBP-NG were compared experimentally. However, since we are using local optimization, we cannot guarantee that every possible solution will be found.

## 2.5 Torsion angles and the search space

Recall that in Section 1.3, we discussed certain DGP instances where not only distance information is captured, but also torsion angles. These instances are very relevant for the branch-and-prune methods. When we compute candidates for a vertex $v$, given reference
vertices $u_1, u_2$ and $u_3$, the problem of intersecting the three spheres (or two spheres and one spherical shell) can be reduced to the problem of computing and exploiting a torsion angle $\tau$. However, in the case that we have a torsion angle $\tau$ in the instance, we can forego computing it. Furthermore, we will have access to the sign of the dihedral angle, which means we may disregard the branch corresponding to the other sign (see Figure 2.3). Including the torsion angles in the input not only discontinues our use for the distance $\delta(u_3, v)$ but also lets us greatly reduce the size of the search space. Such instances allow for the definition of the DDGP$_3$, specific to dimension 3, based on Definition 3 in Section 1.3. We will use the same weighted directed hypergraph $G(V, E, d)$ to represent the instances. Recall that the hyperedge set $E$ is partitioned into $E_\delta$ (a pair of two singletons) and $E_\tau$ (spanning a quadruplet of vertices), where $E_\delta$ relates to the distances and $E_\tau$ corresponds to a torsion angle in the input. We will further partition the edgeset $E_\delta$ such that the edges $E'_\delta \subset E_\delta$ are the edges that relate to exact distances $\delta(u, v)$. The weighting function $d$ maps the edges $E_\delta$ to the distance $\delta$ and the edges $E_\tau$ to a torsion angle interval $\tau$.

**Definition 6** A weighted directed hypergraph $G = (V, E, d)$ represents a DDGP$_3$ instance in dimension 3 if and only if there exists a vertex ordering on $V$ such that the following two assumptions are satisfied:

(a) $G[\{1, 2, 3\}]$ is a clique whose edges are in $E'$;

(b) $\forall v \in \{4, \dots, |V|\}, \text{ there exist 3 vertices } u_1, u_2, u_3 \in V \text{ s.t.}$

(b.1) $u_1 < v, u_2 < v, \dots, u_3 < v$;

(b.2) $\{ (\{u_1\}, \{v\}), (\{u_2\}, \{v\}) \} \subset E'_\delta$ and either $(\{u_3\}, \{v\}) \in E_\delta$ or $(\{u_3, u_2, u_1\}, \{v\}) \in E_\tau$;

(b.3) $\mathcal{V}_S(u_1, u_2, u_3) > 0$

The main difference can be found in assumption (b.2), where we now require either the edge $(\{u_3\}, \{v\}) \in E_\delta$ or the torsion angle formed by the reference vertices $\{u_3, u_2, u_1\}$ and $v$ to be included in the instance $((\{u_3, u_2, u_1\}, \{v\}) \in E_\tau)$. Consider the case where, given the vertex order, for every vertex $v$, we have such a torsion angle $\tau(u_3, u_2, u_1, v)$ available. If all reference distances are exact and the torsion angles are precise values, we will end up with a very narrow tree, with only one solution (given a feasible instance). In case we have an instance including reference interval distances, the search space will only grow as a result of these interval distances, and the lengths of the induced arcs (or the branching parameter $B$, if we branch by sampling). This case is of interest because when dealing
with certain instances related to protein determination, we have access to the torsion angle $\tau(u_3, u_2, u_1, v)$ for (almost) every vertex $v$ the vertex order (see Section 3.3). For these instances, we can greatly reduce the search space when we directly use the torsion angle values in the branching phase. This is done by IBP-NG, but not by the original C-version of MDJEEP, which is purely distance-based. The Java implementation of MDJEEP does use the torsion angles during branching. Of course, the solutions identified by the Java-version of MDJEEP will be a subset of the solutions found by the C-version. However, in practice we often stop execution after a certain number of solutions are identified. These first solutions have a higher chance of being accurate when we use the sign of the torsion angles in the branching phase.
# MOLECULAR STRUCTURE
DETERMINATION

One of the most important applications of the DGP, and the main topic of this thesis, is *molecular structure determination* [27, 33, 37–45]. The work done in the context of this thesis focuses only on *proteins* which are a specific type of molecules. However, a lot of the work may be extended to other types of molecules as well. Proteins are molecules that perform vital functions in our bodies. They make up about 15-20% of the human body [135]. Additionally, proteins are important components of medicines, such as vaccines. How the human body reacts to a medicine relates strongly to the three-dimensional structure of the proteins that make up the treatment. This makes protein structure determination a topic of great practical importance and by extension one of the most researched applications of the Distance Geometry Problem. In this application, the vertices $V$ of the graph $G$ represent atoms of a given molecule, and the distance information $d$ associated with the edges $E$ of the graph reflects inter-atomic distances. These distances may either be derived from the chemical geometry of the molecules or experimental techniques, such as Nuclear Magnetic Resonance (NMR) [2]. The DGP instances that we have in the context of this application satisfy the special conditions discussed in Chapter 2 that let us discretize the search space [6, 8, 26], and use an algorithm that can enumerate all the possible conformations of the protein at hand.

We will start this chapter with an introduction to the problem of protein structure
determination (Section 3.1). Next, in Section 3.2, we review the experimental and machine
learning methods that we can use. In Section 3.4 we look at which distances are available
in the protein instances and we review the use of different DGP and DDGP methods in
the literature. After this, Section 3.5 will cover different experimental works done as part
of this thesis, including experiments with real protein data. Finally, in Section 3.6 we will
summarize the chapter and draw some conclusions.
<table>
  <tbody>
    <tr>
      <td>3.1</td>
      <td>The problem of Protein Structure Determination</td>
      <td>70</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Experimental techniques and machine learning</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>X-ray crystallography</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.2.2</td>
      <td>Cryo-Electron Microscopy</td>
      <td>73</td>
    </tr>
    <tr>
      <td>3.2.3</td>
      <td>Nuclear Magnetic Resonance Spectroscopy</td>
      <td>74</td>
    </tr>
    <tr>
      <td>3.2.4</td>
      <td>Machine learning and ALPHAFOLD</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Inter-atomic distances and torsion angles</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Distance-based approaches</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>Simulated Annealing</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.4.2</td>
      <td>General DGP methods</td>
      <td>82</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>Exploiting the discretizability of the instances</td>
      <td>84</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Experiments and results</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.5.1</td>
      <td>Branch-and-prune with real NMR instances</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.5.2</td>
      <td>Covalent geometry and branch-and-prune</td>
      <td>90</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Discussion</td>
      <td>97</td>
    </tr>
  </tbody>
</table>
**3.1 The problem of Protein Structure Determination**

Proteins are molecules that are part of many vital functions of the human body, from immunity to brain activity. They consist of a chain of smaller molecules, called *amino acids* or *residues*. There are 20 different amino acids that can exist within a protein. Each amino acid has the same core structure, which consists of a central carbon atom (CA) bonded to an amino group (NH2) a carboxyl group (COOH), and a hydrogen atom (HA). This core structure is known as the *backbone* of an amino acid. Aside from the backbone, the residues have a group of atoms that make each type of amino acid unique. This group of atoms is attached to the CA atom and is known as the *side chain*. Amino acids chained within a protein are linked together by *peptide bonds* which are formed when the carboxyl group of one amino acid combines with the amino group of a second. The backbones of the amino acids in combination with the peptide bonds form the backbone of the protein.

Proteins are described by different levels of structure. The most simple level is the primary structure, which refers to the linear sequence of amino acids that make up the protein. The secondary structure of a protein refers to the structural patterns that combinations of amino acids often conform to. For example, groups of amino acids often form secondary structure patterns such as alpha helices and beta sheets. Finally, the tertiary structure of a protein is the three-dimensional conformation of the atoms that make up the molecule. Essentially, the tertiary structure of a protein describes how its different amino acids *fold* into one, complex three-dimensional structure. Figure 3.1 visualises the different structural levels of a protein.

<figure><img src="image_13.png" /><figcaption>Figure 3.1 – Three different levels of protein organization. Left: primary structure (linear sequence of amino acids). Middle: secondary structure (structural patterns, alpha helix, and beta sheets). Right: tertiary structure (three-dimensional confirmation of proteins). Image modified of original provided by National Human Genome Research Institute<fnref n="1" />.</figcaption></figure>

<footnote n="1">1. https://www.genome.gov/genetics-glossary/Protein.</footnote>
If we have an unidentified protein, different techniques may be used to determine its primary structure, such as Edman degradation [136–138] and mass spectrometry [139–141]. Once we know its primary structure, we have identified the protein. However, to know its function we need to predict, compute, or measure its tertiary structure. This is known as protein structure determination. A protein in vivo will tend to fold in such a way that its total energy is minimized. The potential energy of tertiary protein structure x can for example be estimated using a force field f, which are computational models that estimate forces between atoms of a protein:

$$f : \mathbb{R}^{Kn} \rightarrow \mathbb{R}$$

An example force-field function is AMBER [142], which relies on the chemical composition of the protein structure x and force field parameters (such as $r_{eq}$, $\theta_{eq}$, $K_r$ and $K_{\theta}$) to compute the global energy $E_{total}$:

$$E_{\text{total}} = \sum_{\text{bonds}} K_r (r - r_{\text{eq}})^2 + \sum_{\text{angles}} K_{\theta} (\theta - \theta_{\text{eq}})^2 + \\ \sum_{\text{dihedrals}} \frac{V_n}{2} [1 + \cos(n\phi - \gamma)] + \sum_{i<j} \frac{A_{ij}}{R_{ij}^{12}} - \frac{B_{ij}}{R_{ij}^{6}} + \frac{q_i q_j}{\epsilon R_{ij}}$$

More information about this function is found [142], while these types of angles and other force-field parameters used in this function are discussed in Section 3.3. Therefore, protein structure determination (generally in the context of computer simulation) is regarded as a minimization problem [143, 144]:

**Definition 7** Given the primary structure P of a protein, determine a tertiary structure x that respects the chemical composition implied by P such that a force field function $f(x)$ is minimized.

There are different approaches that one can take to try to compute or predict protein structures, these include experimental techniques, machine learning methods, and several distance-based approaches. Some of these methods, in particular the methods that we will use in the experiments in this chapter do not rely on energy terms but will instead define the protein folding problem as a version of the DGP or DDGP. In the next sections we will explore the various methods that one can use to determine protein structures.
## 3.2 Experimental techniques and machine learning

There are various (non-distance-based) approaches to protein structure determination. In this section, we will look at two popular and effective groups: techniques relying on physically analysing the protein and methods that exploit large datasets of known structures using machine learning, mostly by employing neural networks. These approaches are connected because the known structures that make up the datasets that are used for training the neural networks have often been measured using one of the experimental techniques.

### 3.2.1 X-ray crystallography

Most known protein structures have been measured using X-ray crystallography [145–148]. The first step in this process is to purify and the protein and then crystallize it. Obtaining this crystal is not an easy step, as it requires conditions that promote the formation of a regular, repeating crystal lattice. Different techniques are used to grow crystals, such as vapor diffusion [149] and (micro) batch crystallization [150]. Once the crystal has been grown the next step is to mount it on a goniometer inside an X-ray diffractometer. Then, the crystallized protein is exposed to a series of X-ray beams. The crystal causes the incident rays to diffract in many directions onto a screen. These directions, angles, and the intensities of the diffracted beams are then analyzed which allows for the determination of the distribution of the electrons in the protein. Techniques such as Fourier analysis and Patterson synthesis are employed in this step. The resulting map of the electron density can in turn be used to deduce the position of each atom in the crystallized protein, giving us the three-dimensional structure. Once the three-dimensional structure is obtained, it undergoes validation to ensure accuracy and reliability, generally this involves checking its geometric feasibility.

The structures measured with X-ray crystallography are generally of very high quality, but can lack precision, because using the X-ray beams causes a small portion of radiation damage to the proteins. Aside from this, there are several other limitations. Firstly, this method works well only on proteins that have rigid structures that form nice, ordered crystals. It fails to measure the structure of flexible proteins because X-ray crystallography relies on having many molecules aligned in the same orientation. The portions of a protein that do not follow this pattern will be invisible in the electron maps which in turn leads to missing coordinates [151]. Another shortcoming is that the structure that the protein
forms inside the crystal is not necessarily the natural confirmation that it has in other environments.

3.2.2 Cryo-Electron Microscopy

A more recent technique avoids the need for crystallization and X-ray beams. This eliminates the issue of the X-ray beams causing damage to the protein sample. CryoEM [152–154], short for Cryo-Electron Microscopy, relies on cooling the molecules down to cryogenic temperatures, preserving them in a thin layer of non-crystalline ice. To do this, the proteins are applied to a grid covered with a thin layer of support material. An example of a support material that is often used is carbon [155]. The excess solution is blotted away, leaving a thin film of sample across the holes of the grid. The grid is then rapidly plunged into a cryogen (such as liquid ethane or propane) to freeze the sample [156]. Once the sample is frozen, the grid is transferred to a cryo-electron microscope. The microscope operates under vacuum conditions to minimize electron scattering. A beam of electrons is transmitted through the frozen protein, and the interactions of electrons with the sample produce an image. Several such images are taken from different orientations of the sample, resulting in a collection of 2D projection images. The collected 2D images are processed using computational methods to reconstruct a 3D model of the specimen. This process, known as single-particle reconstruction [157], involves aligning and averaging the 2D images to improve signal-to-noise ratio and resolve the 3D structure. Advanced image processing techniques, such as classification and refinement algorithms [158], are employed to sort and refine the images to obtain a high-resolution 3D reconstruction.

The CryoEM technique is more recent and is able to rival [159] and in some cases even surpass [160] the quality of the structures found by X-ray crystallography. CryoEM is still being extensively researched and improved. Despite these advantages, there are several drawbacks [161]. It suffers the same limitation as X-ray crystallography when it comes to flexible proteins: we can only measure to protein structure found in the molecule's frozen state. Furthermore, X-ray crystallography is a much faster technique than CryoEM and thus has a higher throughput. Work is also being done to combine the results of CryoEM and X-ray crystallography in one model. This process is known as *integrative modeling* [162]. Both methods give us access to high-resolution protein structures, albeit structures of the protein in a specific, fixed state. We will see that when we use specific distance-based approaches we can overcome this issue. Applying these experimental techniques is an expensive process requiring large machines and a
lot of power. Therefore, it is desirable to have different methods, which rely only on
computation. One promising avenue of approach is based on machine learning.

3.2.3 Nuclear Magnetic Resonance Spectroscopy

The above two methods are effective for measuring the structure of a protein. However, one main shortcoming is that before the measurements are made, the protein is fixed in space in a way where, either by freezing or by capturing it in a crystal. This means that the structures that we measure using X-ray crystallography or CryoEM do not necessarily correspond to the way a protein will behave in vivo. Moreover, these methods are not well suited for flexible proteins. When we use solution Nuclear Magnetic Resonance (NMR) spectroscopy [2, 163], a protein is dissolved in a solvent (often water), so that it is free to rotate and tumble naturally. NMR spectroscopy relies on the fundamental quantum mechanical property of nuclear spins. The protein (in solution) is placed in a strong and homogeneous magnetic field. In this magnetic field, atoms with nonzero spin numbers are at different energy levels. When a radio-frequency pulse is applied to the sample at a specific frequency, it matches the energy difference between the spin states. Some of the atoms absorb energy and transition between energy states. The amount by which the energy transitions is called the *chemical shift* of the atom. After this process, the resonant frequency of the nuclei of the atoms can be measured, which is known as *resonance*. Using a Fourier transform we can convert these resonances and chemical shifts to a spectrum. Such spectra allow for the assignment of each peak in the spectrum to a specific atom within the protein. In Section 3.2.3 we see an example of such a spectrum, created using the NOE [164]. Such a spectrum is referred to as a NOESY spectrum, and it allows it to identify through-space interactions within the molecule. The Nuclear Overhauser effect is what happens when two atoms of the same type with similar nuclei (e.g. ¹H, ¹³C and ¹⁵N) are close together and engage in *cross-relaxation*, which means that energy is being transferred between their ions. This is visible in the spectrum because, at the intersection of the chemical shifts of the atom pairs, we will see *cross-peaks*, indicating that the offending atoms must be close together in the molecule [165]. This type of spectrum is of great interest to us because it allows us to identify an interval distance between light-weight atoms (such as hydrogen atoms). Aside from this, the chemical shifts in the two-dimensional ¹³C and ¹⁵N spectra give us information about the φ and ψ dihedral angles in the backbone of the proteins. More details about these distance restraints and torsion angles can be found in Section 3.3. The restraints obtained from the NMR experimental
data can then be exploited to generate protein structures. The standard method of doing
this is by using the meta-heuristic Simulated Annealing [166, 167], which is discussed
in Section 3.4.1. However, as we will see in Sections 3.4.2 and 3.4.3, these distances can
be exploited by DG methods as well.

<figure><img src="image_14.png" /><figcaption>Figure 3.2 – An example of a two-dimensional NOESY spectrum for a small molecule, where the chemical shifts on the axes correspond to $^1$H protons (labeled in ppm). At the top of the spectrum, the different frequencies are labeled with the hydrogen that emits it. In the bottom right there is a two-dimensional diagram of the molecule that is being analysed. The red peaks in the diagonal are the self-peaks for each of the hydrogen atoms. The cross-peaks (blue) away from the diagonal are the ones we are interested in. For example, the indicated cross-peaks in the figure tell us that both $H_8$ atoms are in close proximity to the $H_4$ and $H_7$ atoms in the molecule. We see that there is a peak for $H_{8b}$ and the $H_3,H_9$ pair, but no peak between the pair and $H_{8b}$, possibly indicating that $H_{8b}$ is not close to $H_3$ and $H_9$. Figure taken from [168] and modified.</figcaption></figure>
### 3.2.4 Machine learning and AlphaFold

In some cases, when experimental techniques are expensive or otherwise unpractical, other approaches must be used. Some of the most successful methods for protein structure prediction are knowledge-based. Knowledge-based methods involve learning and extracting knowledge from existing solved protein structures (using the methods from the previous section) and generalizing this knowledge to unknown protein structures. Note that we use the word prediction instead of determination because we have no guarantee of the correctness of the output. Machine learning methods automatize this process of knowledge extraction. At the start, a dataset of known structures is partitioned into two parts: a training set and a test set. The training set is used to train the machine learning model so that it becomes capable of recognizing patterns between the primary structure of the proteins and their tertiary structure. Once it has been trained sufficiently, it can now be used to predict the tertiary structure of the proteins in the test set. Afterward, the accuracy of the model is measured by comparing the computed structures to the known structures in the test set. Once we have identified a model with high enough accuracy, we can start to use it on proteins for which the structures are unknown.

The type of machine learning that has gained the most popularity and has attained the best results is *deep learning* [107] using *neural networks*. Recently, ALPHAFOLD [169] by DeepMind has produced structures with very high accuracy. The process of using deep learning for protein structure prediction generally consists of three steps [170]: (1) an input module, which converts the primary structure of a protein into different features, (2) a neural network that transforms these features to spatial information, and finally (3) an output model that uses to spatial information to compute and refine the three-dimensional conformation of the protein. The resulting structure may be further refined in an additional step if necessary. For each of these steps, different approaches have been considered.

1. The features extracted from the primary structures almost always consist of Multiple Sequence Alignments (MSAs) between homologous proteins [171]. Such alignments can be used directly as raw features (ALPHAFOLD) or they can be converted into different features [172–175].

2. The choice of neural network in the second step is of vital importance. The first neural networks used were shallow convolutional neural networks (CNNs). These are a type of neural networks that work very well in the context of computer vision [107].
The motivation for using CNNs was that the features obtained from the MSAs
resemble an image. Modest accuracies were achieved using CNNs [172, 174]. This
is because CNNs have several shortcomings when it comes to protein structure
prediction. For example, CNNs assume translational invariance as to where patterns
appear in the input. This is true for general images, but in the case of proteins, the
position in which a pattern occurs in an amino acid sequence is very relevant to its
final structure. Next, CNNs were enriched with deep residual networks [173, 175],
booking some great improvements. Finally, ALPHAFOLD introduced the use of the
attention networks, which was a major breakthrough [169]. Such attention networks
are powerful architectures that are capable of selectively focusing on certain parts
of the input data that are relevant to the task at hand. Instead of assuming that
local patterns are most important, these networks initially place equal weights on all
possible interactions, and then later refine these weights during the training phase.

3. As output, the first deep learning models produced binary contact maps [172, 173], based on the distances between amino acids. These binary contact maps must then be converted into folded proteins by using protein modeling software such as ROSETTA [133]. Other models improved on this by using discretized distances (using small increments) instead of binary contacts [175, 176]. This improved the quality of the predicted structures. However, an extra step was still necessary to transform these amino acid distance maps into folded protein structures. The final version of ALPHAFOLD improves on these approaches by being able to directly output the three-dimensional structures, based on a method proposed in [177].

**3.3 Inter-atomic distances and torsion angles**

Before we discuss the distance-based methods, we will explore the different types of
inter-atomic distance restraints and dihedral angles that we can deduce from the primary
structure of a protein. As mentioned before, not only the primary structure but also the
spectra obtained from NMR let us derive information about distances and torsion angles.
Some of these distance restraints are precise values, while others give us only a lower and
upper bound on the distance. As mentioned in Section 1.3, these torsion angles may be
converted to distances and included in the DGP instance, so that general DGP methods
may use them. However, as we saw in Section 2.5 if we are using branch-and-prune, we
may directly exploit the signs of the torsion angles during the branching phase. As part
of this thesis, Java code was written that derives all these distance and torsion angle
types directly from STAR-NMR files (see Section 3.5.1). In the Java package the code for
parsing STAR-NMR files can be found in the NMR class (see Figure 1.4). The code that
represents the protein structure determination problem as a DGP is found in the MDGP
class. We will now explore the different distance types and torsion angles.

Type 1: van der Waals radii

The van der Waals radius of an atom $u \in V$ is the radius of an imaginary sphere which represents the closest any other atom $v$ may approach $u$ [178]. This means that the sum of these radii can be used to provide an expected lower bound to the distance between any pair of atoms. The van der Waals radius is specific to the type of atom. This radius allows us to describe the electron cloud around the nucleus of an atom as a sphere. However, due to polarisability, these clouds are never *real* spheres, such that a pair of atoms can actually be closer than the sum of their van der Waals radii. In other words, the spheres in the van der Waals model can be considered as *soft* spheres, meaning that they can slightly penetrate one another, which is why generally only a portion of this sum is used as the lower bound on the distances. On the one hand, van der Waals distances are not very informative because they only carry a lower bound, while their upper bound is generally very large (initially set to infinity, it can be improved by using triangular inequalities involving related distances). On the other hand, these distances are very abundant because they can be defined for every pair of atoms that are not bonded. Furthermore, if these restraints are respected, they make sure that clashes are avoided between atoms within the computed molecular structures. In the experiments conducted as part of this thesis, not all van der Waals distances are included in the instances. This is because only a subset is necessary to avoid clashes, and a very large number of distances will greatly slow down the refinement step in MDJEEP.

Type 2: covalent geometry, force fields, and torsion angles

These distances and torsion angles are based on the *covalent geometry* of the molecules,
which refers to the way the different atoms are bonded together. Because we know the
primary structure of the protein, we know precisely which atoms are bonded together in
which amino acid. This means that we can use information from *force fields* to derive
distances and torsion angles. Alternatively, this information can also be extracted from
protein structures in the Protein Data Bank (PDB), if the structure of the protein at hand
was already measured or computed by a different technique (when generating artificial
instances, see Section 3.5.2). As mentioned in Section 3.1, force field models define a
function for the potential energy of a molecule, where the parameters of this function
encapsulate information about the forces between atoms. One such parameter that is of
interest to us is the bond lengths. For every pair of bonded atoms (u, v) we may derive an
exact distance between u and v, which differs based on the types of atoms involved as well
as where they are located within the amino acid. A second parameter of interest is the bond
angles. For every pair of atoms (u, w) that are bonded to a common atom v, this parameter
describes the angle that the three atoms form. This, in combination with the two bond
lengths that we know, allows us to compute the diagonal distance between u and w. This
bond angle distance is regarded as exact, even though there may be small variations in
the covalent geometry of different proteins. Section 3.5.2 gives more detail on this issue of
variance. Examples of force fields that one can use are AMBER [142], CHARMM [179],
and PARALLHDG [180]. Different force fields may have slightly different parameters. In
the Java package described in Section 1.5, the FORCEFIELD class is responsible for deriving
these distances, and supports the AMBER and CHARMM forcefields. Apart from bond
lengths and bond angles, the primary structure of the protein also lets us deduce values
on specific torsion angles. Firstly, we have information about a torsion angle ω in the
backbone of the protein. This is a backbone torsion angle that is formed by a sequence
of four atoms between two consecutive amino acids. For ω, the sequence of four atoms
consists of the CA and C of the first amino acid and the N and CA of the next amino
acid, describing the rotation around the peptide bond between C and N. The angle ω is
almost always close to 180°, except when the peptide bond is cis, in which case it is close
to 0°. In experiments from the literature, the torsion angles ω are often fixed directly
at 180°, even though in practice there can be quite large variations (see Section 3.5.2).
Another type of torsion angles, known as *improper* torsion angles, can be taken from force
field parameters. These torsion angles regard sequences of atoms (v₁, v₂, v₃, v₄), where we
already know the distance δ(v₁, v₄) from a bond length or bond angle. This means that
for general DGP methods, they do not give any extra information. However, exploiting
the sign of these torsion angles is very valuable for DDGP methods, as they do give a lot
of information about the structure of the proteins at hand (see [181]). For example, the
improper torsion angle created by the atom sequence N, CA, C, and HA (all in the same
amino acid) tells us whether the amino acid is left- or right-handed. They are included in
the instances used in the experiments in this thesis and exploited by Java MDJEEP.
### Type 3: NMR restraints from NOESY spectra

As described before, the NMR spectra obtained from NMR experiments can be used to derive distance restraints. The spectrum that is relevant for this distance type is the ¹H two-dimensional NOESY spectrum, for which an example is shown in Section 3.2.3. This is a spectrum in which each peak corresponds to a chemical shift of a hydrogen atom (or to noise). Once a cross-peak is identified in such a spectrum, we know that the distance between the two hydrogen atoms that correspond to the peaks should be between 1.8 and 5Å [182]. This lets us enrich the distance graph *G* with several such interval distances between hydrogen atoms that are relatively close together.

### Type 4: Distances from $\phi$ and $\psi$ backbone torsion angles

Aside from the distances between hydrogen atoms discussed in the previous paragraph, the NMR spectra can offer us more. The chemical shifts in the two-dimensional ¹³C and ¹⁵N spectra give us information about the $\phi$ and $\psi$ dihedral angles in the backbone of the proteins. In the protein backbone, just like the $\omega$ angle mentioned before, these angles are formed between two consecutive amino acids. The $\phi$ angle corresponds to the dihedral angle formed by the atoms N, CA, and C of the first amino acid and the N of the next amino acid, describing the rotation around the CA - C bond. The $\psi$ angle corresponds to the dihedral angle formed by the atom C of the previous amino acid and the atoms N, CA, and C of the next amino acid, describing the rotation around the N - CA bond. The chemical shifts from the NMR spectra, in combination with machine learning tools such as TALOS+ [183], let us predict the secondary structure of the proteins with high accuracy. Knowing the secondary structure allows for the derivation of bounds on the backbone dihedral angles $\phi$ and $\psi$. The values for the torsion angles can then be exploited in the same way as the $\omega$ dihedral angle.

## 3.4 Distance-based approaches

In this section, we will discuss how we can exploit inter-atomic distances (and torsion angles) to compute protein structures. First we will look at the standard method that is used for generating structures from NMR experimental data: Simulated Annealing. SA works by solving the optimization problem given in Definition 7, while observing the distance and angular restraints given by the NMR data. Next, we will look at methods
that regard the protein structure determination problem as an instance of the Distance Geometry Problem where the dimension $K = 3$. In this case, the vertices $V$ of the graph $G = (V, E)$ are the atoms of the protein and the edges $(u, v) \in E$ together with their weight $\delta(u, v)$ describe inter-atomic distances. This variant of the DGP is often referred to as the Molecular Distance Geometry Problem (MDGP) in the literature [3, 8, 10–12, 22, 33, 94, 95, 97, 102]. Note that the distance graphs for these DGP instances, especially those with interval distances, will likely have many possible embeddings. Many DGP tools will find only one such embedding or several possible embeddings. That is why solving the DGP instances in the context of protein structure determination is often referred to as *conformational sampling* because the methods only sample some of the many possible conformations. Different such methods from the literature are reviewed in Section 3.4.2. Another class of algorithms, present within the branch-and-prune framework, is capable of enumerating all possible conformations. This method is discussed in Section 3.4.3, and is the approach that we used in this thesis for generating protein structures.

### 3.4.1 Simulated Annealing

As mentioned in Section 3.1, a protein in vivo will tend to fold in such a way that its total energy is minimized. Therefore, the goal becomes to compute a protein structure that satisfies the identified distance constraints while minimizing the global energy of the molecule. This optimization problem is not straightforward because one must avoid getting stuck in a local minimum. Simulated Annealing (SA) [166, 167] is used to overcome exactly that. SA is a constrained meta-heuristic search that attempts to approximate a global optimum. It was inspired by annealing in metallurgy, where a metal is heated to a high temperature quickly and then slowly cooled, which increases its strength while making it easier to work with. SA attempts to mimic this as a global optimization approach with a high-temperature phase and a cooling phase.

The SA process starts by generating an initial solution, which will likely have a high global energy. In the context of protein structure determination, different ways of producing initial solutions have been proposed, generally based on random generation [184]. In this initial phase we also initialize a temperature control parameter to a high value. Then we will attempt to minimize the energy of our solutions while gradually lowering the temperature as we attempt to reach a state of minimum energy. At each iteration, a new solution is generated and its energy is measured. In the classic SA approach, this new solution is generated randomly. However, in the context of protein structure deter-
mination, other strategies are employed more often used. For example, in [184], a genetic
algorithms based method is proposed to choose a neighboring solution. Once we have the
next solution, we compute its global energy. If the new energy is lower than the current
energy we accept this new solution as better. However, to avoid getting stuck in a local
optimum, sometimes a solution with a higher energy is accepted as well. This is done us-
ing a specific probability which depends on the temperature control parameter. When the
temperature is high, we accept worse solutions (in terms of energy) with a higher proba-
bility. This way, as the temperature control parameter declines, we converge towards an
(approximate) global minimum.

Several computer programs have been developed based on the above algorithm. Exam-
ples of software tools that are used widely in practice are XPLOR-NIH [185], CYANA [186]
and GENMR [187]. XPLOR-NIH is a powerful tool, which can directly take NOESY
spectra as input, where it can automatically detect crosspeaks, avoiding peaks arising from
noise. It can be used to read input data from X-ray crystallography as well. The software
package CYANA also works directly on NOESY spectra [188]. GENMR was developed
to offer a tool that works as a web service. It uses a combination of XPLOR-NIH and
ROSETTA [133] to generate structures at the back end.

All of these tools provide very good results when tested on real data [189]. These
software tools allow for the integration of X-ray crystallography data as well as in some
cases CryoEm measurements. Combining these experimental techniques with NMR mea-
surements leads to integrative models that can predict protein structures with very high
atomic resolution. Currently, these energy minimization tools are the best approach for
structure determination by NMR if one is interested in sampling the conformational space
of a protein. However, there is a drawback to these methods. While they generate a set of
possible structures with low potential energy, there is no guarantee that these are all the
possible low-energy states the protein could take (because of the use of a meta-heuristic).
Furthermore, proteins can sometimes move from one state to another and shortly be
at a point where their potential energy is not minimal. To model these dynamics other
approaches may need to be used.

3.4.2 General DGP methods

Several general DGP, discussed in Section 1.4, have been tested in the context of
protein structure determination and conformational sampling. However, none of these
methods have been tested on distance constraints derived from real NMR experiments.
Instead, these algorithms were tested on artificially generated instances. Some of these
instances were randomly generated, and tend to follow a general protein-like structure.
For example, several of the meta-heuristics VNS and Multi-start as well as the global
optimization method SBB were tested on small random instances, providing good re-
sults [94]. However, these instances only had exact distances and do not closely resemble
DGP instances derived from NMR experiments at all.

In other experiments, known protein conformations are taken from the Protein Data
Bank (PDB) [190], which are structures that are generally measured using one of the
experimental techniques that were discussed before. Next, distances are generated fol-
lowing a set of rules, which generally aim to mimic instances from real NMR data. The
reason for using these artificial instances is because real NMR data contains a large de-
gree of uncertainty and noise, which cannot be handled by these general DGP methods.
The graph decomposition algorithm ABBIE was tested on such instances, but only ex-
act distances were considered [90]. The BUILD-UP algorithm was tested on similar exact
instances with great results [119]. The simple heuristic Stochastic Proximity Embedding
(SPE) was shown to work well on artificial instances where both intervals and exact
distances were included, providing results that are competitive with other popular confor-
mational sampling methods [191]. GNOMAD was also tested on such instances. In these
instances, only distances between “close” atoms were included, which are all distances
lower than 10Å [98]. Yet another step closer are test instances with interval and exact
distances only between atoms in consecutive amino acids. Both the smoothing-based al-
gorithms, DGSOL and DCA have been tested on such instances with good results [44,
95, 96]. In [103], Monotonic Basin Hopping (MBH) was tested on similar instances and
was shown to outperform DGSOL in most cases. The EMBED algorithm was tested on
instances with exact and interval data as well, where all pairwise distances smaller than
4Å were included. The resulting structures were of good quality. The Alternating Projec-
tions Algorithm was tested on artificial instances that are very close to real data. Here,
each of the distance types was aptly simulated. The algorithm was tested on a specific
protein (bovine pancreatic trypsin inhibitor) and was shown to work well [84].

General DGP methods can work well for artificially generated instances, including
instances with large proteins. However, some of these algorithms work only on instances
with exact distances and are thus not useful in the context of instances derived from NMR
spectra. Others can work on instances with intervals, but the bounds that were used in
the experiments were smaller than they should be and many more distances were included
compared to what is present in real data. It is hard to compare these algorithms fairly, as
they may work well for specific instances and not so well for others. In the next section,
we will look at methods that exploit the discretizability of the instances related to protein
structure determination.

3.4.3 Exploiting the discretizability of the instances

The distance types described in Section 3.3 allow for the discretization of the DGP
instances related to protein structure determination. In particular, for this, we require
the information about bond lengths, bond angles, the ω angle (**type 2**) and the φ, and
ψ torsion angles (**type 4**). These distances and dihedral angles are enough to satisfy
the assumptions given in Definition 6. The other geometric information may be used
as pruning devices. Different branch-and-prune methods have been tested on instances
related to protein structure determination [7, 9–14]. The advantage of this approach is that
BP theoretically will allow us to identify all protein structures that satisfy the geometric
information captured in the DDGP instances. For all the papers cited above, just like for
the general DGP methods, the DDGP instances were artificially generated. As mentioned
in the previous section, one cause of this uncertainty is noise in the NMR data, which
will be further investigated in Section 3.5, where we present experiments with genuine
NMR instances. Secondly, the BP methods regard distances from the covalent geometry
as exact, while in practice they may vary a lot. More information about this is given
in Section 3.5.2, where we experimentally study the effect of the variations in the covalent
geometry.

Initial implementations of branch-and-prune, in 2009, had trouble with small amounts
of noise and required the input distances to be of very high precision [192]. This was
improved on in the same year [10] by including a tolerance parameter for the pruning
step. Test instances were generated by including all distances smaller than 6Å. Different
amounts of noise were then introduced into a subset of these distances, creating interval
distances. The results showed clear correlations between the amount of noise (the range
of the intervals), the number of found solutions, and the quality of the found solutions.
Next, in 2013, experiments were done with instances generated with all distances shorter
than 5Å and a small amount of noise (0.3Å) [11].

In 2010, MDJEEP was first released [7] and was tested on generated instances with
only exact distances. Later releases of MDJEEP [9, 13, 14] could also handle interval data.
In particular, in [9], experiments were done with artificial instances that are getting closer
to genuine NMR instances. The instances were generated with special rules for some of
the distance types discussed in Section 3.3. For **type 2**, the bond length and bond angle
distances, exact values were used with three decimals. For **type 3**, interval distances with
a range of 0.5Å are used. It is important to note here that for these distances from NMR
restraints, only distances between hydrogen atoms closer than 5Å were kept. This is a
clear improvement as compared to previous experiments, where all distances smaller than
5Å are included. Lastly, for **type 4**, intervals with range 0.1Å were created. The methods
were tested on some small proteins and peptides, showing good results. However, note
that while these experiments got closer to genuine data, the ranges on the intervals are
still quite small compared to what they are in real NMR data. In the next section, we will
present results where MDJEEP was tested on instances created from genuine NMR data.

**3.5 Experiments and results**

In this section, we will discuss the several experiments conducted throughout this
thesis, where we used branch-and-prune methods to compute protein structures.

**3.5.1 Branch-and-prune with real NMR instances**

As mentioned in the previous section, the BP implementations are often only tested
on artificially generated data, simulated in a way to resemble real NMR distance data (see
for example [14]). This motivated a study on the different distance types, published and
presented at the Computational Structural Bioinformatics Workshop, which is a satel-
lite workshop of the IEEE Bioinformatics and Biomedicine conference BIBM [27]. In this
study, different experiments were conducted using real data as well as a mix of real and
artificial distance data. The goal was to see which type of distances have the largest ef-
fect on the quality of the embeddings as well as the feasibility of the BP algorithm. The
structures were computed using the C-version of MDJEEP. In this section, we will focus
on the experiments with the instances from real NMR data. Furthermore, the structures
that we will discuss here are computed with the Java implementation of MDJEEP, which
exploits the signs of the dihedral angles (see Table 2.1 for the differences). In all experi-
ments with proteins presented in this thesis, we will only consider atoms in the backbone
of the proteins. The atom names of the atoms considered for each amino acid are: N, CA,
C, O, OXT, H1, H2, H and HA. We only consider the backbone atoms, because once
we have the backbone, placing the side chains becomes a straightforward post-processing step. This final step can be done using protein analysis software such as ROSETTA [133].

In order to be able to handle real NMR data, the Java package created for this thesis is capable of reading a STAR-NMR file. STAR-NMR is an extension of the more general STAR format, which is a text-based file format for storing structured data. STAR-NMR files describe distance constraints and torsion angle constraints that result from the spectra obtained from NMR experiments. Furthermore, parsers were added for two different force fields: AMBER [142] and CHARMM [179]. The NMR data, combined with force field parameters and the primary structure, lets the package generate $DDGP_3$ instances (Definition 6) from real data. Even when there is no NMR data available, or the distance information is sparse, we can still generate discretizable instances. The NMR restraints (type 3) are only used for pruning so are not required for discretization. The NMR torsion (type 4) angles are used for discretization, but when they are missing they can be replaced by full 360° intervals (more about this later in this section). This instance generation tool was extracted from the Java package and can be found at a public GitHub repository<fnref n="2" />. This tool does not only generate instances from NMR data, but also from X-ray crystallography data (see Section 3.5.2). Once we have the instances generated, we can proceed with branch-and-prune to generate solutions (which are saved as PDB files).

Experiments were conducted on seven small proteins and peptides from the Protein Data Bank [190]: 2jmy, 1vm2, 2jp8, 2jta, 6nm2, 6nm3 and 2fbu. The PDB contains known structures of these proteins, which were generated by NMR in combination with one of the simulated annealing methods described in Section 3.2.3. For example, the known structure for the peptide 2jp8 was created with the software XPLOR-NIH [185]. Table 3.1 shows the number of atoms ($|V|$) and the number of inter-atomic distances available in the generated DDGP instances, grouped by distance type (which sum to $|E|$). Note that the values in this table differ from the ones reported in Table 1 in [27]. This is because, for these experiments, we only consider the backbone, while in [27] the side chains were also considered. The table shows the different properties of the proteins. First, we see the size of the proteins in terms of the number of atoms as well as the number of amino acids. Next, the different number of distances are listed. Only a subset of the van der Waals distances are included, only those between a select number of C atoms. This is because including too many distances slows down the refinement step, and we do not need to consider all type 1 distances in order to avoid clashes.

<footnote n="2">2. https://github.com/simonheng/BP_ProteinFileReader</footnote>
<table><caption>Table 3.1 – Different properties of the proteins and peptides considered for these experiments. Includes the number of distances, per type, that can form our DDGP instances.</caption><thead><tr><th>Properties</th><th>2JMY</th><th>1VM2</th><th>2JP8</th><th>2JTA</th><th>6NM2</th><th>6NM3</th><th>2FBU</th></tr></thead><tbody><tr><td>Number of amino acids</td><td>15</td><td>13</td><td>7</td><td>10</td><td>8</td><td>8</td><td>12</td></tr><tr><td>Number of atoms</td><td>91</td><td>78</td><td>43</td><td>61</td><td>48</td><td>48</td><td>73</td></tr><tr><td>Type 1 (van der Waals radii)</td><td>91</td><td>66</td><td>21</td><td>36</td><td>21</td><td>21</td><td>55</td></tr><tr><td>Type 2 (Bonds/bond angles)</td><td>350</td><td>316</td><td>171</td><td>61</td><td>191</td><td>190</td><td>296</td></tr><tr><td>Type 3 (NMR restraints)</td><td>59</td><td>37</td><td>10</td><td>12</td><td>16</td><td>8</td><td>14</td></tr><tr><td>Type 4 (NMR torsion)</td><td>20</td><td>10</td><td>0</td><td>0</td><td>5</td><td>6</td><td>9</td></tr><tr><td>Type 3 Average range (Å)</td><td>2.07</td><td>2.42</td><td>1.72</td><td>3.80</td><td>2.18</td><td>2.60</td><td>1.67</td></tr><tr><td>Type 4 Average range (°)</td><td>68.2</td><td>80.0</td><td>-</td><td>-</td><td>72</td><td>113.7</td><td>65.55</td></tr></tbody></table>

For the first peptide, 2JMY, we see that there are 15 amino acids. This means that there are 14 $\phi$ and 14 $\psi$ backbone dihedral angles for this protein. However, in the table, we see that the NMR chemical shift prediction only gave us information about 20 out of 28 torsion angles. This is a recurring phenomenon for these peptides. For the two peptides 2JP8 and 2JTA we have no dihedral angle information whatsoever. This is logical because these structures have no clear secondary structure present, and thus TALOS+ is not able to predict these backbone dihedral angles. To still discretize despite missing this information, for each of the absent dihedral angles we include an interval from -180° to 180°, covering the full range and both signs. Missing backbone torsion angles is one part of the uncertainty present in these instances. Aside from this lack of information, we have very large degrees of noise, also visible in the information presented in the table. The intervals on the NMR restraints (type 3) can sometimes be quite large. Even more, the backbone torsion angles predicted with TALOS+ [183] tend to have average ranges of
65° up to even 113°. The uncertainty in these backbone torsion angles means that the search space of the tree will be quite large.

For this experiment, we used the force field CHARMM and ran Java MDJEEP with tolerance parameter $\epsilon$ = 0.01Å and the maximum number of solutions was capped at 5000. The resolution parameter $\rho$ was set to 2.0Å. When we decrease $\rho$, we cut the sub-arcs more often, thus increasing the chances of finding even better structures, at the expense of increasing the size of the tree. Recall that the effect of the resolution parameter $\rho$ is twofold, not only does it decide the length of the sub-arcs we create when we branch, but it is also used to decide when two consecutive solutions are too close to each other. In fact, because the uncertainty in the instances causes a very large search space, many of these solutions are pruned using this resolution parameter. The value of 2.0Å was chosen for these proteins after preliminary tests because it gives a good balance between the quality of the best solution found and the total size of the tree. To measure the quality of the accepted solutions, we compute the Root mean square deviation (RMSD) between the computed structures and the known structures, after alignment. This straightforward metric is computed for two realizations $x$ and $y$ as follows:

$$ \text{RMSD}(x, y) = \sqrt{\frac{1}{|V|} \sum_{v \in V} \|x_v - y_v\|^2}, \quad (3.1) $$

where $|| \cdot ||$ is the Euclidean distance norm. For aligning we used a Kabsch alignment algorithm [193]. Our implementation also attempts, together with translations and rotations, to perform a total reflection of the protein models to improve the alignments. The results of the experiments are summarized in Table 3.2.

<table><caption>Table 3.2 – The resulting RMSD scores for the proteins (in Å).</caption><thead><tr><th>Results</th><th>2JMY</th><th>1VM2</th><th>2JP8</th><th>2JTA</th><th>6NM2</th><th>6NM3</th><th>2FBU</th></tr></thead><tbody><tr><td>Best RMSD</td><td>1.63</td><td>2.02</td><td>1.55</td><td>2.00</td><td>1.66</td><td>1.49</td><td>2.52</td></tr><tr><td>Worst RMSD</td><td>3.61</td><td>6.88</td><td>4.67</td><td>5.98</td><td>6.46</td><td>5.06</td><td>5.11</td></tr><tr><td>Mean RMSD</td><td>3.0</td><td>4.43</td><td>2.89</td><td>4.02</td><td>2.78</td><td>3.47</td><td>4.06</td></tr></tbody></table>

Firstly, we see that the best found RMSD scores of the computed structures for each of the proteins are quite low. This means that some of the structures that we compute are close to the structures that are in the PDB, which is a very good sign. Secondly, we can
note that the differences between the worst and best RMSD scores are quite large. This is
no surprise, given the large amount of uncertainty in the instances. This uncertainty leads
to a very large search space, which means there are many possible solutions. Furthermore,
each of these possible solutions satisfy all the distance and angle restraints at hand, given
the tolerance ϵ, which means they are all valid solutions for the DGP at hand. We will
now take a closer look at the computed structures with the best (lowest) RMSD scores.

<figure><img src="image_15.png" /><figcaption>Figure 3.3 – The structures computed by Java MDjeep (blue) with lowest RMSD super-imposed with the structures from the PDB. The structures were not subject to energy minimization or any other post-processing step.</figcaption></figure>
In Figure 3.3 we see the computed structures super-imposed with one of the measured (ground truth) structures from the PDB. Keep in mind that the structures that we see are the raw output of Java MDJEEP and were not subject to energy minimization. Visually, the results for some proteins are better for others, which seem to match the scores reported in the table. For instance, the three proteins 2JMY, 6NM2 and 6NM3 had the lowest RMSD scores, and their computed structures look fairly close to the measured structure. The structure computed for 2FBU has the highest RMSD score among the proteins in the testbed, and this is also clearly visible in the figure. The helix in the computed structure does not look like a proper helix, which means that some of the backbone dihedral angles present in the structure are not correct. One interesting result is that in our computed structure, there seems to be a secondary structure element present in the protein 2JTA (visualized as a thicker ribbon) which is not present in the measured structure of the protein. From these experiments, we can conclude that for small proteins, branch-and-prune (with refinement) is a useful tool and is capable of computing good-quality structures when using real protein data. For larger proteins, there are still some obstacles to be overcome, as we will see in the next section.

### 3.5.2 Covalent geometry and branch-and-prune

In the previous section, we saw that the uncertainty in genuine NMR instances can cause challenges for branch-and-prune. The obstacles discussed relate to interval distances and torsion angles which are **type 3** and **type 4**. However, in practice, there can also be a problem with the distances that are often regarded by BP methods as exact. These are the distances from force field parameters (**type 2**), which are bond lengths, bond angles as well as the $\omega$ torsion angle (which is often set to either 0° or 180°). In practice, these distances and angles vary a lot from what is given by the force field parameter files. This issue was studied in [134], where we compared IBP-NG with an intermediate version of Java MDJEEP. Through the use of the refinement step with SPG, MDJEEP was able to identify solutions for all the proteins in the testbed, while IBP-NG was not. However, IBP-NG on average computed structures with better RMSD scores. Since then, Java MDJEEP was further updated (using torsion angles for branching), which greatly improved the quality of the computed structures and which lead to much better RMSD values. In this section, we will extend the experiments presented in [134], using the newest version of Java MDJEEP. These experiments aim to see the effect of the variations of the covalent geometry on the structures produced by branch-and-prune. In the previous section, we
conducted experiments with proteins for which the structures in the PDB were produced
by NMR and simulated annealing. In this section, we will instead use proteins that were
analyzed with X-ray crystallography. In particular, we selected 435 proteins from the
PISCES server [194]. The criteria for selection were:

- Smaller than 100 residues

- No Proline present (all $\omega$ angles should be close to 180°)

- An identity smaller than 20% between the structures (they are not too similar)

- An X-ray crystallographic resolution better than 1.6Å (the accuracy of the X-ray measurements)

- An R-factor of at least 0.25 (how well the refined structure predicts observed data)

We start by presenting a small data analysis with these proteins which confirms that there are substantial variations present in the bond angles as well as in the $\omega$ torsion angles. This analysis was done on a subset of the proteins in the testbed. More details can be found in [134]. The analysis is summarized in Figure 3.4.

Firstly, we can see that the angles can vary quite strongly between amino acid types, which is not captured by the force field data (dashed line). We also see the deviations based on some Ramachandran regions of the proteins. These regions are defined by combinations of values that the $\phi$ and $\psi$ angles often take [195]. The twelve Ramachandran regions used in this analysis are further illustrated in Figure 3.5, which displays an approximation of the Ramachadran plot. Recall that the $\omega$ angles are generally assumed to be close to 180°. However, we see in Figure 3.4 that for several of these regions, it is especially striking how much the $\omega$ angles can deviate from 180°. In some cases it can even be around 150°. These variations can be a problem for the BP methods. For small proteins, such as those analyzed in the previous section, this may not be a problem. However, when we look at larger proteins, the error introduced by these variations can build up. The proteins that we experiment with in this section are up to ten times larger than the NMR proteins that we analyzed before, and will thus likely exhibit problems with the variations of covalent geometry.
<figure><img src="image_16.png" /><figcaption>Figure 3.4 – The average values of the bond angles between the heavy backbone atoms (N, CA, C, O) as well as on the torsion angle ω of the peptide plane. Np denotes the nitrogen atom N of the next residue, Cm denotes the carboxyl carbon C of the previous residue, and CA denotes the CA atom. The distributions of these angle values are uni-modal, which means that we can use their averages as a suitable descriptor. On the left, we see the averages of the angles according to the amino-acid types, while on the right we see them according to the Ramachandran region of the residue. The horizontal dashed lines correspond to the values of the angles from the force field parameters [180].</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 3.5 – The Ramachandran plot illustrating the 12 different regions used in the data analysis. Note that this figure is only an approximation of the regions based on a 5-degree grid.</figcaption></figure>

To study how these variations can affect branch-and-prune, we ran experiments on two sets of instances, that will be referred to as instances A and B. Both these instances are generated using the same instance generation tool described in the previous section, which can be found on GitHub. The parsing of PDB files is done by the PDB class in the main Java package (see Figure 1.4). Instance A is a completely artificial instance, where we generate distances purely from the PDB files, introducing only a small amount of error. Instance B is an instance where we use genuine (**type 2**) distance information from the force field parameters and set $\omega$ to a small interval around 180°. Note that
since we are dealing with structures from X-ray crystallography, we do not have NMR restraints (**type 3**) in our instances because there are no hydrogen atoms present in the PDB structures. In fact, we should also not have access to torsion angles and distances relating to (**type 4**), because they are generally obtained by making predictions based on the chemical shifts from NMR spectra. However, we require these $\phi$ and $\psi$ torsion angles for the discretization of the instances, which is why they are simulated from the measured structures in the PDB files (with the introduction of some noise). In fact, for instance class B, the only artificial information that is included are these (**type 4**) torsion angles. The instances are further detailed in Table 3.3.

<table><caption>Table 3.3 – Details for the two instance classes that were experimented with. Instances A are completely artificial while instance B mixes artificial with information from force-field parameters.</caption><thead><tr><td>Distance/angle type</td><td>Instances A</td><td>Instance B</td></tr></thead><tbody><tr><td>Type 1 (van der Waals radii)</td><td>60% of sum of radii</td><td>60% of sum of radii</td></tr><tr><td>Type 2 (bond lengths and bond angles)</td><td>Exact distances and ω angle taken from PDB (capturing the variations)</td><td>From force field parameters, ω in a 5 degree interval around 180°</td></tr><tr><td>Type 2 (Improper torsion angles)</td><td>From force field parameters</td><td>From force field parameters</td></tr><tr><td>Type 3 (NMR restraints)</td><td>Not present</td><td>Not present</td></tr><tr><td>Type 4 (TALOS+ torsion angles)</td><td>From PDB placed randomly in 5° intervals</td><td>From PDB placed randomly in 5° intervals</td></tr></tbody></table>

There are some important details to note about these instance types. Recall that the geometric information captured by **type 2** and **type 4** is used for discretization, and that their distances are used as reference distances in branch-and-prune. For these instances, for every amino acid (except the very first one), we have access to the proper torsion angles $\phi$, $\psi$, and $\omega$ and we know the improper dihedral angles from the force fields. This means that whenever we are branching, for almost every candidate atom $v$ we know the dihedral angle $\tau$ ($\{u_3, u_2, u_1, v\} \in E_{\tau}$). Therefore, we will only branch based on the size
of the intervals of these torsion angles, which we know is 5°. Thus, the BP trees that are created for these instances have a small size. For these instances, with our parameters, Java MDJEEP will only identify one solution because the other solutions are too close and are thus rejected. These instances are generated in this way specifically because we want to study the effect of the variations of the covalent geometry.

For each of the 435 proteins in the testbed, we created instances A and B and used Java MDJEEP to find one solution for every instance. The parameters used were $\epsilon = 0.01$, $\rho = 1.0$ and the force field used was CHARMM [179]. Similarly to the experiments in the previous section, to measure the quality of the computed structures we use the RMSD metric (see Equation (3.1)). The results are visualized in Table 3.4 and further explored in Figure 3.6.

<table><caption>Table 3.4 – The resulting RMSD statistics for the 435 proteins for instances A and B.</caption><thead><tr><th>Results</th><th>Instance A</th><th>Instance B</th></tr></thead><tbody><tr><td>Mean</td><td>0.79</td><td>3.45</td></tr><tr><td>Median</td><td>0.65</td><td>3.37</td></tr><tr><td>Standard deviation</td><td>0.62</td><td>1.78</td></tr><tr><td>Worst RMSD score</td><td>5.64</td><td>11.69</td></tr><tr><td>Correlation coefficient size and RMSD</td><td>0.49</td><td>0.80</td></tr></tbody></table>

The results for the instances of class A are very good, with a mean of 0.79Å, which is a rather low score for proteins of these sizes. In fact, as we see in Figure 3.6, most of the proteins are in the first bin with RMSD scores between 0 and 2. Logically, for instance class B, the results are much worse, with much higher RMSD scores with a mean of 3.37Å. This was to be expected because the instances of class A rely mostly on data extracted directly from the X-ray structure that we are comparing with. The last line in Table 3.4 is also of interest, which displays the correlation coefficient between the RMSD scores and the size of the protein. We can see that the correlation between size and RMSD is much stronger for instance class B (0.80) than for A (0.49). This confirms the idea that there is a large build-up of errors when we do not account for variations in the covalent geometry. In instance A, there is still an intermediate correlation between the size of the protein and the RMSD. One of the reasons for this is that the uncertainty in the $\phi$ and $\psi$ angles
will also build up, but to a lesser degree.

<figure><img src="image_18.png" /><figcaption>Figure 3.6 – A histogram with the resulting RMSD scores for the two instance classes, with bins of size 2Å.</figcaption></figure>

In Figure 3.7 we see four example proteins from the testbed. For each of the four proteins, the two generated structures with instances A and B are overlapped with the original, measured PDB structure. We see that as the RMSD scores suggest, the pink structures (instance A) are very close to the original structures, while the blue (instance B) one is not. In particular, for the protein 2ZW2 (which has the worst RMSD score of all instance B proteins), the blue structure deviates quite significantly. For the protein 5YIU, both structures that were generated look quite close to the original. This protein has 45 amino acids and is thus one of the smaller ones in the testbed, which can be one of the reasons why it was easier to solve. The only difference between the instance classes is the covalent geometry (**type 2**), which for instances A captures the variations, and for instances B is static. These much worse RMSD scores confirm that the variations in the covalent geometry pose a problem for branch-and-prune. However, despite these variations, Java MDJEEP was able to find a solution for every protein in the testbed. Considering the small size of the search space, this is a good result, which is most likely due