arXiv:1606.06975v1 [] 22 Jun 2016

Bias Correction in Saupe Tensor Estimation
Y. Khoo∗

A. Singer†

D. Cowburn‡

June 23, 2016

Abstract
Estimation of the Saupe tensor is central to the determination of
molecular structures from residual dipolar couplings (RDC) or chemical shift anisotropies. Assuming a given template structure, the singular value decomposition (SVD) method proposed in [15] has been
used traditionally to estimate the Saupe tensor. Despite its simplicity,
whenever the template structure has large structural noise, the eigenvalues of the estimated tensor have a magnitude systematically smaller
than their actual values. This leads to systematic error when calculating the eigenvalue dependent parameters, magnitude and rhombicity.
We propose here a Monte Carlo simulation method to remove such
bias. We further demonstrate the effectiveness of our method in the
setting when the eigenvalue estimates from multiple template protein
fragments are available and their average is used as an improved eigenvalue estimator. For both synthetic and experimental RDC datasets of
ubiquitin, when using template fragments corrupted by large noise, the
magnitude of our proposed bias-reduced estimator generally reaches at
least 90% of the actual value, whereas the magnitude of SVD estimator
can be shrunk below 80% of the true value.

1

Introduction

1.1

Background

The residual dipolar couplings (RDC) of molecules can be measured when
the molecule ensemble in solution exhibits partial alignment with the mag∗

Department of Physics, Princeton University, Princeton, NJ 08544,USA
(ykhoo@princeton.edu).
†
Department of Mathematics and PACM, Princeton University, Princeton, NJ 08544,
USA (amits@math.princeton.edu).
‡
Department of Biochemistry, Albert Einstein College of Medicine, Bronx, NY 10461,
USA (david.cowburn@einstein.yu.edu).

1

netic field in a nuclear magnetic resonance (NMR) experiment. Due to the
r −3 dependence, such effects can be measured with high precision and provides accurate alignment information of a specific single or multiple bond
vector to the magnetic field. This is in contrast to the Nuclear Overhauser
Effect (NOE) which scales as r −6 . Hence the importance of RDC in obtaining high quality protein structures and studying molecular dynamics has
increased considerably over the last decade. For a detailed survey of RDC
and its applications we refer readers to [14, 1, 25, 19].
We first give a brief summary of RDC. Let vnm be the unit vector denoting the direction of the vector between nuclei n and m. Let b be the unit
vector denoting the direction of the magnetic field. The RDC Dnm due to
the interaction between nuclei n and m is

 T
2
max 3(b vnm ) − 1
Dnm = Dnm
.
(1)
2
t,e
max is a constant depending on the gyromagnetic ratios γ , γ of the two
Dnm
n m
nuclei, bond length rnm , and Planck’s constant h as
max
Dnm
=−

γn γm h
,
3
2π 2 rnm

(2)

and h ·it,e denotes the ensemble and time averaging operator. As presented,
RDC depends on the relative angle between the magnetic field and the
bond. Extracting such angular information from RDC complements NOE
and other measurements for determining the molecular structure. It is conventional to interpret the RDC measurement in the molecular frame. More
precisely here, for this analysis of bias, we treat the molecule as being static
in some coordinate system, and the magnetic field direction being a time
and sample varying vector. In this case the RDC becomes
max T
Dnm = Dnm
vnm Svnm ,

(3)

where the Saupe tensor S is defined as where the Saupe tensor S [20] is
defined as
1
B = bbT t,e .
(4)
S = (3B − I3 ),
2
B is known as the field tensor and I3 denotes the 3 × 3 identity matrix.
We note that S is symmetric and Tr(S) = 0. In order to use RDC for
structural refinement of a protein, S is usually either first determined from
a proposed structure (known vnm ) that is similar to the protein under study,
or estimated from a distribution of RDCs [4] . To satisfy the assumption
2

that the molecule is static in the molecular frame, a rigid fragment of the
known structure has to be selected. S can be determined if the fragment
contains sufficient RDC measurements.

1.2

Notation

We summarize here the notation that is used throughout the paper. For a
3 × 3 matrix A, we use Aij , i, j = x, y, z to denote the nine entries of the
matrix. When A is symmetric, we denote the eigen-decomposition of A by
A = U (A)Λ(A)U (A)T ,
where U (A) is an orthogonal matrix (i.e. U (A)T U (A) = U (A)U (A)T = I3 )
and Λ(A) is a diagonal matrix


λx (A)
0
0
 0
λy (A)
0 
(5)
0
0
λz (A)
that contains the eigenvalues of A on the diagonal in ascending order. For
a matrix A ∈ Rn×n we use
v
uX
u n
A2ij
kAkF = t
i,j=1

to denote the Frobenius norm of the matrix. For a vector v, vi denotes its
i-th entry, and i = 1, . . . , n if v ∈ Rn . In the special case of v ∈ R3 , we use
vx , vy , vz to denote components of the vector v. For a matrix A, we use Ai
to denote its i-th column.

1.3

Previous approach

We review the singular value decomposition (SVD) approach [15] for estimating the Saupe tensor. S is symmetric and Tr(S) = 0, so eqn. (3) can be
rewritten as
max
= (vnm 2y − vnm 2x )Syy + (vnm 2z − vnm 2x )Szz
Dnm /Dnm

