High-order analytic translation matrix elements for real-space six-dimensional polar Fourier correlations

Ritchie, D.W.

doi:10.1107/S002188980502474X

research papers

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 38| Part 5| October 2005| Pages 808-818

doi:10.1107/S002188980502474X

High-order analytic translation matrix elements for real-space six-dimensional polar Fourier correlations

David W. Ritchie ^a ^*

^aDepartment of Computing Science, University of Aberdeen, Aberdeen AB24 3UE, Scotland
^*Correspondence e-mail: dritchie@csd.abdn.ac.uk

(Received 24 May 2005; accepted 2 August 2005)

Analytic expressions are presented for calculating translations of high-order three-dimensional expansions of orthonormal real spherical harmonic and Gaussian-type or exponential-type radial basis functions. When used with real spherical harmonic rotation matrices, the resulting translation matrices provide a fully analytic method of calculating six-dimensional real-space rotational–translational correlations. The correlation algorithm is demonstrated by using an exhaustive search to superpose the steric density functions of a pair of similar globular proteins in a matter of seconds on a contemporary personal computer. It is proposed that the techniques described could be used to accelerate the calculation of e.g. real-space electron density correlations in molecular replacement, docking proteins into electron microscopy density maps, and searching the Protein Data Bank for structural homologues.

Keywords: polar Fourier correlation; 6D correlation; fast Fourier transform; fast Hartley transform; spherical Bessel transform; Gaussian-type orbitals; Slater-type orbitals; exponential-type orbitals; Laguerre polynomials; protein shape; protein docking.

1. Introduction

Solving protein crystal structures by molecular replacement (MR) (Rossmann, 1990 , 2001 ), fitting protein structures into electron microscopy (EM) density maps (Roseman, 2000 ; Rossmann, 2000 ; Frank, 2002 ), and predicting the docking mode of a pair of proteins (Katchalski-Katzir et al., 1992 ; Ritchie & Kemp, 2000 ) are three examples of tasks which typically involve searching for the coordinate transformation that maximizes a six-dimensional correlation. Often, such correlations can be expressed as the overlap between one or more pairs of three-dimensional scalar functions, such as electron densities, and it is primarily this type of correlation that is considered here. For example, if some phase information is available from multiple anomalous dispersion (MAD) or multipe isomorphous replacement (MIR) data, or from a partial model, the phased translation function may be used to correlate a low-resolution electron density map with the density from a model structure to help locate the model in the crystal unit cell (Read & Schierbeek, 1988 ). This type of translational correlation may be accelerated using a three-dimensional fast Fourier transform (FFT). However, the correct solution can sometimes be missed if the rotational solution is in error. This can often be corrected by varying the rotational orientations (Fujinaga & Read, 1987 ) or by translating multiple rotational peaks (Navaza, 2001 ), but both these techniques entail additional computational cost. Hence a more thorough but economical method of correlating electron densities could be beneficial. We note that combined six-dimensional rotational-translational Patterson correlation (PC) approaches have recently been successful in solving difficult MR cases without using phase information (Chang & Lewis, 1997 ; Sheriff et al., 1999 ). However, these approaches are based on point-wise evaluations of the PC in reciprocal space (Chang & Lewis, 1997; Kissinger et al., 1999 ; Jamrog & Zhang, 2003 ) and hence appear not to be amenable to the type of correlation discussed here. Correlating electron densities may also be used to locate secondary structure fragments in density maps (Cowtan, 1998 ), and has been proposed as a novel way to search the Protein Data Bank (PDB) (Sussman et al., 1998 ) for candidate molecular replacement models (Cowtan, 2001 ). Similarly, protein–protein docking algorithms often use a three-dimensional translational FFT to search for good docking orientations (Katchalski-Katzir et al., 1992; Vakser, 1995 ). However, in each of these examples the translational FFTs must be repeated for a large number of rotational orientations in order to cover the remaining three degrees of freedom. Hence, more efficient correlation techniques would be useful in these areas as well.

The FFT may also be used to accelerate rotational searches. For example, the MR fast rotation function represents a Patterson map of self-vectors as an expansion of spherical Bessel and spherical harmonic functions and uses a two-dimensional FFT to orient a search model with respect to the observed diffraction data (Navaza, 1987 ). Vagin & Isupov (2001 ) adapted the fast rotation function to implement a phased rotation function (rotational density correlation) which allows a spherically averaged search model to be located in the unit cell using the phased translation function before solving its rotational orientation. However, this approach still requires rotational and translational correlations to be treated separately. More recently, Kovacs et al. (2003 ) employed an elegant Euler angle factorization in order to fit proteins into EM density maps using a five-dimensional rotational FFT. However, their approach entails calculating certain residual overlap integrals which must be computed numerically.

In our own areas of interest, namely protein docking and shape comparison, we developed a six-dimensional spherical polar Fourier (SPF) correlation approach (Ritchie & Kemp, 2000) to try to address the limitations of the Cartesian FFT docking methods. In the SPF representation, each three-dimensional scalar property of interest is written as an expansion of real spherical harmonic plus orthonormal radial basis functions. This allows the correlation between a pair of proteins to be expressed naturally as a search over one translational and five rotational degrees of freedom. Because the rotational properties of the spherical harmonics have been described extensively elsewhere (see e.g. Rose, 1957 ; Biedenharn & Louck, 1981a ), this article focuses on calculating translations for SPF expansions. Formerly, we calculated translations by numerical integration in the (r, θ) plane (Ritchie & Kemp, 2000). Here, we give analytic expressions for efficiently calculating high-order translation matrix elements for both harmonic oscillator (HO) and orthogonal Coulomb-type radial basis functions using spherical Bessel transforms. The HO functions, sometimes referred to as Gaussian-type orbitals (GTOs), have long been used in molecular orbital calculations because of the ease of calculating overlap integrals between pairs of such functions (Boys, 1950 ). The Coulomb-type functions, sometimes referred to as exponential-type orbitals (ETOs) (Guseinov, 2001 , 2002 ), have an exponential factor in place of the Gaussian and correspond to orthogonalized Slater-type atomic orbitals (STOs) (Barnett & Coulson, 1951 ) and B-functions (Filter & Steinborn, 1978 ), or equivalently Bessel-type orbitals (BTOs) (Ozdogan et al., 2005 ).

In what follows, the abbreviations GTO and ETO will be used to refer to the radial functions used here for consistency with the recent quantum mechanics literature. Despite much prior work, there is still considerable interest in developing efficient and stable methods of calculating integrals over GTOs (Arakane & Matsuoka, 1999 ; Chiu & Moharerrzadeh, 1999 ), and especially STOs (Filter & Steinborn, 1980 ; Weniger & Steinborn, 1983 ; Hierse & Oppeneer, 1994 ; Ozdogan et al., 2005; Rico et al., 2005 ). Nonetheless, to our knowledge, the analytic expressions presented here to translate ETO basis functions and to relate ETOs to BTOs are novel. We also describe non-orthogonal expansions of GTOs and ETOs which give at least a twofold speed-up compared with calculating translations in the orthogonal bases.

Our expressions have been tested for numerical accuracy up to order N = 32, which is sufficient for high-resolution protein–protein docking correlations (Ritchie, 2003 ). We find it is necessary to use an extended-precision arithmetic library to calculate high-order translation matrix elements analytically, although all subsequent calculations may be performed in ordinary double-precision arithmetic. However, because much of our docking algorithm has been described previously (Ritchie & Kemp, 2000), the calculations are illustrated here by superposing a pair of similar superantigen proteins. In contrast to multi-dimensional FFT approaches, our correlation is calculated as a combinatorial search over pairs of rotated and translated coefficient vectors. With low-order polynomial expansions there is no advantage in using an FFT, although for high-order correlations our algorithm can be accelerated using a one-dimensional real fast Hartley transform (FHT) (Bracewell, 1999 ). In the example given, a good superposition is achieved using low-order expansions to N = 6, with a total calculation time of under 8 s for an exhaustive six-dimensional search on a 2 GHz Pentium Xeon processor. This corresponds to evaluating approximately 6 × 10⁶ trial orientations s⁻¹. This demonstrates that the SPF correlation technique offers a novel and powerful approach for problems that involve calculating real-space six-dimensional rotational–translational correlations.

2. Methods

2.1. Polar Fourier expansions

In the SPF approach, each scalar property of interest, A(r), is represented as an infinite expansion, truncated to order N as

$[A({\bf r}) = \sum_{n = 1}^{N} \sum_{l = 0}^{n-1} \sum_{m = -l}^{l} a_{nlm} R_{nl}(r) y_{lm}(\theta,\phi), \eqno (1)]$

where a_nlm are the expansion coefficients, $[y_{lm}(\theta,\phi)]$ = $[\vartheta_{l|m|} (\theta) \varphi_m(\phi)]$ are real normalized spherical harmonics, and R_nl(r) are orthonormal GTO or ETO radial basis functions. Generally, we use GTOs to represent steric shape and the more diffuse ETOs to represent electrostatic properties. We follow the quantum chemistry convention in which the radial index n, or principal quantum number, counts from unity. Hence the highest harmonic order and highest polynomial power in any individual coordinate is L = N − 1. An expansion to order N involves N(N + 1)(2N + 1)/6 coefficients. These are calculated just once for each property as described previously (Ritchie & Kemp, 2000).

