NMR assignment through linear programming

Abstract

Nuclear Magnetic Resonance (NMR) Spectroscopy is the second most used technique (after X-ray crystallography) for structural determination of proteins. A computational challenge in this technique involves solving a discrete optimization problem that assigns the resonance frequency to each atom in the protein. This paper introduces LIAN (LInear programming Assignment for NMR), a novel linear programming formulation of the problem which yields state-of-the-art results in simulated and experimental datasets.

Appendices

A: Grouping peaks

As we mentioned in Sect. 3.1.1, grouping consistent peaks together is a crucial step in the graph creation process for ${\mathcal {G}}=({\mathcal {V}},{\mathcal {E}})$. One would wish the enumeration of valid assignments to be as thorough as possible. We can effectively enumerate peak groupings to construct nodes in ${\mathcal {G}}$ by matching measured and expected peaks in a self-consistent way. In particular, we expect a specific set of peaks due to N–H${}^{N}$ from residue k (see Fig. 10 for a standard example with three experiments) where the values of these peaks in ${\mathbb {R}}^3$ along certain dimensions are consistent. If there are n residues, we should have n sets of such expected peaks. Therefore, each layer in ${\mathcal {G}}=({\mathcal {V}},{\mathcal {E}})$ in principle should have n nodes, although in practice there are more nodes due ambiguities.

Fig. 10

With three NMR experiments (often HSQC, HNCACB, and HN(CO)CACB) we generally expect 7 distinct peaks for each base N–H${}^{N}$ pair in a residue, k. These peaks must be consistent—that is, the frequencies assigned to the same atom by two different peaks must be approximately the same up to some experimental tolerance. In principle, there should be n sets of such 7 peaks, one for each residue

The notion of consistency can help significantly simplify the enumeration process (which would otherwise result in an exponential number of nodes). In order to efficiently enumerate consistent peak groupings, we do the following. Let $\mathcal {S}_1, \ldots , \mathcal {S}_{L}$ be collections of measured peak lists corresponding to different heteronuclear experiments, i.e. $\cup _{l=1}^L\mathcal {S}_l:=[p_1, \ldots , p_{m_2}]$. In the case of Fig. 10, $L=3$, as we have peaks from three experiments. Now from these $m_2$ experimental peaks we form all combinations of seven peaks that each consists of one peak from ${\mathcal {S}}_1$, two peaks from ${\mathcal {S}}_2$, and four peaks from ${\mathcal {S}}_3$ using the following criteria.

• For any pair of $p_u, p_v$ in a combination of seven peaks,

  $$\begin{aligned} \vert p_u(1)-p_v(1)\vert&\le \delta _1 \\ \vert p_u(2)-p_v(2)\vert&\le \delta _2. \end{aligned}$$

  This means that the frequencies of the seven peaks in the N–H${}^{N}$ dimension have to coincide up to tolerance $\delta _1,\delta _2$.

• Furthermore, for a combination of seven peaks, let $p_u, p_v$ be the two peaks in ${\mathcal {S}}_2$. These peaks should coincide with two of the peaks in ${\mathcal {S}}_3$ (denoted $p_i,p_j$) up to tolerance $\delta _3$, i.e.

  $$\begin{aligned} \vert p_u(3)-p_i(3)\vert&\le \delta _3 \\ \vert p_v(3)-p_j(3)\vert&\le \delta _3 \end{aligned}$$

  along the $\text {C}$ dimension.

B: Atom cost

Recall that we defined the cost of an atom, a, under a given set of assigned observations, $\{x_l\}_{l=1}^{o_a}$ as

Definition 3

(Atom cost) The cost associated with atom a, with a normally distributed prior $\mathcal {N}(\mu _a, \sigma _a)$, and $o_a$ observations $\{x_l^a\}_{l=1}^{o_a}$ defined by the peak grouping, also assumed to be normally distributed around the true frequency, $\mu$, according to $\mathcal {N}(\mu , \sigma _l)$ is defined as

$$\begin{aligned} \text {cost}\left( a, \{x_l^a\}_{l=1}^{o_a}\right) \triangleq -\log {\mathbb {E}}_{\mu \sim {\mathcal {N}}(\mu _a, \sigma _a)}\left[ \prod _{l=1}^{o_a}f(x_l^a\mid \mu , \sigma _l)\right] . \end{aligned}$$

(18)

where $f(\cdot \mid u, v)$ is the Gaussian density with mean u and standard deviation v.