+ 2vnmx vnmy Sxy + 2vnmx vnmz Sxz + 2vnmy vnmz Syz (6)
where vnmi , i = x, y, z are the different components of vnm in the molecular
max , which are the RDC measurements. When
frame. We let dnm = Dnm /Dnm
3

there are M RDC measurements, eq. (6) results in M linear equations in
five unknowns (Syy , Szz , Sxy , Sxz and Syz ), that can be written in matrix
form as


Syy


dn1 m1
Szz 




 ∈ R5 ,
As = d,
s=
S
(7)
d =  ...  ∈ RM ,
xy


Sxz 
dnM mM
Syz
and A ∈ RM ×5 .
Let the SVD of matrix A be
A = U ΣV T ,

(8)

where U ∈ RM ×5 is a column orthogonal matrix (i.e. U T U = I5 ), V ∈ R5×5
is an orthogonal matrix, and Σ ∈ R5×5 is a positive diagonal matrix. We
assume that M ≥ 5 and that A has full rank for otherwise there is no unique
solution to the linear system (7). The estimator of the Saupe tensor entries
s proposed in [15] is
ŝ = V Σ−1 U T d.
(9)
This is equivalent to the ordinary least squares (OLS) solution to the linear
system (7), given by
ŝ = (AT A)−1 AT d.
(10)
For this reason, we will refer to this SVD method for Saupe tensor estimation as the OLS method. The computational aspects of employing the
expressions in (9) and (10) are discussed in [29]. Notice that the Saupe
tensor estimator given by (9) and (10), denoted Ŝ, is the solution to the
optimization problem
min
S

M
X

|dni mi − vnTi mi Svni mi |2

such that S is symmetric, Tr(S) = 0. (11)

i=1

As such, the OLS estimator is also the maximum likelihood estimator when
the error on dnm is assumed to be white Gaussian noise.
Although the SVD procedure ensures Tr(S) = 0 and S symmetric, it does
not ensure that S is obligatorily derived from a specific linear combination
of the field tensor and the identity matrix as expressed in eqn. (4). In
contrast, we can use a positive semidefinite matrix description of the field
tensor B to exactly characterize S. A matrix is positive semidefinite (PSD)
if it is symmetric and has nonnegative eigenvalues. The field tensor B can
be characterized by the following observation:
4

Observation 1. B = hbbT it,e where b is a unit vector ⇔ B is positive
semidefinite and Tr(B) = 1.
This follows from the convexity of both the set of PSD matrices and the
set of unit trace matrices. A set is convex if and only if any weighted average
of the elements in the set belongs to the set. Since for any unit vector b, bbT
is PSD and Tr(bbT ) = Tr(bT b) = 1, the time and ensemble average of such
matrices is PSD with unit trace. Using this observation, the set of physical
Saupe tensors can be characterized as


1
(12)
S = S = (3B − I3 ) | B is PSD, Tr(B) = 1 .
2
However, since B = (1/3)(2S + I3 ), when Saupe tensor S has small entries
it is dominated by I3 and B is always positive semidefinite. This is often the
case in practice, hence the SVD method suffices for Saupe tensor estimation
at this approximation. We note that should we ever need to estimate S with
large entries of magnitude around 10−1 , we can solve
min
S∈S

M
X

|dni mi − vnTi mi Svni mi |2 ,

(13)

i=1

using semidefinite programming toolboxes available, e.g., in CVX [10] so that
the derived Saupe tensor remains physically reasonable.

1.4

An Alternative Approach

Here, we illustrate how the noise on the bond vectors vnm leads to bias
in the OLS estimation of Saupe tensor parameters. We call such type of
noise structural noise. In particular, we consider the situation when using
noisy template structures differing from the true structure of the molecule
for Saupe tensor fitting. This situation may arise when using homologous
structure in Saupe tensor fitting, or when the protein of interest has small
conformation changes due to the dynamic nature of protein in solution. Our
simulation consisted of adding modest noise to the backbone torsion angles
of the protein backbone. This results in the magnitude of the estimated
Saupe tensor eigenvalues being typically smaller than their true value, as
demonstrated in Fig. 1. Our observation corroborates with the simulation
results reported in [31], in which independently and identically distributed
(i.i.d.) random noise [28] is added to each bond vector instead. In linear regression, such decrease in magnitude of the estimator in the presence of noise
5

on the regressor is commonly known as attenuation [3]. While the focus of
[31] is mainly to use Monte Carlo simulation to evaluate the uncertainty of
estimated alignment magnitude and rhombicity, we here focus on using it
to correct the attenuation effect in the OLS Saupe tensor eigenvalues estimator. The method we propose bears similarity with the statistical method
simulation extrapolation (SIMEX) [5, 23] that is frequently used to correct
for the attenuation effect. Typically this type of methods are parametric
and require noise variance as input. We show that an estimator of the noise
magnitude can be obtained from the root mean square (RMS) of the residual
of OLS estimator. We further demonstrate the usefulness of removing such
bias when estimating the Saupe tensor eigenvalue from homology fragments
of ubiquitin, using RDCs measured in two different alignment medias. We
note that there are other approaches to improve the estimation of Saupe
tensor in the presence of structural noise by studying local bond orientations using multiple alignment media [16, 24, 17, 18]. However here, we seek
to remove the bias in the Saupe tensor eigenvalues in a single alignment
media, when multiple Saupe tensor estimates is available from a collection
of predetermined molecular fragments.