2.2. Overlap Integrals as translation matrix elements

First, we consider the correlation between a fixed `body' A, or scalar function A(r), and a moving body B, or function B(r), under an active translation of B by T = (R, 0, 0) along the positive z axis:

$[C({\bf T}) = \int A({\bf r}) B({\bf T}^{-1} {\bf r}) \, {\rm d}V. \eqno (2)]$

Substituting the expansions of A(r) and B(r) gives

$[\eqalignno{ C({\bf T}) =\hskip.2em& \sum_{nlm} \sum_{n'l'm'} a_{nlm} b_{n'l'm'} \cr&\!\times \int R_{nl}(r) y_{lm}(\theta,\phi) R_{n'l'}(r') y_{l'm'}(\theta',\phi') \, {\rm d}V ,&(3)}]$

where the shorthand notation $[\textstyle\sum_{nlm}]$ etc. is used to indicate summation over the subscript ranges given in equation (1), and where r = $[(r,\theta,\phi)]$ and r′ = r − T = $[ (r',\theta',\phi')]$ . In this case, $[\phi]$ and $[\phi']$ remain coincident, so the circular functions, $[\varphi_m(\phi)]$ , may be integrated out and equation (3) reduces to a sum over two-dimensional integrals in the $[(r,\theta)]$ plane. Because these integrals clearly depend only on the distance R, we write

$[C(R) = \sum_{nlm} \sum_{n'l'm'} a_{nlm} b_{n'l'm'} T^{(|m|)}_{nl,n'l'}(R) \delta_{mm'} \eqno (4)]$

and interpret each $[T^{(|m|)}_{nl,n'l'}(R)]$ as a matrix element of the translation operator. For example, from equation (4) it can be seen that the two sums

$[b^{R}_{nlm} = \sum_{n'l'} T^{(|m|)}_{nl,n'l'}(R) b_{n'l'm }\eqno (5)]$

and (after re-labeling the subscripts)

$[a^R_{nlm} = \sum_{n'l'} T^{(|m|)}_{n'l',nl}(R) a_{n'l'm} \eqno (6)]$

represent a positive translation of the body B, or equivalently a negative translation of the body A, respectively. The translation matrices are obviously five-dimensional quantities. However, because they do not depend on the sign of m (see below), it is useful to consider each matrix as being composed of $[\sum_{m = 0}^{N-1} ]$ = N two-dimensional arrays, each indexed by $[\sum_{nl}]$ = N(N + 1)/2 possible values for each pair of nl subscripts. The matrix elements vanish trivially where |m| > l. The notation used here is intended to be consistent with the usual convention for the complex and real spherical harmonic rotation matrix elements, $[D^{(l)}_{m'm}(\alpha,\beta,\gamma)]$ and $[R^{(l)}_{m'm}(\alpha,\beta,\gamma)]$ , respectively, in the sense that a positive z translation of the basis functions is expressed as

$[R_{nl}(r')y_{lm}(\theta',\phi) = \sum_{n' = 1}^{\infty} \sum_{l' = 0}^{n'-1} T^{(|m|)}_{n'l',nl}(R) R_{n'l'}(r)y_{l'm}(\theta,\phi). \eqno (7)]$

2.3. Spherical Bessel transform method of calculation

If the spherical Bessel transform of R_nl(r) is defined as

$[\tilde{R}_{nl}(\beta) = (2/\pi)^{1/2} \int\limits_0^{\infty} R_{nl}(r) j_{l}(\beta r) r^2 \,{\rm d}r, \eqno (8)]$

where j_l(z) is the spherical Bessel function, then the inverse transform is given by (Hochstadt, 1971 ):

$[{R}_{nl}(r) = (2/\pi)^{1/2} \int\limits_0^{\infty} \tilde{R}_{nl}(\beta) j_{l}(\beta r) \beta^2 \,{\rm d}\beta. \eqno (9)]$

It is shown in Appendix A that the translation matrix elements for SPF basis functions may be calculated as a sum over one-dimensional inverse Bessel transforms

$[T_{n'l',nl}^{(|m|)}(R) = \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \int\limits_0^{\infty} \tilde{R}_{n'l'}(\beta) \tilde{R}_{nl}(\beta) j_{k}(\beta R) \beta^2 \,{\rm d}\beta ,\eqno (10)]$

where the coefficients $[A_k^{(ll'|m|)}]$ are given by

$[\eqalignno{A_k^{(ll'|m|)} =\hskip.2em & (-1)^{ (k+l'-l)/2 + m} (2k+1) \big [{(2l+1)(2l'+1)} \big]^{1/2} \cr & \!\times \bigg({l \, l' \, k \atop 0 \, 0 \, 0 } \bigg) \bigg({l \ l' \ k \atop m \, \overline{m} \ 0 } \bigg). &(11)}]$

From the permutational symmetries of the second 3–j symbol, the right-hand side is independent of the sign of m, hence justifying the use of |m| to label the matrix elements. Because the first 3–j symbol vanishes whenever l + l′ + k is odd, it can be seen that the non-vanishing coefficients are always real and that the summation in equation (10) need only be calculated for even increments: k = |l − l′|, |l − l′| + 2, …, l + l′. By similar arguments (see Appendix A), it can also be shown that

$[T_{nl,n'l'}^{(|m|)}(R) = T_{n'l',nl}^{(|m|)}(-R) = (-1)^{l'-l} T_{n'l',nl}^{(|m|)}(R). \eqno (12)]$

Consequently, nearly half of all matrix elements can be found by symmetry. Given that the original basis functions form a complete orthonormal set, it is straightforward to show that the translation matrices are also orthonormal in the sense that

$[\sum_{n' = 1}^{\infty} \sum_{l' = 0}^{n'-1} T_{n'l',nl}^{(|m|)}(R) T_{n'l',n''l''}^{(|m|)}(R) = \delta_{n n''} \delta_{l l''}. \eqno (13)]$

Evaluating this equation provides a convenient way to verify the following calculations.

2.4. GTO translation matrix elements

The normalized GTO radial functions are given by (Biedenharn & Louck, 1981b )

$[R_{nl}(r) = \bigg [{{2}\over{\lambda^{3/2}\pi^{1/2}}} {{(n-l-1)!}\over{(1/2)_{n}}} \bigg]^{1/2} \exp(-\rho^2/2) \rho^{l} L_{n-l-1}^{(l+1/2)}(\rho^2), \eqno (14)]$

where $[\rho^2]$ = r²/ $[\lambda]$ with scale factor $[\lambda]$ . For protein shape representations we set $[\lambda]$ = 20. Throughout this article (x)_n = x(x + 1)…(x + n − 1) = $[\Gamma]$ (x + n)/ $[\Gamma]$ (x) denotes a rising factorial and the substitution $[\Gamma]$ (n + 1/2) = $[\pi^{1/2}]$ (1/2)_n is frequently employed for numerical accuracy. The generalized Laguerre polynomials, $[L_{k}^{(\alpha)}(x)]$ , are defined by the binomial expansion (Erdelyi et al., 1953a )

$[L_{k}^{(\alpha)}(x) = \sum_{j = 0}^{k} \left({ k + \alpha \atop k - j} \right) { {{(-x)\,^j}\over{j!}}}. \eqno (15)]$

High-order polynomials may be calculated efficiently using the stable recursion

$[(k+1) L_{k+1}^{(\alpha)}(x) = (2k+\alpha+1-x) L_{k}^{(\alpha)}(x) - (k+\alpha) L_{k-1}^{(\alpha)}(x), \eqno (16)]$

along with the identities $[L_{0}^{(\alpha)}(x)]$ = 1 and $[L_{1}^{(\alpha)}(x)]$ = $[\alpha]$ + 1 − x. The GTO functions are eigenfunctions of the spherical Bessel transform (derived from Erdelyi et al., 1953b , p. 42, equation 3 therein):

$[\eqalignno{\tilde{R}_{nl}(\beta) =\hskip.2em& (-1)^{n-l-1} \left [{ 2 \lambda^{3/2} \over \pi^{1/2} } { (n-l-1)! \over (1/2)_{n} } \right] ^{1/2} \cr&\!\times\exp(-x^2/2) x^l L_{n-l-1}^{(l+1/2)} (x^2), & (17)}]$

where x² = $[\lambda\beta^{2}]$ . Here, it is convenient to use equation (15) to expand equation (17) as a power series:

$[\tilde{R}_{nl}(\beta) = \bigg [{ 4 \lambda^{3/2} \over \pi^{1/2} } \bigg]^{1/2} \sum_j X_{nlj} \exp(-x^2/2) x^{2j+l} \eqno (18)]$

where $[\textstyle\sum_j]$ serves as shorthand notation for $[\textstyle\sum_{j = 0}^{n-l-1}]$ and where the coefficients X_nlj are given by

$[X_{nlj} = \left [{ (n-l-1)! (1/2)_{n} \over 2 } \right]^{1/2} { (-1)^{n-l-j-1} \over j! (n-l-j-1) ! (1/2)_{l+j+1} }. \eqno (19)]$

Substituting equation (18) twice into equation (10) and collecting coefficients of x^2k using

$[C_k^{(nl,n'l')} = \sum_{j} \sum_{j'} \delta_{k,j+j'} X_{nlj} X_{n'l'j'} \eqno (20)]$

(where $[\delta_{kj}]$ is the Kronecker delta) gives for the GTO translation matrix elements:

$[\eqalignno{ T_{n'l',nl}^{(|m|)}(R) =\hskip.2em& \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \sum_{j = 0}^{n-l+n'-l'-2} C_j^{(nl,n'l')} \cr &\!\times {4 \over \pi^{1/2} }\int\limits_0^\infty x^{2j+l+l'} j_k(x R/\lambda^{1/2}) x^2\, {\rm d}x . &(21)}]$

Applying the relation (from Erdelyi et al., 1953b, p. 30, equation 13 therein)

$[\eqalignno{&{4 \over \pi^{1/2} } \int\limits_0^{\infty} \exp(-x^2) x^{2m+k} j_k(xy) x^2 \,{\rm d}x \cr&\quad= m! \exp(-y^2/4) (y^2 / 4)^{k/2} L_{m}^{(k+1/2)}(y^2 / 4),&(22)}]$

then gives the final analytic result

$[\eqalignno{ T_{n'l',nl}^{(|m|)}(R) =\hskip.2em& \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \sum_{j = 0}^{n-l+n'-l'-2} C_j^{(nl,n'l')} M! \cr &\!\times \exp(-R^2/4\lambda) (R^2 / 4\lambda)^{k/2} L_M^{(k+1/2)}(R^2/4\lambda) ,\cr &&(23)}]$

where M = j + (l + l′ − k)/2.

2.5. ETO translation matrix elements

The normalized ETO radial functions are given by (Biedenharn & Louck, 1981b)

$[S_{nl}(r) = \left [(2\Lambda)^3 {(n-l-1)! \over (n+l+1)!} \right]^{1/2} \exp(-\rho/2) \rho^l L_{n-l-1}^{(2l+2)}(\rho), \eqno (24)]$

where $[\rho]$ = 2 $[\Lambda]$ r with scale factor $[\Lambda]$ . We set $[\Lambda]$ = 1/2 for protein–protein electrostatic calculations (Ritchie & Kemp, 2000). Using an argument based on orthogonality, Keister & Polyzou (1997 ) recently proved that the spherical Bessel transform of these functions may be written in terms of the Jacobi polynomials, $[P_{k}^{(\mu,\nu)}(t)]$ :

$[\eqalignno{ {\tilde S}_nl(\beta) =\hskip.2em& {{2}\over{(1/2)_n}} \left [{{(n-l-1)! (n+l+1)!}\over{\pi \Lambda^3}} \right] ^{1/2} \cr &\!\times {{s^l}\over{(s^2+1)^{l+2}}} P_{n-l-1}^{(l+3/2,l+1/2)} \bigg({{s^2-1}\over{s^2+1}} \bigg) , & (25)}]$

where s = $[\beta]$ / $[\Lambda]$ . Following a similar treatment to the GTO case, the shifted series expansion (Keister & Polyzou, 1997)

$[\eqalignno{&P_{k}^{(\mu,\nu)}(t) \cr&\,\,\,\,= \sum_{j = 0}^k {{\Gamma(k+\mu+1)}\over{ \Gamma(k+\mu+\nu+1)}} {{(-1)^j \Gamma(k+j+\mu+\nu+1)}\over{j! (k-j)! \Gamma(j+\mu+1)}} \left({{1-t}\over{2}} \right)^j \cr&&(26)}]$

may be used to collect factors of 1/(s² + 1) = (1 − t)/2 to write equation (25) as a power series:

$[\tilde{S}_{nl}(\beta) = \left( {{2}\over{\pi \Lambda^3}} \right) ^{1/2} \sum_j Y_{nlj} {{s^l}\over{(s^2+1)^{l+j+2}}} ,\eqno (27)]$

where

$[Y_{nlj} = \left [{{1}\over{2}} {{(n-l-1)!}\over{(n+l+1)!}} \right] ^{1/2} {{(-1)^{j} (2n+1) (n+l+j+1)!}\over{j! (n-l-j-1)! (1/2)_{l+j+2}}}. \eqno (28)]$

Substituting equation (27) twice into equation (10) and collecting coefficients of 1/(s² + 1)^k using

$[D_{k}^{(nl,n'l')} = \sum_{j} \sum_{j'} \delta_{k,j+j'} Y_{nlj} Y_{n'l'j'} ,\eqno (29)]$

gives for the ETO translation matrix elements:

$[\eqalignno{ U_{n'l',nl}^{(|m|)}(R) =\hskip.2em& \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \sum_{j = 0}^{n-l+n'-l'-2} D_j^{(nl,n'l')} \cr &\!\times {{2}\over{\pi}} \int\limits_0^\infty {{s^{2M+k}}\over{(s^2+1)^{J+2}}} j_k(s \Lambda R) s^2 \,{\rm d}s , & (30)}]$

where M = (l + l′ − k)/2 and J = j + l + l′ + 2. It is shown in Appendix B that the remaining integral may be calculated as

$[\eqalignno{& {{2}\over{\pi}} \int\limits_0^{\infty} {{s^{2M+k}}\over{(s^2+1)^{J+2}}} j_{k}(s \Lambda R) s^2\, {\rm d} s \cr&\quad= \sum_{q = 0}^{M} \left({M \atop q} \right) {{ (-1)^{M+q}}\over{2^{J+1-q} (J+1-q)!}} (\Lambda R)^k \hat{k}_{J-k-q + 1/2}(\Lambda R) ,\cr&& (31)}]$

where $[\hat{k}_{\sigma}(z)]$ is a reduced Bessel function of the second kind. For half-integral degree, these functions may be calculated using the recurrence relations (Weniger & Steinborn, 1983):

$[{\hat k}_{1/2}(z) = \exp({-z}), \eqno (32)]$

$[{\hat k}_{3/2}(z) = (1 + z) \exp({-z}), \eqno (33)]$

and

$[{\hat k}_{n+3/2}(z) = (2n+1) {\hat k}_{n + 1/2}(z) + z^{2}{\hat k}_{n-1/2}(z). \eqno (34)]$

Thus, the ETO translation matrix elements may also be calculated analytically, although compared with the GTO basis an additional inner summation is necessary.

2.6. Non-orthogonal translation matrices

Translations of SPF expansions in both the GTO and ETO bases can be computed more economically by eliminating the inner summation on the subscript j in equations (23) and (30). This is equivalent to calculating overlap integrals that correspond to expansions of non-orthogonal radial basis functions. For example, substituting equation (18) into equation (10) and applying equation (22) directly gives the factorization

$[T_{n'l',nl}^{(|m|)}(R) = \sum_{j'} \sum_{j} X_{n'l'j'} \bar{T}_{j'l',jl}^{(|m|)}(R) X_{nlj}, \eqno (35)]$

where each $[\bar{T}_{j'l',jl}^{(|m|)}(R)]$ is an overlap integral in a non-orthogonal basis,

$[\eqalignno{ \bar{T}_{j'l',jl}^{(|m|)}(R) =\hskip.2em& \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} M!\exp({-R^2/4\lambda}) \cr &\!\times (R^2/4\lambda)^{k/2} L_{M}^{(k+1/2)}(R^2/4\lambda), &(36)}]$

now with M = j + j′ + (l + l′ − k)/2. This corresponds to expanding R_nl(r) as a sum of non-orthogonal functions, $[\bar{R}_{nl}(r)]$ :

$[R_{nl}(r) = (-1)^{n-l-1} \sum_{j} X_{nlj} \bar{R}_{jl}(r). \eqno (37)]$

It can be shown that these functions correspond to the non-orthogonal Laguerre–Gaussian basis proposed by Chiu & Moharerrzadeh (1999). With this factorization, translated expansion coefficients, a^R_nlm, in the original orthogonal basis may be calculated using the sequence

$[{\bar a}_{jlm} = \sum_n \, X_{nlj} a_{nlm}, \eqno (38)]$

$[{\bar a} ^{R}_{jlm} = \sum_{j'l'} {\bar T}_{jl,j'l'}^{(|m|)}(R) \,{\bar a}_{j'l'm}, \eqno (39)]$

$[{\bar a}^{R}_{nlm} = \sum_{j} X_{nlj} {\bar a}^{R}_{jlm}. \eqno (40)]$

In a similar manner, ETO translations may be calculated in a non-orthogonal basis by substituting equation (27) into equation (10) and collecting powers of 1/(s² + 1) directly to give

$[U_{n'l',nl}^{(|m|)}(R) = \sum_{j'} \sum_j Y_{n'l'j'} \bar{U}_{j'l',jl}^{(|m|)}(R) Y_{nlj}, \eqno (41)]$

where the non-orthogonal matrix elements, $[\bar{U}_{j'l',jl}^{(|m|)}(R)]$ , are given by

$[\eqalignno{ {\bar U}_{j'l',jl}^{(|m|)}(R) =\hskip.2em& \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \sum_{q = 0}^{M} \left ({{M}\atop{q}} \right) {{(-1)^{M+q}}\over{2^{J+1-q} (J+1-q)!}} \cr &\!\times (\Lambda R)^k \hat{k}_{J-k-q+1/2}(\Lambda R) & (42)}]$

with M = (l + l′ - k)/2 and now J = j + j′ + l + l′ + 2.

Finally, using equation (31) with M = 0 and J = l + j, it is straightforward to apply the inverse Bessel transform to equation (27) to give

$[\eqalignno{ S_{nl}(r) =\hskip.2em& \Lambda^{3/2} \sum_j Y_{nlj} {{1}\over{2^{l+j+1} (l+j+1)!}} (\Lambda r)^l \hat{k}_{j+1/2}(\Lambda r) \cr =\hskip.2em& \Lambda^{3/2} \sum_j Y_{nlj} B_{jl}(\Lambda r) , &(43)}]$

where the $[B_{jl}(\Lambda r)]$ are the well known B-functions, or BTOs (Filter & Steinborn, 1978). In other words, equation (43) shows that ETOs correspond to orthogonalized BTOs.

2.7. Correlation algorithm and complexity analysis

To illustrate an application of our approach, we demonstrate the superposition of a pair of similar proteins, A and B, with each protein initially centred at the origin and represented as a single three-dimensional function [equation (1)], e.g. steric density or electrostatic charge. Letting $[A({\bf r})]$ and $[B({\bf r})]$ represent the pair of properties to be correlated, it is then convenient to express the correlation as

$[\eqalignno{ C(\beta_1,\gamma_1,R,\alpha_2,\beta_2,\gamma_2) = \hskip.2em&\int [\hat{T}_z(-R) \hat{R}(0,\beta_1,\gamma_1) A({\bf r})] \cr &\!\times [\hat{R}(\alpha_2,\beta_2,\gamma_2) B({\bf r})] \, {\rm d}V, &(44)}]$

where $[\hat{R}(\alpha,\beta,\gamma)]$ represents an Euler rotation operator, $[\hat{R}(\alpha,\beta,\gamma)]$ $[\equiv]$ $[\hat{R}_z(\alpha)\hat{R}_y(\beta)\hat{R}_z(\gamma)]$ , and where $[\hat{T}_z(-R)]$ translates the rotated $[A({\bf r})]$ by R along the negative z axis. We could have chosen instead to translate the rotated $[B({\bf r})]$ along the positive z axis. We use an icosahedral tessellation of the sphere to generate near-regular $[(\beta,\gamma)]$ rotational samples for each molecule (Ritchie & Kemp, 1999 ). If the proteins are similar, the translational part of the search will normally be small. In any case, we seek to find the maximum of the expansion

$[\eqalignno{ C(\beta_1,\gamma_1,R,\alpha_2,\beta_2,\gamma_2) =\hskip.2em& \sum_{n l n' l' m m' m''} T^{(|m'|)}_{n'l',nl}(R) R^{(l')}_{mm'}(0,\beta_1,\gamma_1) \cr &\!\times a_{n'l'm'} R^{(l)}_{mm''}(\alpha_2,\beta_2,\gamma_2) b_{nlm''}. &(45)}]$

The real rotation matrix elements, $[R^{(l)}_{mm'}(\alpha,\beta,\gamma)]$ , may be calculated efficiently using recursion formulae (Navaza, 1990 ; Ritchie & Kemp, 1999). However, the cost of computing equation (45) as it stands scales in the order of O(N⁷) operations for each trial orientation. A much more efficient strategy is to compute the sum in stages using precalculated rotation and translation matrix elements. For example, vectors of coefficients representing different three-dimensional $[(\beta_1,\gamma_1,R)]$ orientations of molecule A may be calculated in O(N³) × [O(N) + O(N²)] = O(N⁵) operations per vector (orientation) using

$[a_{nlm}^{\beta_1\gamma_1} = \sum_{m'} R^{(l)}_{mm'}(0,\beta_1,\gamma_1) a_{nlm'} \eqno (46)]$

and

$[a_{nlm}^{R\beta_1\gamma_1} = \sum_{n'l'} T^{(|m|)}_{nl,n'l'}(-R) a_{n'l'm}^{\beta_1\gamma_1}. \eqno (47)]$

Similarly, vectors of rotated instances of molecule B may be calculated in O(N³) × O(N) = O(N⁴) operations per orientation using

$[b_{nlm}^{\beta_2\gamma_2} = \sum_{m'} R^{(l)}_{mm'}(0,\beta_2,\gamma_2) b_{nlm'}. \eqno (48)]$

The final degree of freedom is a twist rotation about the z axis, which could be applied to molecule B as

$[b_{nlm}^{\alpha_2\beta_2\gamma_2} = \sum_{m'} R^{(l)}_{mm'}(\alpha_2,0,0) b_{nlm'}^{\beta_2\gamma_2} \eqno (49)]$

$[\hphantom{b_{nlm}^{\alpha_2\beta_2\gamma_2}} = \sum_{m'} b_{nlm'}^{\beta_2\gamma_2} \cos m' \alpha_2 + b_{nl{\bar m}'}^{\beta_2\gamma_2} \sin {\bar m}' \alpha_2. \eqno (50)]$

However, it is more efficient to complete the correlation by iterating combinatorially over all pairs of the above-calculated A and B orientations in O(N) × O(N²) = O(N³) operations per iteration using

$[P_m^{\beta_1 \gamma_1 R \beta_2 \gamma_2} = \sum_{nl} a_{nlm}^{R\beta_1\gamma_1}b_{nlm}^{\beta_2\gamma_2} \eqno (51)]$

and

$[Q_m^{\beta_1 \gamma_1 R \beta_2 \gamma_2} = \sum_{nl} a_{nlm}^{R\beta_1\gamma_1} b_{nl\overline{m}}^{\beta_2\gamma_2}, \eqno (52)]$

followed by an inner iteration over $[\alpha_2]$ using

$[\eqalignno{ C(\beta_1,\gamma_1,R,\alpha_2,\beta_2,\gamma_2) = \hskip.2em&\sum_{m} P_m^{\beta_1 \gamma_1 R \beta_2 \gamma_2} \cos m \alpha_2 \cr&\! + Q_m^{\beta_1 \gamma_1 R \beta_2 \gamma_2} \sin {\bar m} \alpha_2 , & (53)}]$

which clearly scales as O(N) per trial orientation. The calculation of equation (53) for multiple angular samples, M, where M $[\ge]$ N, may be performed in O(M logM) operations by using a one-dimensional FFT. However, because all quantities here are real, when N $[\ge]$ 16, we use a single real FHT instead of the complex FFT to obtain around a 10% speed-up. Otherwise, it is faster to compute equation (53) explicitly for each value of $[\alpha_2]$ in O(MN) time. Thus, despite the relatively high cost of rotating [O(N⁴)] and translating [O(N⁵)] individual three-dimensional coefficient vectors, it can be seen that the above scheme reduces the cost of the combinatorial part of a six-dimensional correlation to just two inexpensive nested iterations, with the inner cycle costing only O(N) operations or less per trial orientation.

3. Results

3.1. Verification and numerical precision of expressions

In order to verify our implementation, numerical results from the analytic translation matrix element expressions were compared with those from a two-dimensional numerical integration of equation (3) and a one-dimensional integration of equation (10). The two-dimensional integration was calculated over a regular 200 × 200 grid in the $[[r,\cos(\theta)]]$ plane with r ranging from 0 to 50 Å. The one-dimensional integration used 200 steps in $[\beta]$ in a log-numerical scheme using

$[\eqalignno{& \int\limits_0^{\infty} \tilde{R}_{nl}(\beta) \tilde{R}_{n'l'}(\beta) j_{k}(\beta R) \beta^2 \,{\rm d}\beta \cr&\quad\simeq \sum_{\beta_{1} = 10^{-4}}^{\beta_{200} = 5\times 10^{8} } \tilde{R}_{nl}(\beta_i) \tilde{R}_{n'l'}(\beta_i) j_{k}(\beta_i R) \beta_{i}^{2} \exp({\gamma_i}) \Delta \gamma_i & (54)}]$

with $[\gamma_i]$ = $[\ln \beta_i]$ and $[\Delta \gamma_i]$ = $[\gamma_i-\gamma_{i-1}]$ . In both cases, increasing the number of integration steps or integration range gave a negligible effect on the numerical results. When calculating the angular coefficients $[A_k^{(ll'|m|)}]$ [equation (11)] for the one-dimensional integration, we found it became essential to use extended-precision arithmetic for N $[\ge]$ 16 due to the many large factorials in Wigner's formula for the 3–j symbols (Biedenharn & Louck, 1981a). Hence an array of coefficients was precalculated using 256-bit arithmetic using the GMP arbitrary precision math library (http://www.swox.com/gmp/ ). We note that Tuzan et al. (1998 ) describe a procedure for calculating 3–j symbols using only single-precision arithmetic, but this requires significantly more programming than our current approach. Nonetheless, the remaining calculations do not require such high-precision arithmetic. For example, when using 128-bit arithmetic, the analytic calculations agreed with the numerical integration results to within three decimal digits or more for all matrix elements up to N = 32 for all translations up to R = 9 Å in the GTO basis, and up to R = 90 Å in the ETO basis. As an additional check, the orthonormality of both the GTO and ETO translation matrices was verified by evaluating the orthogonality formula equation (13) at a reduced distance of $[\rho]$ = 0.1. This gave zero or unity, as appropriate, to within at least eight decimal digits up to N = 20, and to within at least two decimal digits for all translation matrices up to N = 32. Hence we are confident that our formulae and implementation are correct.

However, the above tests indicate that numerical precision can fall off significantly at high order. Table 1 presents further details on the overall computational cost and numerical precision for calculating GTO translations. Table 2 gives the corresponding results for the ETO basis. These tables show that log-numerical integration is the fastest method and that this gives better numerical precision than using 64-bit (i.e. Fortran or C `double precision') arithmetic to evaluate the analytic expressions for high-order expansions. The tables also show that the numerical accuracy of the analytic calculations falls off markedly with high-order expansions unless high-precision arithmetic is used. This is presumably due to summing many terms of very different magnitudes. Overall, it can be seen that the ETO basis gives somewhat longer calculation times and slightly larger rounding errors than the corresponding GTO functions. The final columns of each table show that the non-orthogonal expansion method is at least twice as fast as using the orthogonal analytic formulae and gives equal or better precision. To achieve double-precision programming accuracy for calculations to N = 32, Tables 1 and 2 indicate that translation matrix elements should be calculated using 160-bit and 192-bit arithmetic for the GTO and ETO bases, respectively. Our correlation algorithm has the option to use either non-orthogonal expansions for `on the fly' translations or to calculate and store orthogonal translation matrices. In the latter case, the matrix elements may be stored and subsequently used in ordinary double-precision arithmetic.