This is Definition 1 in the main text. Note that the term inside the expectation is a product of $o_a$ univariate Gaussian probability density functions. Furthermore, expanding the expectation, we note that

$$\begin{aligned} {\mathbb {E}}_\mu \left[ \prod _{l=1}^{o_a}f(x_l^a\mid \mu , \sigma _l)\right]&= \int _{-\infty }^{+\infty }f(\mu \mid \mu _a, \sigma _a)\prod _l^{o_a}f(x_l^a\mid \mu , \sigma _l)d\mu \end{aligned}$$

(19)

$$\begin{aligned}&=\int _{-\infty }^{+\infty }f(\mu \mid \mu _a, \sigma _a)\prod _l^{o_a}f(\mu \mid x_l^a, \sigma _l)d\mu \end{aligned}$$

(20)

by symmetry. Using a standard result regarding the product of univariate Gaussian PDFs (see, e.g., [11]), we can write

$$\begin{aligned} {\mathbb {E}}_\mu \left[ \prod _{l=1}^{o_a}f(x_l^a\mid \mu , \sigma _l)\right]&=\int _{-\infty }^{+\infty }f(\mu \mid \mu _a, \sigma _a)\prod _l^{o_a}f(\mu \mid x_l^a, \sigma _l)d\mu \end{aligned}$$

(21)

$$\begin{aligned}&=\int _{-\infty }^{+\infty }Z_af(\mu \mid M_a, \Sigma _a)d\mu \end{aligned}$$

(22)

$$\begin{aligned}&=Z_a \end{aligned}$$

(23)

where

$$\begin{aligned} \Sigma _a&= \left( \frac{1}{\sigma _a^2}+\sum _{l=1}^{o_a} \frac{1}{\sigma _l^2}\right) ^{-1/2} \end{aligned}$$

(24)

$$\begin{aligned} M_a&= \left( \frac{\mu _a}{\sigma _a^2}+\sum _{l=1}^{o_a}\frac{x_l}{\sigma _l^2}\right) \Sigma ^2_{a} \end{aligned}$$

(25)

$$\begin{aligned} Z_a&=\frac{1}{(2\pi )^{o_a/2}}\sqrt{\frac{\Sigma _a^2}{\sigma _a^2\prod _{l=1}^{o_a}\sigma _l^2}}\exp \left[ -\frac{1}{2}\left( \frac{\mu _a^2}{\sigma _a^2}+\sum _{l=1}^{o_a}\frac{x_l^2}{\sigma _l^2}-\frac{M_a^2}{\Sigma _a^2}\right) \right] . \end{aligned}$$

(26)

We see that this choice of cost function is therefore computationally advantageous, as the desired expectation is a simple function of the observations, $\{x_l\}_{l=1}^{o_a}$ and of the distributional parameters of the prior, $(\mu _a, \sigma _a)$ and experiments, $\{\sigma _l\}_{l=1}^{o_a}$. That said, it is certainly not the only cost function that one could use. As an example, we could instead solve a maximum likelihood problem for each peak grouping that would assign the highest likelihood frequency to each atom, given the prior and the observations. The exploration of alternative cost functions is left for future work.

C: Statistical Typing

Statistical typing is a process that happens both during the node and edge creation steps. In particular, we want to avoid the creation of nodes and edges which are too unlikely to constitute a valid assignment. The way we action on this notion is to define a threshold below which we would rather have a null assignment than the assignment induced by the relevant nodes. This threshold also determines the cost of the edges to (and from) the dummy nodes, which are therefore the highest cost edges in the graph.

For all simulations in this paper, we use the following definition:

Definition 4

(Atom cost threshold) The maximum allowable cost associated with atom a, with an expected frequency, $\mu$, distributed according to the normally distributed prior $\mathcal {N}(\mu _a, \sigma _a)$, and a total of $o_a$ expected observations is given by:

$$\begin{aligned} \text {threshold}\left( a\right) \triangleq \text {cost}(a, \{w_l^a\}_{l=1}^{o_a}) \end{aligned}$$

(27)

where

$$\begin{aligned} w_l^a = \mu _a +\delta \sigma _a + (-1)^{l+1}\delta \sigma _l. \end{aligned}$$

(28)

That is, we define the maximum allowable cost for atom a by setting $\{x^a_l\}_{l=1}^{o_a}$ in Definition 1 to $\{w^a_l\}_{l=1}^{o_a}$, which constitute an adversarial realization of the observations. In this realization, the mean of the observations is $\approx \delta$ standard deviations away from the prior mean, and the observations are split into two clusters, $2\delta$ experimental standard deviations apart.