2

Method

We now introduce a Monte Carlo method for correcting the bias in the
eigenvalues of the OLS estimator arising from structural noise of backbone
torsion angles. For a protein with N + 1 peptide planes, we assume the
{φi , ψi }N
i=1 torsion angles fully determine the backbone conformation, i.e.
variations of {ωi }N
i=1 are minimal. The template structure’s torsion angles
t
t
φi , ψi ’s are related to the true structure via
φti = φi + σαi ,

ψit = ψi + σβi ,

i = 1, . . . , N

(14)

where αi , βi ’s are i.i.d. random normal variables with mean 0 and variance 1,
and σ is the level of noise on the torsion angles. Henceforth for a variable θ,
we make explicit the dependence on the torsion angles and noise by writing
θ as θ(φti , ψit ). We also assume that the normalized dipolar coupling d is
noiseless, i.e. d = A(φi , ψi )s, where s corresponds to the entries of the
”ground truth” Saupe tensor S. The validity of this assumption is discussed
below in section 4. Our method consists of the following steps:
(1) Compute
ŝ(φti , ψit ) = (A(φti , ψit )T A(φti , ψit ))−1 A(φti , ψit )T d
6

.
(2) Generate n1 copies of Asim = A(φti + σαi , ψit + σβi ) by adding i.i.d.
Gaussian noise with variance σ 2 to the torsion angles of the template structure.
(3) Find
ŝsim = ŝ(φti + σαi , ψit + σβi ) = (ATsim Asim )−1 ATsim d.
(4) Let Ŝ and Ŝsim be the Saupe tensor estimators corresponding to ŝ
and ŝsim . Let
d = hΛ(Ŝsim )isim − Λ(Ŝ)
Bias

denote the bias estimate for the eigenvalues of the OLS estimator Ŝ, where
h·isim denotes the averaging over n1 simulated template structures. We
propose deriving Λ (Eqn 5) from
d = 2Λ(Ŝ) − hΛ(Ŝsim )isim
Λ̃ = Λ(Ŝ) − Bias

as an estimator with less bias.
The rationale relies on the notion that upon adding noise of similar
magnitude to the linear system (7), the eigenvalues of the OLS estimator
for the simulated samples should be biased away from Λ(Ŝ) by an amount
similar to the difference between Λ(Ŝ) and the true Λ(S) This is also the
intuition behind twicing [26, 9], and related bootstrapped [8] biased reduced
estimators. Alternatively, one can understand this procedure from the viewpoint of the SIMEX technique [5] for correcting bias resulting from regressor noise. Under the SIMEX estimation framework one would simulate
Asim = A(φti + kσαi , ψit + kσβi ) with noise magnitudes of kσ for various
positive k to find out the dependency of Λ(Ŝsim ) on k. The k = 0 point corresponds to the case when no additional simulated noise is added, i.e. when
the eigenvalue estimator is Λ(Ŝ). From the extrapolation of the relation
between Λ(Ŝsim ) and k one can obtain a debiased estimator at k = −1. Our
method corresponds to the special case of SIMEX where we only add simulated noise with magnitude kσ where k = 1. Our numerical results shows
that this suffices for the application of Saupe tensor eigenvalue estimation.

2.1

Estimating σ

We note that there is a general caveat when using any parametric Monte
Carlo method, in that it requires knowledge of the noise magnitude σ. Let
the residual of the OLS estimator be defined as
r ≡ d − Aŝ.
7

2
is added to
In the simple case when additive noise with variance σadd
the normalized dipolar couplings d, and A has no structural noise, i.e.
A = A(φi , ψi ), the dependence between the RMS of the residual, denoted
RMS(r) and the noise magnitude can be readily calculated. In particular,
2
an unbiased estimator of σadd
is given by [11]
2
σd
add =

M
RMS(r)2 .
M −5

where M is the number of linear equations.
Now in the case when there is noise on the design matrix A = A(φti , ψit )
due to noise on the torsion angles (14), we show that there exists a linear
dependence of RMS(r) on σ. We define A0 = A(φi , ψi ), and A(φti , ψit ) =
A0 + E. In this notation, normalized RDC d = A0 s. Then
krk22 = kd − Aŝk22
= kA0 s − A(AT A)−1 AT (A0 s)k22
= sT AT0 (IM − A(AT A)−1 AT )A0 s.

(15)

The second equality follows from the fact that IM − A(AT A)−1 AT is a projection matrix. From
=
=
=
=

AT0 (IM − A(AT A)−1 AT )A0
AT0 A0 − (A − E)T A(AT A)−1 AT (A − E)
AT0 A0 − AT A + E T A + AT E − E T A(AT A)−1 AT E
AT0 A0 − (A − E)T (A − E) + E T E − E T A(AT A)−1 AT E
E T (IM − A(AT A)−1 AT )E,

(16)

we get
krk22 = sT E T (IM − A(AT A)−1 AT )Es
≈ sT E T (IM − A0 (AT0 A0 )−1 AT0 )Es
= sT E T P Es

(17)

where P = IM − A0 (AT0 A0 )−1 AT0 is a projection operator projecting vectors
in RM to RM −5 . We drop the terms involving entries of E raised to the
power greater than 2 to obtain the approximation in (17). Using Taylor
expansion,
Eij