Table 1
GTO basis vector translation times and accuracy

This table lists the total number of distinct non-trivial matrix elements (integrals) for GTO translation matrices to order N, along with the time to calculate each matrix and translate a unit coefficient vector by 1 Å along the positive z axis. The remaining columns are labeled as 2D64: matrix elements calculated by two-dimensional numerical integration of equation (3) using a 200 × 200 grid in the $[(r,\theta)]$ plane, working in 64-bit arithmetic; 1D64: one-dimensional log-numerical integration of equation (54) using 200 steps in $[\ln \beta]$ and 64-bit arithmetic; A64: analytic integration [equation (23)] using 64-bit arithmetic, etc; N128: non-orthogonal analytic integration [equaions (36) and (38)–(40)] in 128-bit arithmetic, etc. All times are in CPU s on a 2 GHz Pentium Xeon processor. Numbers in brackets give a measure of the overall precision of the result. For example, an entry of (16) means that the root mean squared (RMS) error in the translated vector is no more than O(10⁻¹⁶). This is derived from the RMS deviation between a vector translated by a matrix calculated using the given number of bits of arithmetic and one for which 512-bit (effectively infinite precision) arithmetic was used.

N	Integrals	2D64	1D64	A64	A128	N128	A160	N160
6	504	0.3 (4)	0.00 (12)	0.03 (16)	0.03 (16)	0.01 (16)	0.03 (16)	0.01 (16)
9	2871	1.2 (3)	0.01 (12)	0.07 (12)	0.10 (16)	0.02 (16)	0.07 (16)	0.03 (16)
12	10374	3.6 (3)	0.04 (12)	0.21 (12)	0.19 (16)	0.10 (16)	0.20 (16)	0.09 (16)
16	38760	12.2 (3)	0.15 (11)	0.58 (9)	0.77 (16)	0.38 (16)	0.82 (16)	0.39 (16)
20	109802	37.2 (2)	0.45 (8)	1.90 (5)	2.54 (16)	1.19 (16)	2.73 (16)	1.26 (16)
25	315315	102.4 (2)	1.36 (5)	6.13 (1)	8.47 (16)	3.97 (16)	9.30 (16)	4.24 (16)
30	752928	233.1 (2)	3.65 (3)	17.41 (0)	23.38 (15)	10.83 (15)	25.78 (16)	11.72 (16)
32	1026256	315.1 (2)	5.39 (2)	24.49 (0)	33.49 (13)	15.57 (14)	36.82 (16)	16.74 (16)

