Solving the molecular distance geometry problem with inaccurate distance data

The knowledge of the protein structure is very important to understand its function and to analyze possible interactions with other proteins. Different methods can be applied to acquire protein structural information. Until 1984, the X-ray crystallography was the ultimate tool for obtaining information about protein structures, but the introduction of nuclear magnetic resonance (NMR) as a technique to obtain protein structures made it possible to obtain data with high precision in an aqueous environment much closer to the natural surroundings of living organism than the crystals used in crystallography [1].

The NMR technique provides a set of inter-atomic distances for certain pairs of atoms of a given protein. The molecular distance geometry problem (MDGP) arises in NMR analysis context. The MDGP consists of finding one set of atomic coordinates such that a given list of geometric constraints are satisfied [2]. Formally, the molecular distance geometry problem can be defined as the problem of finding Cartesian coordinates $x_{\mathsf{\text{1}}},...,x_n\in {ℝ}^{\mathsf{\text{3}}}$ of atoms of a molecule such that l _ij ≤ ||x _i - x _j || ≤ u _ij , ∀(i, j) ∈ E, where the bounds l _ij and u _ij for the Euclidean distances of pairs of atoms (i, j) ∈ E are given a priori [3].

Clearly, $x=\left(x_{\mathsf{\text{1}}},...,x_n\right)\in {ℝ}^{\mathsf{\text{3}}n}$ solves the MDGP if, and only if, x is a global minimizer of f and f(x) = 0.

An overview on methods applied to the MDGP is given in [4] and a very recent survey on distance geometry is given in [5].

where $i\in \left\{i_1,i_2,i_3,i_4\right\}\subset \mathcal{B}$ and d _ij = ||x _j - x _i ||. The indexes i₁, i₂, i₃, i₄ can be chosen in an arbitrarily way, allowing us to choose another base subset when calculating the coordinate of the next x _j . At each iteration, the index j of the new coordinate x _j is inserted in the set $\mathcal{B}$ increasing the number of subsets {i₁, i₂, i₃, i₄} used as anchors to fix the remaining unset coordinates.

Unfortunately, in practice, the NMR experiments just provide a subset of distances between atoms spatially close and the data accuracy is limited. Thus in the real scenario, the set E is sparse and l _ij < u _ij . So, we just have bounds to some of the entries of the distance matrix D. In this situation, neither the singular value decomposition nor the buildup algorithm can be applied directly because they are both designed to deal with exact distances. In fact, the inaccurate and sparse instances of MDGP, where l _ij <u _ij , are much harder to solve as pointed by Moré and Wu who showed that the MDGP with inaccurate distances belongs to the NP-hard class of problems [7].

Our contribution is a new algorithm that can handle with inaccurate and sparse distance data. We propose an iterative method based on simple ideas: generate an approximated distance matrix D, take as base a clique in the graph that has D as a connectivity matrix, solve the system (1) and refine the solution using a nonlinear least-squares method. It needs to be pointed that the authors of the buildup algorithm and coworkers have done some modifications in the original form of the algorithm in order to handle inaccurate data [8, 9]. However, the main advantage of our proposal is its simplicity and robustness. We have been able to find solutions with acceptable quality to instances of MDGP with inaccurate and sparse data, considering up to thousands of atoms.

We have implemented our algorithm lsbuild in Matlab and tested it with a set of model problems on an Intel Core 2 Quad CPU Q9550 2.83 GHz, 4GB of RAM and Linux OS-32 bits. In all experiments the parameters of the function φ _λ,τ of the algorithm lsbuild were set at λ = 1.0 and at τ = 0.01.

The algorithm dgsol starts with an approximated solution and, given a sequence of smoothing parameters λ₀ > λ₁ > ... > λ _p = 0, it determines a minimizer x_k+1of 〈f〉_λ. The algorithm dgsol uses the previous minimizer x _k as the starting point for the search. In this manner a sequence of minimizers x₁, ..., x_p+1is generated, with the x_p+1a minimizer of f and the candidate for the global minimizer. In our experiments, we used the implementation of the algorithm dgsol encoded in language C and downloaded from [23].

with variables $x_i=\left(x_i^1,x_i^2,x_i^3\right),x_j=\left(x_j^1,x_j^2,x_j^3\right)\in {ℝ}^3$ and B being the set formed by the index k and the indexes of all neighbors of x _k in the current base set. The parameters d _kj are the given distances between the node x _k and its neighbors x _j in the base and, for the nodes x _j and x _i already in the base, if the distance between them is unknown, we consider d _ij = ||x _i - x _j ||. Once the substructure is obtained, it is inserted in the original structure by an appropriated rotation and translation and the point x _k is included in the base. This process is repeated until all nodes are included in the base. We have implemented the buildup algorithm in Matlab.

Our decision to compare the lsbuild with the algorithms dgsol and buildup is mainly motivated by theirs similarities with our proposal. In fact, the algorithm dgsol uses a smooth technique in order to avoid the local minimizers and the algorithm buildup solves a sequence of systems which produce partial solutions and iteratively try to construct a candidate to global solution. Our algorithm combines some variations of these two ideas. We use a hyperbolic smooth technique to insert differentiability in the problem and a divide-and-conquer approach based in sucessive solutions of overdetermined linear systems in order to construct a candidate to global solution.

In our experiments, the distance data were derived from the real structural data from the Protein Data Bank (PDB) [24]. It needs to be pointed that each of the algorithms considered has a level of randomness, the algorithm dgsol takes random start point and the algorithms lsbuild and buildup starts with an incomplete random matrix D = [d _ij ] where l _ij ≤ d _ij ≤ u _ij . So, in order to do a fair comparison, we run each test 30 times.

where R(k) represents the k-th residue.

where ${\widehat{d}}_{ij}$ is the real distance between the nodes x _i and x _j in the known structure x* of protein PDB:1GPV. In this way, all distances between atoms in the same residue or neighboring residues were considered. We generated two instances for each fragment by taking ε equals to 0.00 and 0.08.

is the error associated to the constraint l _ij ≤ ||x _i - x _j || ≤ u _ij : We also measured the deviation

with $h\in {ℝ}^{n\times 3}$ and $Q\in {ℝ}^{3\times 3}$ orthogonal.

where $d_{ij}^{*}$ is the true distance between atom i and atom j and ${\bar{\epsilon }}_{ij},{\underset{-}{\epsilon }}_{ij}~\mathcal{N}\left(0,{\sigma }_{ij}^2\right)$ (normal distribution). With this model, we generate a sparse set of constraints and introduce a noise in the distances that are not so simple as the one used in the instances proposed by Moré and Wu.

Finally, the results of both set of instances indicate that our algorithm lsbuild based on the combination of the resolution of linear systems, derived from the approximated EDM matrices, and the refinement process based on hyperbolic smoothing penalty is a very effective strategy to solve MDGP instances with sparse and inaccurate data.

We presented a new algorithm to solve molecular distance geometry problems with inaccurate distance data. These problems are related to molecular structure calculations using data provided by NMR experiments which, in fact, are not precise. Our algorithm combines the divide-and-conquer framework and a variation of the hyperbolic smoothing technique. The computational results show that the proposed algorithm is an effective strategy to handle uncertainty in the data.