= Aij − A0ij
N
X
∂Aij (φtk , ψkt )
≈
∂φtk
k=1

σαk +
φk ,ψk

8

∂Aij (φtk , ψkt )
∂ψkt

σβk
φk ,ψk

= Fij σ.

(18)

Plugging this into (17), it is clear that krk22 depends linearly on σ 2 and
hRMS(r)2 iαi ,βi ≈

1 T T
hs F P F siαi ,βi σ 2
M

(19)

in the small noise regime. We therefore use
r
M
RMS(r)
σ̂ =
T
s FTPFs
as the approximate noise magnitude when using the Monte Carlo method
for bias reduction. Although we do not have the parameters s, F and P
derived from the ground truth Saupe tensor and conformations, we can use
ŝ as surrogate of s, and use the noisy structure to derive an approximation
of F and P .

3

Numerical results

We first demonstrate that σ̂ obtained through the method described in section 2.1 is a good estimate of σ. For simulation purposes, we use a segment
of ubiquitin with seven peptide planes (residue 1-8) containing 21 N − H,
C − CA and C − N bonds. We note that in all the simulations, we do not
consider the RDC of the nonexistent N − H bond for proline. In Fig. 2(left),
we plot σ̂ v.s. σ. The simulation shows a close agreement between σ̂ and σ,
especially when the angular noise is less than 12 degrees.
We next show that the SIMEX-like method proposed in section 2 is able
to reduce the bias in eigenvalue estimation, where the bias of an estimator
θ̂ of parameter θ is defined to be
Bias(θ̂) = hθ̂i − θ.
h·i denotes averaging over the distribution of data. For this simulation, we
use a specific ground truth Saupe tensor and the aforementioned ubiquitin
fragment to generate precise RDC measurements. From the fragment, 200
realizations of noisy conformation are generated with σ = 20◦ . To obtain
Λ̃, we set n1 = 8000 when simulating Ab in step (2) of the Monte Carlo
procedure. In Fig. 2(right), we see that the values of hΛ̃iαi ,βi (Red dotted
line) obtained from averaging over 200 samples are almost the same as the
eigenvalues of S (Black line), while there is a clear bias in the estimator
Λ(Ŝ) (Blue dotted line).
9

1
0.9

Normalized Eigenvalues

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

5

10

15

σ

20

25

30

Figure 1 Plot of the eigenvalues of the OLS estimator Ŝ normalized by
the eigenvalues of S v.s. σ. Increasing the noise level biases the eigenvalues
towards zero. A fragment of ubiquitin composed of 7 peptide planes (residue
1-8) and a specific Saupe tensor S is used for the simulation and each point
in the plot is computed from 200 different realizations of αi , βi ’s. ]

10

40
35

σ estimate

30
25
20
15
10
5
0
0

5

10

15

−0.9 −0.7 −0.5 −0.3

−0.6 −0.4 −0.2
−3

0

0.6

70

70

60

60

60

50

50

50

40

40

40

30

30

30

20

20

20

10

10

10

Λ (S)
xx

0

1

Λ (S)
yy

0

1.2
−3

× 10

70

0

0.8

25

−3

× 10

Frequency

20

σ (Degree)

× 10

Λ (S)
zz

Figure 2 Above: Plot of σ̂ v.s. σ. For a given noise level σ, σ̂ is averaged
over 200 different realizations of αi , βi ’s. Below: Histograms of the diagonal
11
entries of Λ(Ŝ) and Λ̃ obtained from 200 fragment conformations with 20◦
noise on the torsion angles. The values of Λ(S), hΛ(Ŝ)iαi ,βi and hΛ̃iαi ,βi are
denoted by black, blue and red line respectively.

3.1

Estimation of Saupe tensor eigenvalues from multiple
molecular fragments

While the proposed eigenvalue estimator Λ̃ has less bias, this does not obligate that Λ̃ has a lower mean squared error (MSE).This can be understood
from the bias-variance decomposition, which is a classical way in statistics
to decompose the MSE of an estimator θ̂ [28]. The MSE of an estimator θ̂
admits the following decomposition
MSE(θ̂) =
=
=
=

h(θ − θ̂)2 i
h(θ − hθ̂i + hθ̂i − θ̂)2 i
Bias(θ̂)2 + Var(θ̂) + 2(θ − hθ̂i)hhθ̂i − θ̂i
Bias(θ̂)2 + Var(θ̂)

(20)

Var(θ̂) denotes the variance of θ̂. Although we achieve less bias with the
estimator Λ̃, we pay the price of having larger variance due to bias estimation involved in obtaining Λ̃. This increase in variance can lead to Λ̃
having higher MSE than Λ(Ŝ). From this point of view, when estimating
the Saupe tensor eigenvalue using a single template fragment, the Monte
Carlo method for debiasing may seem unnecessary or even disadvantageous.
However, when multiple template fragments are available, the average of Λ̃
over these fragments, denoted Λ̃ave , enjoys variance reduction proportional
to the number of fragments. Therefore in the case when there are many
fragments, it is worth paying the price of increased variance because the systematic bias error cannot be reduced via averaging. In the rest of the section,
we use Λave (Ŝ) to denote the average of Λ(Ŝ) over multiple fragments.
We now demonstrate the usefulness of our method under the setting of
Molecular Fragment Replacement (MFR) approach [7, 12]. When RDCs
are measured in two different alignment medias for a protein of unknown
structure, the MFR method can construct its structure by combining short
homologous fragments obtained from chemical shift and dipolar homology
database mining. Typically for every protein fragment of seven residues,
10 homologous structures are searched based on the similarity of chemical
shifts and the goodness of Saupe tensor fit to the observed RDC. When
OLS is used to fit the Saupe tensor with design matrix A constructed from
homologous structures, one can average all OLS eigenvalue estimated to
obtain improved estimators of the parameters such as alignment magnitude
and rhombicity that depend on the eigenvalues [12]. These parameters can
in turn be used in a simulated annealing procedure such as XPLOR-NIH
[22] to refine the structure.
We first use synthetic data to demonstrate our method. We generate
12