Table 2
ETO basis vector translation times and accuracy

This table lists the computation times in s and gives in brackets the number of digits of precision for calculating and using ETO translation matrices to order N, as described in the legend for Table 1. For the final two columns, higher precision arithmetic was used than in Table 1.

N	Integrals	2D64	1D64	A64	A128	N128	A192	N192
6	504	0.3 (3)	0.01 (7)	0.03 (14)	0.04 (16)	0.01 (16)	0.04 (16)	0.01 (16)
9	2871	1.4 (3)	0.02 (6)	0.07 (11)	0.08 (16)	0.03 (16)	0.08 (16)	0.04 (16)
12	10374	4.3 (2)	0.05 (6)	0.17 (6)	0.21 (16)	0.10 (16)	0.24 (16)	0.12 (16)
16	38760	14.8 (3)	0.17 (5)	0.61 (1)	0.79 (16)	0.41 (16)	0.98 (16)	0.47 (16)
20	109802	45.1 (2)	0.46 (5)	1.85 (0)	2.47 (14)	1.29 (14)	3.23 (16)	1.50 (16)
25	315315	123.3 (2)	1.40 (5)	6.23 (0)	8.43 (7)	4.20 (7)	10.87 (16)	4.92 (16)
30	752928	287.5 (2)	3.89 (5)	17.17 (0)	22.80 (0)	11.24 (0)	30.35 (16)	13.42 (16)
32	1026256	391.6 (1)	5.89 (4)	24.62 (0)	33.11 (0)	16.12 (0)	44.27 (15)	19.00 (15)

3.2. Protein superposition correlation example

Fig. 1(A) shows some GTO steric density representations of a pair of globular protein domains, namely the superantigens Streptococcal pyrogenic exotoxin A1 (PDB code 1B1Z) (Papageorgiou et al., 1995 ) and the Staphylococcus aureus exotoxin SEC3 (PDB code 1 JCK) (Fields et al., 1996 ). These globular proteins have a relatively low sequence identity of 46%, but share a highly similar fold. Hence they may be superposed well by conventional least-squares fitting of conserved C_α coordinates. However, in this illustration the superposition was performed by maximizing the overlap between the respective GTO steric density expansions using correlations to N = 6. A near-identical superposition (not shown) was also achieved by correlating electrostatic charge densities in the ETO basis. Fig. 1(B) shows the corresponding backbone traces in the calculated superposition, along with the molecular surfaces from which the GTO expansions were derived. It can be seen that the high-order expansions capture the detailed shape of each protein remarkably well, although the low-order expansions still encode sufficient steric information to allow a very good global superposition to be calculated. The superposition shown was calculated by searching over some 21 × 10⁶ trial orientations generated from 162 $[(\beta,\gamma)]$ icosahedral angular samples for each protein, 128 twist samples in $[\alpha_2]$ , and 40 distance steps of 0.25 Å. After each cycle over the twist angle, the best overlap score was saved along with the corresponding coordinate values in a block of memory sufficient to store around 10⁶ orientations. As the search progressed, this memory block was periodically sorted and culled. Hence there is essentially no practical limit on the number of trial orientations that may be evaluated and subsequently stored. In this example, the total memory usage never exceeded 15 Mbyte and the overall correlation time was under 8 s on a 2 GHz Pentium Xeon processor. This corresponds to a rate of 6 × 10⁶ trial orientations s⁻¹, or just a few hundred arithmetic operations per orientation.

Figure 1
Illustration of the GTO shape representation and superposition of a pair of globular proteins, the superantigens SpeA and SEC3. (A) From top left to bottom right: the steric density functions of SpeA and SEC3 shown at expansion orders N = 6, 12, 16, 20, 25 and 30. Each pair is in the same superposed orientation, separated horizontally for clarity, with SpeA on the left and SEC3 on the right. For visualization, each surface shape was contoured from the three-dimensional density function as described previously (Ritchie, 2003

). Expansions to N = 32 are visually almost indistinguishable from the N = 30 expansions, and are not shown here. (B) The corresponding backbone traces of the superposition with SpeA in yellow and SEC3 in cyan (top), the separated backbones with SpeA on the left and SEC3 on the right (centre), and the original molecular surfaces from which the GTO density representations were derived (bottom).

4. Discussion

Conceptually, our approach follows in the spirit of Crowther's original spherical harmonic plus spherical Bessel method of correlating Patterson maps (Crowther, 1972 ). However, by replacing the Bessel functions with a set of orthonormal radial basis functions which are eigenfunctions of the spherical Bessel transform, our SPF approach extends this rotation-centric representation in a way which also allows the effect of translations to be calculated analytically. We have presented efficient spherical Bessel transform methods of calculating translation matrix elements for both GTO and ETO radial basis functions. Although it is well known that the Laguerre–Gaussian functions are eigenfunctions of the spherical Bessel transform, we believe the methods of calculation and the expressions presented here, especially for the ETO basis, are straightforward yet novel. However, it is necessary to call upon a significant body of special function theory to derive analytic forms for the SPF translation matrices, and the resulting expressions require considerably more programming to implement than the FFT. Nonetheless, the computational `pay-off' that follows is that once an initial transform into the SPF basis has been calculated for each molecular property of interest, a full six-dimensional correlation search may be performed using only the initial expansion coefficients. Because all subsequent expressions are entirely real, no inverse transform is required. On the other hand, computing high-order matrix elements is relatively expensive. However, this overhead can be alleviated by calculating non-orthogonal expansions or by precalculating and storing orthogonal translation matrices for subsequent use.