12 random Saupe tensors, by sampling two eigenvalues from the uniform
distribution on [−10−3 , 0] and [0, 10−3 ] respectively, and extract the third
eigenvalue by requiring Λ(S)xx + Λ(S)yy + Λ(S)zz = 0. The orthogonal
matrix U (S) is sampled uniformly from the group of 3 × 3 orthogonal matrices, by computing the orthogonal factor in the polar decomposition of a
3 × 3 Gaussian random matrix [2]. After obtaining the RDC dnm ’s from the
clean structure and the ground truth Saupe tensor, under each simulated
alignment condition we add structural noise of magnitude σ to every fragment of seven peptide planes of the ubiquitin structure obtained from X-ray
crystallography (PDB ID 1UBQ). We only consider the first 71 residues of
the 76 residues of ubiquitin, as there are few RDC reported for the last five
residues of ubiquitin. This gives a total of 64 fragments. We evaluate the
estimators of the Saupe tensor eigenvalues Λave (Ŝ) and Λ̃ave computed from
the average of Λ(Ŝ) and Λ̃ of all fragments, by comparing their fractional
errors averaged over the 12 different Saupe tensors and torsion angle noise
realizations in Fig. 3. The fractional error is defined as
kΛave (Ŝ) − Λ(S)kF
kΛ(S)kF

and

kΛ̃ave − Λ(S)kF
.
kΛ(S)kF

In this simulation, the fractional error of Λave (Ŝ) is at least three times
larger than Λ̃ave .
We finally apply this method to estimate the Saupe tensor of ubiquitin in two different alignment medias using the experimental RDC data in
[6]. From 600 homologous structures returned by MFR homology search,
each containing seven residues, we obtain 600 Saupe tensor estimates using
OLS. Since we expect our method to have a significant effect for fragments
severely corrupted by structural noise, we average the fragments with residual RMS above a certain threshold and plot Λave (Ŝ) and Λ̃ave normalized by
Λ(S) v.s. RMS thresholds. To get an estimated (and approximate) ground
truth Saupe tensor S, we use the high resolution ubiquitin structure 1UBQ
obtained from X-ray crystallography [27] to fit the RDC data. We demonstrate the results in Fig. 4. Other than the estimators for Λ(S)yy of the
second alignment media which has a large percent error due to the relatively
small magnitude of Λ(S)yy , Λ̃ave typically achieves 0.9 of the ground truth
value, whereas Λave (Ŝ) can shrink to 0.8 of the value of Λ(S) when only
the fragments of high RMS are used in averaging. We therefore recommend
the use of our proposed bias removing method when estimating eigenvalues
from multiple noisy fragments.

13

0.14

0.12

SVD (Debiased)
SVD

Fractional error

0.1

0.08

0.06

0.04

0.02

0
0

5

10

σ (Degree)

15

20

Figure 3 Plot of the fractional error of Λave (Ŝ) and Λ̃ave v.s. σ. Each data
point is averaged over 12 different Saupe tensor and noise realizations for
1UBQ. The plot shows a clear advantage of the bias reduced estimator over
the OLS estimator.

14

Normalized Eigenvalue

Λ (S) = −5.34e−04
xx

1
0.95
0.9
0.85
0.8
0.75
0.7
0

0.5

1

Λ (S) = −3.18e−04

1
0.95
0.9
0.85
0.8
0.75
0.7
1.5
0

0.5

Normalized Eigenvalue

RMS x 10−4
Λxx(S) = −1.24e−03
1
0.95
0.9
0.85
0.8
0.75
0.7
0

1

2

Λ (S)=8.53e−04

yy

1

zz

1
0.95
0.9
0.85
0.8
0.75
0.7
1.5
0

RMS x 10−4
Λyy(S) = −0.19e−03

1
0.9
0.8
0.7
0.6
0.5
0.4
3
0

1

RMS x 10−4

2

0.5

1

1.5

RMS x 10−4
Λzz(S)=1.44e−03

1
0.95
0.9
0.85
0.8
0.75
0.7
3
0

RMS x 10−4

1

2

3

RMS x 10−4

Figure 4 Plot of the eigenvalue estimators normalized by Λ(S) v.s. residual
RMS thresholds. Estimators are obtained from experimental RDC measurements in two different alignment medias. While the magnitude of Λ̃ave (Red
curves) and Λave (Ŝ) (Blue dotted curves) both decrease as low quality (high
RMS) fragments are solely used in averaging, Λ̃ave in general is within 90%
of the ground truth value but Λ(S) drops to 80% of Λave (Ŝ). The value of
Λ(S) for both alignment medias are indicated in the plot title.