There are also close parallels between our technique and the complex five-dimensional FFT-based EM density fitting approach of Kovacs et al. (2003). However, their approach is currently limited to relatively low-order (L = 16) spherical harmonic expansions in order to stay within the memory limits of contemporary 32-bit computers, and their use of discrete radial envelope functions requires the translational component of the correlation to be computed by two-dimensional numerical integration (Kovacs et al., 2003). Unlike the translation matrix elements used here, these integrals cannot be pre-computed as a one-off cost because they depend intrinsically on the shapes to be compared.

In general, FFT-based correlation algorithms use a power of two for both the polynomial order of the representation and the number of search increments in each degree of freedom in the correlation. For example, for angular searches, a rotational FFT step size of at least 64 is desirable in order to give reasonably small search increments (e.g. 360°/64 $[\simeq]$ 5.6°). Similarly, Cartesian FFT docking algorithms often use grids of 64³ or 128³ elements in order to accommodate all possible translations of one protein about the other with a sub-Å step size (Katchalski-Katzir et al., 1992). Hence, unless a truncation technique is employed (Segal & Eisenstein, 2005 ), most FFT docking algorithms implicitly use a very high-order polynomial representation for each protein. In contrast, in the SPF approach the polynomial order and search step size may easily be varied independently. In fact, because our polynomial order, N, is generally significantly less than the corresponding FFT order (say M) for any given coordinate, we typically have O(MN) ≲ O(M lnM). Thus, in practice, it is not a disadvantage that our correlation algorithm uses an FFT (actually a real FHT) in at most just one of the search coordinates. Indeed, by calculating and reusing vectors of rotated and translated expansion coefficients, the complexity analysis given above shows that the inner iteration of the combinatorial part of a six-dimensional search scales as O(N) or better per trial orientation.

Given the modest memory overhead of under 15 Mbyte of our low-order N = 6 superposition calculation (even our most exhaustive high-order N = 32 docking correlations use under 200 Mbyte), and considering that our superposition correlation time compares very favourably with the timings in Table I of Kovacs et al. (2003), it would seem that the SPF approach could offer significant computational advantages to both EM and MR applications. The straightforward analytic method of calculating ETO translation matrix elements presented here could also provide a useful new alternative to using BTO expansions to calculate STO overlap integrals in quantum mechanics. In principle, GTO and ETO expansions should give comparable results if suitable scale factors are chosen and if a sufficient number of terms are used. However, in practice, the ETOs decay much more gradually with distance than the GTOs. In our experience, the GTO basis works well for both superposing and docking globular protein domains, whereas the ETO basis is more appropriate for rapidly calculating electrostatic interactions (Ritchie & Kemp, 2000). However, ETOs might also be a good choice for correlating very large low-resolution macromolecular structures. The results here indicate that either type could be used to calculate e.g. real-space electron density correlations. Although the example presented here shows that SPF expansions have good computational properties, it is not yet clear how the algorithm might perform when the translational component is large or when very high angular resolution is required. It would also be interesting to test the algorithm on a wider range of examples. Hence a more extensive study of the utility of the approach is planned. Developing a fast and robust six-dimensional correlation algorithm would be useful in several areas including e.g. molecular replacement, docking proteins into electron microscopy density maps, and searching the PDB for structural homologues.

5. Conclusion