15

4

Effect of additive noise on Saupe tensor estimation

So far we have been neglecting the presence of additive noise on dnm , which
is considered by [15]. We define the noisy RDC measurements corrupted by
additive noise as
dadd = d + σadd ε
(21)
where entries of the column vector ε are i.i.d random variables with mean
zero. In this section, we show using perturbation theory that this type of
additive noise biases the eigenvalue magnitude positively, therefore it cannot
explain the magnitude shrinkage we see when fitting the Saupe tensor to real
RDC data (Fig. 4). Moreover, the order of magnitude of this positive bias
is not sufficient to explain the error between Ŝ and S. This has been noted
by the authors of [15] that in order to account for the size of the OLS misfit,
an uncertainly of 2-3 Hz for the RDC measurements is required although
the experimental uncertainty is only about 0.2-0.5 Hz. This is the reason
why in this paper we focus on removing the bias that arises from structural
noise.
Let
S = U (S)Λ(S)U (S)T
(22)
be the eigendecomposition of S. Assuming the eigenvalues of S are nondegenerate, the second order perturbation theory [13] states
λj (Ŝ) ≈ λj (S) + U (S)Tj (Ŝ − S)U (S)j
X (U (S)T (Ŝ − S)U (S)j )2
k
.
+
λj (S) − λk (S)

(23)

k=x,y,z,
k6=j

Averaging the perturbation expansion over the distribution of ε, we get
λj (Ŝ) ≈ λj (S) +

X

k=x,y,z,
k6=j

h(U (S)Tk (Ŝ − S)U (S)j )2 iε
.
λj (S) − λk (S)

(24)

Here we use the fact that hŜ − Siε = 0 since
hŝ − siε = h(AT A)−1 AT (As + ε) − siε = h(AT A)−1 AT εiε = 0
The expression in (24) reveals that in the presence of noise, the largest
eigenvalue of Ŝ is always greater than the largest eigenvalue of S, while the
16

smallest eigenvalue behaves in the exact opposite manner. Such effect of bias
of pushing the extreme eigenvalues outwards is also commonly seen in the
context of estimating the extreme eigenvalues of covariance matrices [21].
For this type of bias we now give an estimate of its order of magnitude.
First we bound the numerator in the second order correction term in (24):
(U (S)Tk (Ŝ − S)U (S)j )2 = Tr((Ŝ − S)U (S)j U (S)Tk )2
≤ kŜ − Sk2F kU (S)j U (S)Tk k2F
≤ 3kŝ − sk22 .

(25)

The first inequality results from Cauchy-Schwarz inequality, and the second
inequality relies on the fact that kU (S)j U (S)Tk kF = 1 and kŜ − Sk2F ≤
3kŝ − sk22 , which can be verified easily. It is a classical result [11] that the
OLS estimator has covariance matrix
2
(AT A)−1 ,
h(ŝ − s)(ŝ − s)T iε = σadd

(26)

2
Tr((AT A)−1 ).
hkŝ − sk22 iε = σadd

(27)

therefore
Using (27), (25) we obtain an upperbound for the bias in (24). Taking
λz (Ŝ) for example:
λz (Ŝ) − λz (S) .

3Tr((AT A)−1 ) 2
σ
λz (S) − λy (S) add

(28)

We now give an estimate of the order of magnitude of the bias. Since the
magnitude of the extreme eigenvalues of the Saupe tensor is around 10−3 ,
for example for the two RDC datasets acquired for ubiquitin, we simply
assume λz (S) − λy (S) ∼ 10−4 . The typical experimental uncertainty for
RDC measurements is about 0.2 Hz - 0.5 Hz, and the dipolar coupling
max for e.g. N − H bonds, is about 23 kHz, therefore the noise
constant Dnm
magnitude σadd of the additive noise on the normalized dipolar coupling is
about 0.5/(23 × 103 ) ≈ two × 10−5 . For ubiquitin, the average value of
Tr((AT A)−1 ) for fragments containing seven peptide planes is about 1.35.
Using these numbers in (28), we get
λz (Ŝ) − λz (S) . 1.6 × 10−5 ,
which amounts to 1-2% error when λz (S) ∼ 10−3 . This cannot explain
the 10% or larger error in fitting Saupe tensor to real RDC datasets using
homology fragments in the previous section.
17

We present a simulation to illustrate the bias in OLS eigenvalues estimation in the presence of additive noise. We use the Saupe tensor eigenvalues
for ubiquitin in the first alignment media presented in Fig. 4, and a ubiquitin
fragment consisting of seven peptide planes (residue 1-8) for this simulation.
We generate noisy datasets using the noise model
dadd = As + σadd ε.
For every noise level, we average Λ(Ŝ) normalized by Λ(S) over 500 different
realizations of s and ε where entries of ε are i.i.d. random normal variables.
The different realization of s are generated from S = U (S)Λ(S)U (S)T where
Λ(S) is fixed but U (S) is sampled uniformly from the orthogonal group in
R3 . We vary the orientation of the Saupe tensor since it is clear from (24)
that the bias depends on U (S). We change σadd from 0 to 10% of λz (S)
and present the results in Fig. 5. We note again from previous calculations,
σadd ∼ 2 × 10−5 , which amounts to 2-3% of the λz (S) = 0.85 × 10−3 considered. As shown in the simulation and our crude estimate, such magnitude
of noise gives rise to bias error of about 1%. Even in the case of having very
noisy RDC (having noise magnitude 10% of λz (S)), the bias error caused
by additive noise is around 3%. Whereas in a typical MFR search with
torsion angle tolerance being set to ±20◦ − 30◦ [30], the simulation in Fig. 1
suggests structural noise can cause bias error sometimes much greater than
10%. Therefore in this paper we focus on removing the bias that arises from
structural noise. In the case when accurate template structure is available
and the additive noise is a concern, we refer readers to the appendix for the
removal of such bias using an analytic expression derived from perturbation
theory.

5

Conclusion

We observe a negative bias when estimating the Saupe tensor eigenvalues
through the classical SVD method, in the presence of structural noise on the
template structure due to torsion angle noise. We present a Monte Carlo
method that simulates noise on the template structure by perturbing the
torsion angles and use the simulated structure to estimate the bias in the
eigenvalues. We demonstrate the effectiveness of our method in reducing the
error arising from bias when estimating Saupe tensor eigenvalues from multiple protein fragments, which is a natural setting to consider when building
protein structure from homologous substructures.

18

1.035

Normalized Eigenvalues

1.03
1.025
1.02
1.015
1.01
1.005
1
0

0.02

0.04
σ

0.06
/λ (S)

add

0.08

0.1

z

Figure 5 Plot of the three eigenvalues of the OLS estimator Ŝ normalized
by the eigenvalues of S v.s. σadd under noise model (21). Each point is
averaged over 2000 noise and Saupe tensor realizations. Increasing the noise
level biases the eigenvalues positively, unlike the case for structural noise.
At 10% noise level, the bias is about 3%.

19

6

Acknowledgement

The research of AS was partially supported by award R01GM090200 from
the NIGMS, by awards FA9550-12-1-0317 and FA9550-13-1-0076 from AFOSR,
by award LTR DTD 06-05-2012 from the Simons Foundation, and the Moore
Foundation Data Driven Discovery Investigator award.

7
7.1

Appendix
Removing bias from additive noise

Define a linear operator L : R5 → R3×3 that forms a Saupe tensor S from
the vector s as


−s(1) − s(2) s(3) s(4)
L(s) = 
s(3)
s(1) s(5) , s ∈ R5 .
(29)
s(4)
s(5) s(2)
For the additive noise model (21) we have
Ŝ = L(ŝ) = L((AT A)−1 AT dadd ) = S + L((AT A)−1 AT ε).

(30)

We also define the adjoint operator of L, L∗ : R3×3 → R5 through the
relation
Tr(X T L(y)) = L∗ (X)T y,
(31)
for every y ∈ R5 and X ∈ R3×3 . To obtain the form of L∗ , we let y ∈
{e1 , . . . , e5 }, where ei (i) = 1 and ei (j) = 0 if j 6= i. Plugging such y into Eq.
(31), we get
L∗ (X)1
L∗ (X)2
L∗ (X)3
L∗ (X)4
L∗ (X)5

=
=
=
=
=

−Xxx + Xyy
−Xxx + Xzz
Xxy + Xyx
Xxz + Xzx
Xyz + Xzy

20

(32)

Using such notion of the adjoint operator, the perturbation series in (24)
can be written as
hλj (Ŝ)iε ≈ λj (S) +

X

k=x,y,z,
k6=j

= λj (S) + Tr

 X

h((ŝ − s)T L∗ (U (S)k U (S)Tj ))2 iε
λj (S) − λk (S)

L∗ (U (S)k U (S)Tj )L∗ (U (S)k U (S)Tj )T
λj (S) − λk (S)

k=x,y,z,
k6=j




Var(ŝ)
(33)

where
2
Var(ŝ) ≡ h(ŝ − s)(ŝ − s)T iε = (AT A)−1 σadd
.

(34)

Therefore we can subtract the second order term in (33) to correct for the
bias in the eigenvalues. Although we do not know the eigenvectors and
eigenvalues of S, we can replace them with the eigenvectors and eigenvalues
Ŝ. This change will only affect on the higher order terms in the perturbation
series.

References
[1] Martin Blackledge, Recent progress in the study of biomolecular structure and dynamics in solution from residual dipolar couplings, Progress
in Nuclear Magnetic Resonance Spectroscopy 46 (2005), no. 1, 23–61.
[2] Gordon Blower, Random matrices: high dimensional phenomena, vol.
367, Cambridge University Press, 2009.
[3] Raymond J Carroll, David Ruppert, Leonard A Stefanski, and
Ciprian M Crainiceanu, Measurement error in nonlinear models: a modern perspective, CRC press, 2006.
[4] G Clore, A Gronenborn, and A Bax, A robust method for determining the magnitude of the fully asymmetric alignment tensor of oriented
macromolecules in the absence of structural information, Journal of
magnetic resonance 133 (1998), no. 1, 216–21.
[5] John R Cook and Leonard A Stefanski, Simulation-extrapolation estimation in parametric measurement error models, Journal of the American Statistical Association 89 (1994), no. 428, 1314–1328.
21