Spherical Bessel transform formulae have been given for the analytic calculation of translation matrices for SPF expansions using GTO and ETO radial basis functions. To our knowledge, the formulae presented to translate ETO basis functions and to relate ETOs to BTOs are novel. It has been shown that translation matrices for both the GTO and ETO bases may be calculated accurately up to N = 32 using an extended-precision arithmetic library. The calculation may be accelerated by over a factor of two by using non-orthogonal translation matrix elements, which give equal or better precision. Although our SPF correlation algorithm uses an FFT in only one of the search coordinates, this is not a practical disadvantage for relatively low-resolution superposition and docking correlations. On the contrary, the SPF approach is very economical in both CPU time and computer memory. By ordering the sub-expressions according to computational cost, an exhaustive six-dimensional superposition correlation of a pair of protein shapes can be performed in a matter of seconds on a contemporary personal computer. This demonstrates the utility of the SPF approach. It is proposed that the techniques described here, which have been implemented in the protein docking and superposition program Hex (http://www.csd.abdn.ac.uk/hex/ ), may be of value in other fields, such as MR and EM, that need to calculate real-space six-dimensional correlations.

APPENDIX A

Derivation of equation (10)

The two-coordinate systems r = $[(r,\theta,\phi)]$ and r′ = r − T = $[(r',\theta',\phi)]$ may be related functionally by multiplying the vector equation

$[{\bf r} = {\bf T} + {\bf r}' \eqno (55)]$

by an arbitrary complex vector $[i {\bf k}]$ and by exponentiating each side to give

$[\exp({ i {\bf k}\cdot {\bf r} }) = \exp({ i {\bf k}\cdot {\bf T} } )\,\exp({ i {\bf k}\cdot {\bf r}' }). \eqno (56)]$

Then, substituting k = $[(\beta,\Theta,\Phi)]$ , r = $[(r,\theta,\phi)]$ , r′ = $[(r',\theta',\phi')]$ , and T = $[(R,\gamma,\delta)]$ into Raleigh's plane wave equation (Bransden & Joachain, 1997 )

$[\exp({ i {\bf k}\cdot{\bf r} }) = 4 \pi\sum_{l = 0}^{\infty}\sum_{m = -l}^{l} i^l j_{l}(\beta r) Y_{lm}(\Theta,\Phi)^* Y_{lm}(\theta,\phi), \eqno (57)]$

gives

$[\eqalignno{& \sum_{pq} i^p j_{p}(\beta r) Y_{pq}(\Theta,\Phi)^* Y_{pq}(\theta,\phi) \cr&\quad= 4 \pi \sum_{kj} \sum_{st} i^{k + s} j_{k}(\beta R) Y_{kj}(\Theta,\Phi)^* Y_{kj}(\gamma,\delta) j_{s}(\beta r') \cr&\qquad\times Y_{st}(\Theta,\Phi)^* Y_{st}(\theta',\phi') ,&(58)}]$

where $[Y_{pq}(\theta,\phi)]$ etc. are complex spherical harmonics. Here, each summation has an infinite range, subject to $[|q| \,\lt\, p]$ , etc. Using the identity $[Y_{st}(\Theta,\Phi)^*]$ = $[(-1)^t Y_{s\overline{t}}(\Theta,\Phi)]$ , multiplying both sides by $[Y_{lm}(\Theta,\Phi)]$ , and integrating over $[(\Theta,\Phi)]$ gives

$[\eqalignno{ &j_{l}(\beta r) Y_{lm}(\theta,\phi) \cr&\quad= 4 \pi \sum_{kj} \sum_{st} i^{k+s-l} j_{k}(\beta R) Y_{kj}(\gamma,\delta) j_{s}(\beta r') Y_{st}(\theta',\phi') (-1)^{t} \cr &\qquad\times \int Y_{lm}(\Theta,\Phi) Y_{s\overline{t}}(\Theta,\Phi) Y_{kj}(\Theta,\Phi)^* \,{\rm d} \Omega . & (59)}]$

Substituting the standard formula for Gaunt's integral on the right gives

$[\eqalignno{ &j_{l}(\beta r) Y_{lm}(\theta,\phi)\cr&\quad = 4 \pi \sum_{kj} \sum_{st} i^{k+s-l} j_{k}(\beta R) Y_{kj}(\gamma,\delta) j_{s}(\beta r') Y_{st}(\theta',\phi') (-1)^{t+j}\cr &\qquad\times \bigg [{{(2l+1)(2s+1) (2k+1)}\over{4 \pi}} \bigg] ^{1/2} \bigg ({l \, s \, k \atop 0 \, 0 \, 0 } \bigg) \bigg ({l \, s \, k \atop m \, \overline{t} \, \overline{j} } \bigg). & (60)}]$

Now, the first 3–j symbol vanishes when l + s + k is odd, so the phase factor for the surviving terms is always real. Furthermore, the second 3–j symbol vanishes unless m + $[\overline{t}]$ + $[\overline{j}]$ = 0. Hence the expression may be reduced to a triple infinite sum by substituting j = m − t:

$[\eqalignno{ &j_{l}(\beta r) Y_{lm}(\theta,\phi) \cr&= 4 \pi \!\!\sum_{k}\!\!\sum_{st}\! (-1)^{(k+s-l)/2} j_{k}(\beta R) Y_{k,m-t}(\gamma,\delta) j_{s}(\beta r') Y_{st}(\theta',\phi') (-1)^m \cr &\quad\,\times\bigg [{{(2l+1)(2s+1) (2k+1)}\over{4 \pi}} \bigg] ^{1/2}\! \bigg ({l \, s \, k \atop 0 \, 0 \, 0 } \bigg)\! \bigg ({l \quad s \quad k \atop m \, \overline{t} \ t-m } \bigg) . & (61)}]$

When the translation is only in the z direction, $[\delta]$ = 0 and either $[\gamma]$ = 0 or $[\gamma]$ = $[\pi]$ , and $[\phi']$ = $[\phi]$ , which entails t = m. This allows the summation over t to be eliminated, and the expression becomes valid for real as well as complex spherical harmonics. For a positive z translation, the harmonic Y_k,0( $[\gamma]$ = 0, 0) reduces to

$[Y_{k,0} (0,0) = \bigg [{{2k+1}\over{4 \pi}} \bigg] ^{1/2}. \eqno (62)]$

Then, collecting the angular factors, relabelling $[s\rightarrow l']$ , and restricting the range of k according to the triangle rule for the 3–j symbols, gives the plane wave addition theorem:

$[j_{l}(\beta r) Y_{lm}(\theta,\phi) = \sum_{l'} \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} j_{k}(\beta R) j_{l'}(\beta r') Y_{l'm}(\theta',\phi), \eqno (63)]$

where the angular coefficient $[A_k^{(ll'|m|)}]$ is given by

$[\eqalignno{ A_k^{(ll'|m|)} = \hskip.2em&(-1)^{ (k+l'-l)/2 + m} (2k+1) \big [{(2l+1)(2l'+1)} \big] ^{1/2} \cr&\times\bigg ({l \, l' \, k \atop 0 \, 0 \, 0 } \bigg) \bigg ({l \ l' \ k \atop m \, \overline{m} \ 0 } \bigg) . & (64)}]$

For a negative translation, an additional factor of $[(-1)^{l-l'}]$ appears in the above expression because Y_k,0( $[\gamma]$ = $[\pi]$ , 0) = (-1)^k Y_k,0 (0,0) and because [k+l'-l] is even.

Now, multiplying both sides of equation (63) by $[(2/\pi)^{1/2} \tilde{R}_{nl}(\beta) \beta^2]$ and integrating over $[\beta]$ [i.e. applying the inverse Bessel transform, equation (9)] gives

$[\eqalignno{ R_{nl}(r) Y_{lm}(\theta,\phi) =\hskip.2em& (2 / \pi)^{1/2} \sum_{l'} \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \cr &\times \int\limits_0^{\infty} \tilde{R}_{nl}(\beta) j_{k}(\beta R) j_{l'}(\beta r') \beta^2 \,{\rm d}\beta \, Y_{l'm}(\theta',\phi).\cr&&(65)}]$

Then, multiplying each side by $[R_{n'j'}(r') Y_{j'm'}(\theta',\phi)]$ and integrating over all space in the corresponding variables gives

$[\eqalignno{ T_{n'j',nl}^{(|m|)}(R) =\hskip.2em& (2 / \pi)^{1/2} \sum_{l'} \delta_{j'l'} \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \cr &\times \int\limits_0^{\infty} \int\limits_0^{\infty} R_{n'l'}(r') j_{l'}(\beta r') \tilde{R}_{nl}(\beta) j_{k}(\beta R) \beta^2 \,{\rm d}\beta \, r^{\prime 2} \,{\rm d}r' . \cr&&(66)}]$

Finally, recognizing the integral in [r'] as the spherical Bessel transform of $[R_{n'l'}(r')]$ , the result reduces to

$[T_{n'l',nl}^{(|m|)}(R) = \sum_{k = |l-l'|}^{l+l'} A_k^{(ll'|m|)} \int\limits_0^{\infty} \tilde{R}_{n'l'}(\beta) \tilde{R}_{nl}(\beta) j_{k}(\beta R) \beta^2 \,{\rm d}\beta. \eqno (67)]$

This generalizes the expression given by Danos & Maximon (1965 ) for translating multipole expansions to the more general case for arbitrary orthonormal radial functions, R_nl(r).

APPENDIX B

Derivation of equation (31)

Following Weniger & Steinborn (1983), but correcting a mistake in their working [here, equation (75) onwards], a closed-form expression for the integral

$[I = \int\limits_0^{\infty} {{x^{2m+k}\over (x^2+1)^{n+2}}}\,\,j_k(x y) x^2 \,{\rm d} x \eqno (68)]$

may be obtained with the help of the basic relation (Erdelyi et al., 1953b, p. 24, equation 20 therein):

$[ \int\limits_0^{\infty} {{x^{\nu+1}}\over{(x^2+a^2)^{\mu+1}}} J_{\nu}(xy)\, {\rm d} x = {{a^{\nu-\mu} y^{\mu}}\over{2^{\mu} \Gamma(\mu+1)}} K_{\nu-\mu}(ay), \eqno (69)]$

where $[J_{\nu}(z)]$ is a general Bessel function of the first kind, and $[K_{\sigma}(z)]$ is a modified Bessel function of the second kind. When the degree $[\sigma]$ is half-integral, which is the case here, $[K_{\sigma}(z)]$ has a closed form:

$[K_{n+1/2}(z) = (\pi / 2z)^{1/2} \exp({-z}) \sum_{m = 0}^{n} \bigg ({{1}\over {2z}} \bigg)^{m} {{(n+m)!}\over{m! (n-m)!}}. \eqno (70)]$

However, this function has a singularity at the origin, so it is useful to use instead the reduced Bessel function

$[\hat{k}_{n+1/2}(z) = (2 z / \pi)^{1/2} z^{n} K_{n+1/2}(z), \eqno (71)]$

which is well behaved for all z. Equation (69) can then be cast in the desired form using the standard relation (Hochstadt, 1971):

$[(1/z\, \partial/{\partial z})^m z^{\sigma} J_{\sigma}(tz) = t^m z^{\sigma-m} J_{\sigma-m}(tz). \eqno (72)]$

Hence putting $[\sigma]$ = $[\nu]$ + m and applying $[(1/y\, \partial/{\partial y})^m y^{\nu+m}]$ to each side of equation (69) gives

$[\eqalignno{&\int\limits_0^{\infty} {{x^{\nu+2m+1}}\over{(x^2+a^2)^{\mu+1}}} J_{\nu}(xy)\, {\rm d} x \cr&\quad= {{a^{\nu+m-\mu}}\over{y^{\nu} 2^{\mu} \Gamma(\mu+1)}} (1/y\, \partial/{\partial y})^m y^{\nu+m+\mu} K_{\nu+m-\mu}(ay),\cr& & (73)}]$

or in terms of the reduced Bessel function [equation (71)]

$[\eqalignno{& \int\limits_0^{\infty} {{x^{\nu+2m+1}}\over{(x^2+a^2)^{\mu+1}}} J_{\nu}(xy)\, {\rm d} x \cr&\quad= (\pi / 2)^{1/2} {{1}\over{a^{2\mu} y^{\nu} 2^{\mu} \Gamma(\mu+1)}} (1/y\, \partial/{\partial y})^m (ay)^{2\mu} \hat{k}_{\nu+m-\mu}(ay). \cr&& (74)}]$

It can be shown (Weniger & Steinborn, 1983) that

$[\eqalignno{& (1/z\, \partial/{\partial z})^m z^{\lambda} \hat{k}_{\sigma}(z) \cr&\quad= (-1)^m \sum_{q = 0}^{m} \bigg ({m \atop q} \bigg) 2^q (-\lambda/2)_q z^{\lambda-2q} \hat{k}_{\sigma-m+q}(z) .& (75)}]$

Hence setting $[\lambda]$ = 2 $[\mu]$ and substituting equation (75) into equation (74) gives

$[\eqalignno{& \int\limits_0^{\infty} {{x^{\nu+2m+1}}\over{(x^2+a^2)^{\mu+1}}} \,\,J_{\nu}(xy)\, {\rm d} x \cr&\,= (\pi / 2)^{1/2} {{a^{2m-2\mu} (-1)^m }\over{2^{\mu} \Gamma(\mu+1)}} \sum_{q = 0}^{m} \bigg ({m \atop q} \bigg) 2^q (-\mu)_q y^{2\mu-2q-\nu} \hat{k}_{\nu-\mu+q}(ay) . \cr&& (76)}]$

Here, $[\mu]$ is an integer so the rising factorial may be recast as

$[(-\mu)_q = (-1)^q \mu! / (\mu-q)!. \eqno (77)]$

We also have $[\mu\, \gt\, \nu]$ + 2m and so to keep the degree of the Bessel function positive it is convenient to use the identity

$[\hat{k}_{-\sigma}(z) = z^{-2\sigma} \hat{k}_{\sigma}(z). \eqno (78)]$

Some further working then gives

$[\eqalignno{& \int\limits_0^{\infty} {{x^{\nu+2m+1}}\over{(x^2+a^2)^{\mu+1}}} \,\,J_{\nu}(xy) \,{\rm d} x \cr&\quad= (\pi / 2)^{1/2} {{a^{2\nu+2m-2\mu}}\over{2^{\mu}}} \sum_{q = 0}^{m} \bigg ({m \atop q} \bigg) {{(-1)^{m+q} 2^q}\over{ (\mu-q)!}} y^{\nu} \hat{k}_{\mu-\nu-q}(ay) .\cr&& (79)}]$

Finally, putting $[\mu]$ = n + 1, $[\nu]$ = k + 1/2, a = 1, and replacing J_k+1/2(xy) by the corresponding spherical Bessel function gives

$[\eqalignno{& \int\limits_0^{\infty} {{x^{2m+k}}\over{(x^2+1)^{n+2}}} \,\,j_{k}(x y) x^2 \,{\rm d} x \cr&\quad= {{\pi}\over{2}} \sum_{q = 0}^{m} \bigg ({m \atop q} \bigg) {{ (-1)^{m+q}}\over{2^{n+1-q} (n+1-q)!}} y^k \hat{k}_{n-k-q+1/2}(y) .\cr&& (80)}]$

Acknowledgements

Part of this work was funded by the BBSRC (Grant Ref. 1/B10454).

References

Arakane, F. & Matsuoka, O. (1999). Int. J. Quant. Chem. 66, 273–279. Web of Science CrossRef Google Scholar
Barnett, M. P. & Coulson, C. A. (1951). Philos. Trans. R. Soc. London A, 243, 221–249. CrossRef Google Scholar
Biedenharn, L. C. & Louck, J. C. (1981a). Angular Momentum in Quantum Physics. Reading, MA: Addison-Wesley. Google Scholar
Biedenharn, L. C. & Louck, J. C. (1981b). The Racah–Wigner Algebra in Quantum Theory. Reading, MA: Addison-Wesley. Google Scholar
Boys, S. F. (1950). Proc. R. Soc. London A, 200, 542–554. CrossRef CAS Google Scholar
Bracewell, R. N. (1999). The Fourier Transform and its Applications, 3rd ed. New York: McGraw-Hill. Google Scholar
Bransden, B. H. & Joachain, C. J. (1997). Introduction to Quantum Mechanics. Harlow, UK: Addison Wesley Longman. Google Scholar
Chang, G. & Lewis, M. (1997). Acta Cryst. D53, 279–289. CrossRef CAS Web of Science IUCr Journals Google Scholar
Chiu, L. Y. C. & Moharerrzadeh, M. (1999). Int. J. Quant. Chem. 73, 265–273. Web of Science CrossRef CAS Google Scholar
Cowtan, K. (1998). Acta Cryst. D54, 750–756. Web of Science CrossRef CAS IUCr Journals Google Scholar
Cowtan, K. (2001). Acta Cryst. D57, 1435–1444. Web of Science CrossRef CAS IUCr Journals Google Scholar
Crowther, R. A. (1972). The Molecular Replacement Method, edited by M. G. Rossmann, pp. 173–178. New York: Gordon and Breach. Google Scholar
Danos, M. & Maximon, L. C. (1965). J. Math. Phys. 6, 766–778. CrossRef Web of Science Google Scholar
Erdelyi, A., Magnus, W., Oberhettinger, F. & Tricomi, F. G. (1953a). Higher Transcendental Functions, Vol 2. New York: McGraw-Hill. Google Scholar
Erdelyi, A., Magnus, W., Oberhettinger, F. & Tricomi, F. G. (1953b). Tables of Integral Transforms, Vol. 2. New York: McGraw-Hill. Google Scholar
Fields, B. A., Malchiodi, E. L., Li, H., Ysern, X., Stauffacher, C. V., Schlievert, P. M., Karjalainen, K. & Mariuzza, R. A. (1996). Nature (London), 384, 188–192. CrossRef CAS PubMed Web of Science Google Scholar
Filter, E. & Steinborn, O. (1978). Phys. Rev. A, 18, 1–11. CrossRef CAS Web of Science Google Scholar
Filter, E. & Steinborn, O. (1980). J. Math. Phys. 21, 2725–2736. CrossRef Web of Science Google Scholar
Frank, J. (2002). Ann. Rev. Biophys. Biolmol. Struct. 31, 303–319. Web of Science CrossRef CAS Google Scholar
Fujinaga, M. & Read, R. J. (1987). J. Appl. Cryst. 20, 517–521. CrossRef Web of Science IUCr Journals Google Scholar
Guseinov, I. I. (2001). Int. J. Quant. Chem. 81, 126–129. CrossRef CAS Google Scholar
Guseinov, I. I. (2002). Int. J. Quant. Chem. 90, 114–118. Web of Science CrossRef CAS Google Scholar
Hierse, W. & Oppeneer, P. M. (1994). Int. J. Quant. Chem. 52, 1249–1265. CrossRef CAS Web of Science Google Scholar
Hochstadt, H. (1971). The Functions of Mathematical Physics. New York: Wiley. Google Scholar
Jamrog, D. C., Zhang, Y. & Philips, G. N. Jr (2003). Acta Cryst. D59, 304–314. Web of Science CrossRef CAS IUCr Journals Google Scholar
Katchalski-Katzir, E., Shariv, I., Eisenstein, M., Friesem, A. A. & Aflalo, C. (1992). Proc. Natl Acad. Sci. 89, 2195–2199. CrossRef PubMed CAS Web of Science Google Scholar
Keister, B. D. & Polyzou, W. N. (1997). J. Comput. Phys. 134, 231–235. CrossRef CAS Web of Science Google Scholar
Kissinger, C. R., Gelhaar, D. K. & Fogel, D. B. (1999). Acta Cryst. D55, 484–491. Web of Science CrossRef CAS IUCr Journals Google Scholar
Kovacs, J. A., Chacon, P., Cong, Y., Metwally, E. & Wriggers, W. (2003). Acta Cryst. D59, 1371–1376. Web of Science CrossRef CAS IUCr Journals Google Scholar
Navaza, J. (1987). Acta Cryst. A43, 645–653. CrossRef Web of Science IUCr Journals Google Scholar
Navaza, J. (1990). Acta Cryst. A46, 619–620. CrossRef Web of Science IUCr Journals Google Scholar
Navaza, J. (2001). Acta Cryst. D57, 1367–1372. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ozdogan, T., Orbay, M. & Degirmenci, S. (2005). J. Math. Chem. 37, 27–35. Web of Science CrossRef Google Scholar
Papageorgiou, A. C., Collins, C. M., Gutman, D. M., Kline, J. B., O'Brien, S. M., Tranter, H. S. & Acharya, K. R. (1995). EMBO J. 18, 9–21. Web of Science CrossRef Google Scholar
Read, R. J. & Schierbeek, A. J. (1988). J. Appl. Cryst. 21, 490–495. CrossRef CAS Web of Science IUCr Journals Google Scholar
Rico, J. F., López, R., Ema, I. & Ramírez, G. (2005). J. Comput. Chem. 26, 846–855. Web of Science CrossRef PubMed CAS Google Scholar
Ritchie, D. W. (2003). Proteins Struct. Func. Genet. 52, 98–106. Web of Science CrossRef CAS Google Scholar
Ritchie, D. W. & Kemp, G. J. L. (1999). J. Comput. Chem. 20, 383–395. Web of Science CrossRef CAS Google Scholar
Ritchie, D. W. & Kemp, G. J. L. (2000). Proteins Struct. Func. Genet. 39, 178–194. CrossRef CAS Google Scholar
Rose, M. E. (1957). Elementary Theory of Angular Momentum. New York: Wiley. Google Scholar
Roseman, A. M. (2000). Acta Cryst. D56, 1332–1340. Web of Science CrossRef CAS IUCr Journals Google Scholar
Rossmann, M. G. (1990). Acta Cryst. A46, 73–82. CrossRef CAS Web of Science IUCr Journals Google Scholar
Rossmann, M. G. (2000). Acta Cryst. D56, 1341–1349. Web of Science CrossRef CAS IUCr Journals Google Scholar
Rossmann, M. G. (2001). Acta Cryst. D57 1360–1366. Web of Science CrossRef CAS IUCr Journals Google Scholar
Segal, D. & Eisenstein, M. (2005). Proteins Struct. Func. Bioinf. 59, 580–591. Web of Science CrossRef CAS Google Scholar
Sheriff, S., Klei, H. E. & Davis, M. E. (1999). J. Appl. Cryst. 32, 98–101. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sussman, J. L., Lin, D., Jiang, J., Manning, N. O., Prilusky, J., Ritter, O. & Abola, E. E. (1998). Acta Cryst. D54, 1078–1084. Web of Science CrossRef CAS IUCr Journals Google Scholar
Tuzan, R. E., Burkhardt, P. & Secrest, D. (1998). Comput. Phys. Commun. 112, 112–148. Web of Science CrossRef Google Scholar
Vagin, A. A. & Isupov, M. N. (2001). Acta Cryst. D57, 1451–1456. Web of Science CrossRef CAS IUCr Journals Google Scholar
Vakser, I. A. (1995). Protein Eng. 8, 371–377. CrossRef CAS PubMed Web of Science Google Scholar
Weniger, E. J. & Steinborn, E. O. (1983). J. Chem. Phys. 78, 6121–6132. CrossRef CAS Web of Science Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 38| Part 5| October 2005| Pages 808-818

doi:10.1107/S002188980502474X