[6] Gabriel Cornilescu, John L Marquardt, Marcel Ottiger, and Ad Bax,
Validation of protein structure from anisotropic carbonyl chemical shifts
in a dilute liquid crystalline phase, Journal of the American Chemical
Society 120 (1998), no. 27, 6836–6837.
[7] Frank Delaglio, Georg Kontaxis, and Ad Bax, Protein structure determination using molecular fragment replacement and NMR dipolar
couplings, Journal of the American Chemical Society 122 (2000), no. 9,
2142–2143.
[8] Bradley Efron and Robert J Tibshirani, An introduction to the bootstrap, CRC press, 1994.
[9] Theo Gasser, HG Muller, and Volker Mammitzsch, Kernels for nonparametric curve estimation, Journal of the Royal Statistical Society.
Series B (Methodological) (1985), 238–252.
[10] Michael Grant and Stephen Boyd, CVX: Matlab software for disciplined
convex programming, version 2.1, http://cvxr.com/cvx, March 2014.
[11] Jürgen Groß, Linear regression, vol. 175, Springer Science & Business
Media, 2003.
[12] Georg Kontaxis, Frank Delaglio, and Ad Bax, Molecular fragment replacement approach to protein structure determination by chemical shift
and dipolar homology database mining, Methods in Enzymology 394
(2005), 42–78.
[13] Lev Davidovich Landau and Evgenii Mikhailovich Lifshitz, Quantum
mechanics: non-relativistic theory, vol. 3, Elsevier, 2013.
[14] Rebecca S Lipsitz and Nico Tjandra, Residual dipolar couplings in NMR
structure analysis, Annu. Rev. Biophys. Biomol. Struct. 33 (2004), 387–
413.
[15] Judit A Losonczi, Michael Andrec, Mark WF Fischer, and James H
Prestegard, Order matrix analysis of residual dipolar couplings using singular value decomposition, Journal of Magnetic Resonance 138
(1999), no. 2, 334–342.
[16] Jens Meiler, Jeanine J Prompers, Wolfgang Peti, Christian Griesinger,
and Rafael Brüschweiler, Model-free approach to the dynamic interpretation of residual dipolar couplings in globular proteins, Journal of the
American Chemical Society 123 (2001), no. 25, 6098–6107.
22

[17] Eva Meirovitch, Donghan Lee, Korvin FA Walter, and Christian
Griesinger, Standard tensorial analysis of local ordering in proteins from
residual dipolar couplings, The Journal of Physical Chemistry B 116
(2012), no. 21, 6106–6117.
[18] T Michael Sabo, Colin A Smith, David Ban, Adam Mazur, Donghan Lee, and Christian Griesinger, Orium: Optimized rdc-based iterative and unified model-free analysis, Journal of biomolecular NMR 58
(2014), no. 4, 287–301.
[19] L. Salmon and M. Blackledge, Investigating protein conformational energy landscapes and atomic resolution dynamics from NMR dipolar couplings: a review, Rep Prog Phys 78 (2015), no. 12, 126601.
[20] AZ. Saupe, Kernresonanzen in kristallinen flssigkeiten und in
kristallinflssigen lsungen. teil i, Naturforsch. 19 (1964), no. a, 161–171.
[21] Juliane Schäfer and Korbinian Strimmer, A shrinkage approach to largescale covariance matrix estimation and implications for functional genomics, Statistical applications in genetics and molecular biology 4
(2005), no. 1.
[22] Charles D Schwieters, John J Kuszewski, Nico Tjandra, and G Marius
Clore, The Xplor-NIH NMR molecular structure determination package,
Journal of Magnetic Resonance 160 (2003), no. 1, 65–73.
[23] LA Stefanski and JR Cook, Simulation-extrapolation: the measurement
error jackknife, Journal of the American Statistical Association 90
(1995), no. 432, 1247–1256.
[24] Joel R Tolman, A novel approach to the retrieval of structural and
dynamic information from residual dipolar couplings using several oriented media in biomolecular nmr spectroscopy, Journal of the American
Chemical Society 124 (2002), no. 40, 12020–12030.
[25] Joel R Tolman and Ke Ruan, NMR residual dipolar couplings as probes
of biomolecular dynamics, Chemical Reviews 106 (2006), no. 5, 1720–
1736.
[26] John W Tukey, Exploratory data analysis, Addison-Wesley Series in
Behavioral Science: Quantitative Methods, Reading, Mass.: AddisonWesley, 1977, vol. 1, 1977.

23

[27] Senadhi Vijay-Kumar, Charles E Bugg, and William J Cook, Structure
of ubiquitin refined at 1.8Å resolution, Journal of Molecular Biology
194 (1987), no. 3, 531–544.
[28] Larry Wasserman, All of statistics: a concise course in statistical inference, Springer Science & Business Media, 2013.
[29] Lukas N Wirz and Jane R Allison, Fitting alignment tensor components
to experimental RDCs, CSAs and RQCs, Journal of Biomolecular NMR
(2015), 1–5.
[30] Zhengrong Wu, Frank Delaglio, Keith Wyatt, Graeme Wistow, and
Ad Bax, Solution structure of γs-crystallin by molecular fragment replacement NMR, Protein science 14 (2005), no. 12, 3101–3114.
[31] Markus Zweckstetter and Ad Bax, Evaluation of uncertainty in alignment tensors obtained from dipolar couplings, Journal of Biomolecular
NMR 23 (2002), no. 2, 127–137.

24

