Determination of absolute structure using Bayesian statistics on Bijvoet differences

A description is given of a maximum-likelihood approach to absolute structure determinations of biologically active molecules.


Introduction
Bijvoet, Peerdeman and van Bommel were the first to demonstrate that the absolute configuration of a chiral molecule could be determined by X-ray crystallography (Bijvoet et al., 1951). Their method was based on the complex resonant scattering contributions to the atomic scattering factors that make the intensities of Friedel-related reflections (or their symmetry equivalents) different. This difference in intensity (the 'Bijvoet difference') depends both on the atom types present in the molecule and the wavelength of the radiation used (Flack & Shmueli, 2007). The concept of 'absolute configuration' has since been generalized to 'absolute structure' to include cases where the polarity of the structure rather than the absolute configuration is determined (Jones, 1984;Glazer & Stadnicka, 1989).
Traditionally absolute structure determination was based on analysis of Cu K data collected on a diffractometer with a point detector for compounds containing atom types heavier than phosphorus. Currently, most small-molecule structure determinations are based on data collected on diffractometers equipped with CCD detectors using Mo K radiation. The impact of this change is that often a more accurate, highly redundant and complete data set is obtained, which, however, often contains a weaker resonant scattering signal.
There exists a significant interest in the determination of the absolute configuration of biologically active molecules (van der Helm & Hossain, 1987). Unfortunately, many molecules of interest do not contain atoms heavier than sulfur. In the past, this problem was solved with the introduction of a heavier atom in the structure, e.g. with the addition of HBr (Spek, 1976). The current trend is to attempt absolute structure determination on the native compound, even when no atoms heavier than oxygen are present.
Over time a number of methods for the determination of the absolute structure have been proposed.
The most straightforward way of establishing the absolute structure of a small enantiopure molecule is to refine both enantiomers separately, subsequently select the absolute structure with the lowest crystallographic R factor and test for the statistical significance of the R-factor difference. The latter is commonly done with the Hamilton test (Hamilton, 1965).
A much more sensitive method (Zachariasen, 1965;Engel, 1972) is to select a subset of reflections from the measured data that are most sensitive to the absolute structure (relatively weak reflections with a large Bijvoet difference), and compare the calculated Bijvoet differences with the observed differences. Just by comparing the signs of these differences, the absolute structure can often be established even if the difference in R factor is inconclusive. Although the absolute structure can be determined using this method, it is not easy to quantify the degree of certainty of the assignment. Le Page et al. (1990) present a method to accompany an absolute structure determined in this way by a calculation of the probability that the absolute structure should be inverted. For this calculation use a binomial distribution. This method has not found widespread use, and therefore its performance is difficult to assess.
Another variation on this method is used by the Bijvoet program in the DIRDIF program suite (Beurskens et al., 1980;Beurskens et al., 1999). This program uses a weighted average of the signs of the Bijvoet difference (B). This method can be very successful, but it needs a carefully selected subset of Bijvoet pairs to be effective. The absolute structure assignment using this calculation is accompanied by a standard uncertainty, but it is very hard to establish the statistical correctness of this value as it relies on distributions being Gaussian, and disregards the careful selection of the reflection subset. Also, no difference can be seen in this calculation between a racemic twin and a weak resonant scattering signal: both will result in smaller absolute values of B and larger standard uncertainties. Rogers (1981) was the first to introduce a parameter that can be refined as part of the least-squares refinement. This parameter encodes the 'strength' and sign of the measured resonant scattering signal measured in units of f 00 , the imaginary component of the complex atomic scattering factor.
The Rogers parameter was soon superseded by the Flack x parameter (Flack, 1983). The Flack x parameter encodes the relative abundance of the two components in an inversion twin. The value of the Flack x parameter can be determined using a full-matrix least-squares procedure [e.g. with the TWIN/BASF instructions in SHELXL97 (Sheldrick, 1997)]. A reasonable estimate of the Flack x parameter can be obtained by determining the parameter separately; this is automatically performed for all non-centrosymmetric structures in the SHELX97 package. Since the Flack x parameter can correlate with the atomic coordinates, especially for structures in space groups that do not have a fixed origin, the estimate can be inaccurate if its value deviates significantly from zero .
Since the value of the Flack x parameter is the result of a least-squares refinement, its standard uncertainty can be derived from the covariance matrix. This standard uncertainty can be used to quantify the degree of confidence in the proposed absolute structure. Flack & Bernardinelli (2000) discuss criteria for the reliability of the absolute structure assignment based on the standard uncertainty in the Flack x parameter value. Their analysis, starting from only the standard uncertainty, has to assume that the distribution is normal also in its tails. The paper does not distinguish between the probability of obtaining the absolute structure given the observations and the probability of the observations given the absolute structure. The Bayesian prior relating the two probabilities is ignored. The Flack & Bernardinelli (2000) method does not result in a quantitative statement about the absolute structure assignment. Parsons & Flack (2004) recently introduced a variation of the refinement of the Flack x parameter. Their method relies on the careful determination of a few selected Bijvoet differences [as the ratio ðI þ À I À Þ=ðI þ þ I À Þ] which can either be obtained directly from a good redundant data set or by carefully adding some extra observations. Parsons & Flack (2004) show that this method increases the sensitivity of the absolute structure determination. Dittrich et al. (2006) recently reported advances in absolute structure determinations made using 'invarioms'. Invarioms are aspherical scattering factors that take into account electron density deformations. Using invarioms instead of the normal spherical scattering factors can improve figures of merit as well as the standard uncertainties of all refined parameters. Experience has shown that by using invarioms, the standard uncertainty in the value of the Flack x parameter can be significantly reduced, and that the calculated value of the Flack x parameter is frequently closer to 0.0.
In the next section, we introduce a new way to determine the absolute structure. The method applies Bayesian statistics to the Bijvoet differences. The result of this approach is a series of probabilities for different hypotheses of the absolute structure. The solution with the highest probability can be determined, and this can be used to map the results to a value with standard uncertainty that can be directly compared with the value of the Flack x parameter.

Theory and methods
For each Bijvoet pair of reflections, we can define Here, F 2 c are calculated intensities and F 2 o are observed intensities. If we assume completely independent observations of the two reflection intensities, Now, we can define a variable z as follows: If the absolute structure of the model for calculation of the structure factors is correct, and we assume the calculated intensities to be correct (i.e. they do not carry a standard uncertainty) the probability distribution of z is a standard normal Gaussian 1 Based on all pairs of observations, we can now calculate the probability of the measured data, given the fact that the absolute structure is correctly specified in the model (the correct absolute structure is noted as the condition y = 0, it will become clear later in the paper why this notation is chosen): In our case, we can use this theorem to invert our probability given above: This value cannot be computed, as pðobservationsÞ (the probability of obtaining the current observations) is unknown. But to be able to make the absolute structure assignment, we would like to calculate the ratio of pðy ¼ 0 j observationsÞ and the similar term for the opposite absolute structure, designated as pðy ¼ 1 j observationsÞ. The term of pðobservationsÞ disappears in the calculation of this ratio: And if no prior knowledge about the absolute structure exists [i.e. pðy ¼ 0Þ To be able to do this, we define a value q analogous to z: This value represents z h for the inverted structure. Now Hence If the correct absolute structure and wrong absolute structure hypotheses are the only two possibilities for a certain structure, this would be sufficient. However, in practice a structure may be a twin consisting of two inverses (so-called inversion twins), and a more general probability model is desired to express this. The twinning can be described as a linear combination of the two structure factors of the pure enantiomeric structures. For each Bijvoet pair, This linear combination is analogous to the definition of the Flack x parameter (Flack, 1983). Since ÁF 2 c ðy ¼ 1Þ ¼ ÀÁF 2 c ðy ¼ 0Þ, this equation can be simplified to ÁF 2 c ðtwinÞ ¼ ð1 À 2xÞÁF 2 c ðy ¼ 0Þ ÁF 2 c ðy ¼ 0Þ: ð16Þ We refer to the variable as the 'generalized absolute structure'. For the correct absolute structure = 1.0 (and x = 0), and for the wrong absolute structure = À1.0 (and x = 1). With the help of this parameter, we can now introduce for each reflection h the function xðÞ: It can be seen easily that z h is equal to x h ð ¼ 1Þ and q h is equal to x h ð ¼ À1Þ. Note that this computation is also allowed with physically impossible values of jj > 1:0. With this generalization, the probability distribution becomes and from Bayes' theorem, We can now avoid the need to calculate pðobservationsÞ in two ways: we can either have a discrete set of possible hypotheses for the value of , or we can study the continuum of all possible values.
In the case of the discrete set 1 , 2 , . . . , n , we can normalize easily: In most cases, all priors pð i Þ will be set to 1=n. Two useful discrete sets of hypotheses that can be treated this way are the two-membered set the absolute structure is either correct or wrong and the three-membered set the absolute structure can be correct or wrong, but the sample may also be a 50/50 inversion twin.
If we want to consider the whole continuum of possible values of , the normalization of the probability function is less meaningful. In the case of a continuum, only ratios of different probabilities should be used, and these ratios do not depend on the normalization. However, all of the probability numbers are exceedingly small. For numerical stability reasons, it is advisable to bring all relative probabilities that we want to use in calculations to a reasonable size. To achieve this we can simply divide by a high value of the probability function. For this goal, we chose to use pðÞ with = 0 . We call the 'incompletely normalized' result p u : research papers 98 Hooft, Straver and Spek Absolute structure determination Since this approach is most useful if no prior knowledge is assumed at all (note that we always have prior information, namely À1 1, but here we explicitly choose to ignore this), we simplify it to It is observed in practice that p u ðÞ (in the second definition) is a reasonably well behaved Gaussian-like function. We can therefore calculate 2 a quantity G: Using this definition, G is the best approximation of for the structure based purely on the observations and not using any prior knowledge (not even the physical restriction that must be in the interval [À1, 1]). Since in our practical experience p u ðÞ looks very much like a symmetric Gaussian distribution, G will also be very close to the most probable value of . Like the value of the Rogers parameter, the value of G will be close to 1.0 for structures for which the absolute structure of the model is correct, and close to À1.0 for structures for which the absolute structure of the model is incorrect. Continuing along this path, we can calculate the variance of the distribution using This can be used to estimate an uncertainty in the obtained value of G.
The concept of the unrestricted absolute structure parameter G follows naturally from the comparison of the definitions of z and q. This is, however, a new concept. With a simple change of parameter expression we can cast our result in a form comparable with the Flack x parameter: and With this definition, y behaves like the Flack x parameter in that it will have a value of 0.0 for the correct absolute structure model, and 1.0 for the inverted model. Hooft, Straver and Spek Absolute structure determination 99 Table 1 Samples studied.

Test calculations
Conditions are as follows.
(1) Measured on a Bruker AXS SMART APEX system with an Mo K sealed-tube X-ray generator at room temperature.
(2) Measured on a Nonius KappaCCD system with a Nonius Mo K rotating-anode X-ray generator at 100 K.
(3) Measured on a Bruker AXS SMART 6000 system with a Cu K sealed-tube X-ray generator with graphite monochromator at 100 K. (4) Measured on a Bruker AXS SMART 6000 system with a Siemens Cu K rotatinganode tube and focusing multilayer optics at 100 K. (a) Data integrated using EvalCCD (Duisenberg et al., 2003) and scaled using SADABS (Sheldrick, 1996). (b) Data integrated with DENZO (Otwinowski, 1993) and scaled using SCALEPACK (Otwinowski, 1993). (c) Data integrated using SAINT (Bruker, 2004) and scaled using SADABS. R1 = P ðjF o j À jF c jÞ= P jF o j. AMBI is ammonium bitartrate, M048A is threonine, M049A is glutamic acid, M050A is ammonium bitartrate, M051A is alanine. The rest of the data were supplied to us by a pharmaceutical company.

Sample
Conditions interest at that time; others were specially collected to test the statistical methods introduced in this paper. All of the structures have weak resonant scattering signals. Roughly half of the data sets were collected using Mo K radiation. The theoretical resonant scattering signal at 2 = 0 was estimated for each of the data sets from Both summations run over atom types i, N i is the number of atoms of type i in the structure, f is the scattering factor of the atom type, and f 00 the imaginary part of the resonant scattering factor (Weiss et al., 2001). For Mo K radiation, f 00 ðOÞ = 0.0060, f 00 ðNÞ = 0.0033 and f 00 ðCÞ = 0.0016. For Cu K radiation, f 00 ðOÞ = 0.0322, f 00 ðNÞ = 0.018 and f 00 ðCÞ = 0.0091. ÁF=F is called the 'signal'. This does assume a random distribution of atoms in the cell; locations of resonant scatterers close to symmetry elements can cause weakening of the signal. On the other hand, this formula can be a pessimistic guess since f will decrease for increasing diffraction angles 2, whereas the resonant scattering factor f 00 is nearly independent of the diffraction angle. All structures were refined using SHELXL97. After refinement, the observed Bijvoet pairs listed in the FCF output file of SHELXL were used for an analysis of the value of y. Care was taken not to use FCF files produced by SHELXL run using the TWIN/BASF instructions, as in such a case the calculated structure factors already have the Flack x calculation embedded and this would invalidate the analysis. Where available, the absolute configuration assignment was crosschecked with prior information; in other cases the structure for which the value of y was closest to 0.0 was chosen. Results of the analyses are given in Table 2.

Results and discussion
The power of the introduced method comes from the fact that it weights each observed Bijvoet difference based on its expected accuracy directly, rather than relying on the weight of the reflection intensities. Calculating these proper weights for a least-squares procedure is very difficult, but proper weighting can be rather easily accomplished with the derived maximum-likelihood procedure instead of using least-squares. Bijvoet differences can be much smaller than the residual differences between the observed and calculated intensities, and the calculated differences are accurate as long as the resonant scatterers have been accurately positioned.
It is essential to measure Bijvoet pairs for the calculation of y, where the Flack x parameter can be determined even if the data set covers at least the asymmetric unit corresponding to the space group with an added inversion centre.
We only tested our methods on data sets with close to 100% coverage of Bijvoet pairs. When prior information is given, e.g. that the sample must be either the structure or its inverse, the method presented can be used to calculate probabilities of the two possible hypotheses [p2ðokÞ and p2ðwrongÞ]. These probabilities can be surprisingly decisive, even when the resonant scattering signal is very weak. For the test data sets measured using Cu K radiation, the chirality of all structures can be proven beyond reasonable doubt if it is assumed (prior knowledge) that the original compound was enantiopure. For the three-hypotheses model where the additional possibility of a 50% inversion twin cannot be ruled out, the distinction given by the probabilities [we call these p3ðokÞ, p3ðtwinÞ and p3ðwrongÞ] is less pronounced, but even in that case many of the determined values for the Cu K data sets would satisfy the most stringent pharmaceutical requirements. The least surprising results are obtained when the whole continuum of inversion twin structure compositions must be considered. In this case, the estimate that is obtained as the value of y has a smaller standard uncertainty than the value of the Flack x parameter. In most cases, the value of y is also closer to zero than the value of the  Table 2 Absolute structure analyses.
For some data sets, calculations were performed both on the correct and on the inverted model, refined in SHELXL. For the AMBI data set, the p3ðwrongÞ value of 1.64 Â 10 À10 increases to p3ðokÞ = 2.3 Â 10 À10 when the inverted structure is refined. The small difference between these values shows that the inverted refinement cannot absorb more than a small fraction of the resonant scattering signal into the other refined structural parameters. Comparing the equivalent numbers for the M006C data set, which has a weaker resonant scattering signal and which does not have a fixed origin, shows a similarly sized relative increase from p3ðwrongÞ = 0.031 to p3ðokÞ = 0.043 for the inverted structure. The magnitude of this structural bias is largely insignificant for the absolute structure determination of pharmaceutically active compounds. It may, however, be significant for accurate determination of the twin ratio of inversion twins; this has not been the subject of our study.
There are two assumptions in the derivation of the probabilistic model: firstly, that the standard uncertainty of the two reflections that form each Bijvoet pairs are independent; secondly, that the standard uncertainties of the individual reflections are accurate.
Both of these conditions are necessary conditions for xðÞ to follow a standard normal distribution. These assumptions can be verified by making a normal probability plot (Abrahams & Keve, 1971) from all values xð ¼ 1:0Þ. Such normal probability plots, made for the data sets above, show that the observed distribution of xð ¼ 1:0Þ for most data sets indeed follows a Gaussian distribution (the correlation coefficient of the normal probability plot is 0.999) but with < 1:0. Two possible reasons can be suggested. (i) The used scaling programs overestimate the errors in the reflection intensities. This is highly unlikely. (ii) The measurement error in the Bijvoet difference is smaller than could be expected if the two errors in the reflection intensities were independent. The errors are in fact positively correlated, and the error in the Bijvoet difference is really smaller. 3 The second hypothesis is most likely. Even without knowing the source of the smaller standard uncertainty, it is possible to use the information obtained from the normal probability plot to scale the standard uncertainties in the Bijvoet differences, thereby obtaining a corrected xðÞ. This correction scales down the standard uncertainties in in all but two of the cases that were examined for this paper.
The validity of such a downscaling of the errors can be confirmed by studying the result for a group of independent structure determinations and determining the value Z = y= y for each of them. If all standard uncertainties have been determined correctly, the values of Z from a random population of structure determinations should form a standard normal distribution. For the structure determinations using Mo K radiation given in this paper, the average absolute value of Z is 0.61 (the expected value is 0.85) and the root mean square (r.m.s.) value of Z is 0.67 (expected 1.0) ( Table 3). These results suggest that the error is indeed systematically overestimated. After applying the slope from the normal probability plot to correct the estimated standard uncertainties in the observed Bijvoet differences, the average absolute value of Z is 0.66 and the r.m.s. value of Z is 0.72. These values are still smaller than the expected values. The current benchmark set is too small for this to be considered proof of the merits of the downscaling procedure.  and  investigated the value of the Flack x parameter for a set of centrosymmetric structures that were refined in a noncentrosymmetric space group. Looking at the definition of the Flack x parameter,

Centrosymmetric structures
it can be clearly seen that for the correct model, x is indeterminate since the two terms F calc h and F calc Àh are equal. The determination of x in these cases is therefore based purely on the random incorrect differences between the two 'halfstructures' in the refinement. In this light, it is at first sight surprising and discomforting that the values observed have such small standard uncertainties. It is clear that for the Flack x parameter the assumption that the off-diagonal elements of the covariance matrix may be ignored is wrong. The assumption that all other parameters have been determined correctly by the least-squares refinement has been violated.
We have attempted a non-centrosymmetric solution and refinement of a centrosymmetric ruthenium-containing compound (Hotze et al., 2005) Table 3 Correction of the calculated value of y for the error in the standard uncertainties as derived from a normal probability plot.
Z is the deviation of y from the enantiopure value expressed in units of y .  (3). Both values are close to 0.5 with a relatively small standard uncertainty. A detailed analysis of the data set indicated that the small standard uncertainty is due to a few reflections for which the differences between the two half-structures create a significant Bijvoet difference ÁF c while, as expected for a centrosymmetric structure, the ÁF o value is statistically insignificant. Such pairs are normally indicative of twinning by inversion. The only statistical difference in reciprocal space between a real inversion twin and a wrongly refined centrosymmetrical structure is that the calculated Bijvoet differences are much smaller than for a normal non-centrosymmetric structure with the same elemental composition. This is due to the fact that the configuration of the atoms is almost centrosymmetric (with respect to a suitably chosen origin, the phases of many reflections are close to 0 and and the phases of the resonant scattering contributions are close to =2 and 3=2) and hence the resonant scattering contribution to the scattering factors only results in relatively small scattering amplitude differences. It is difficult to determine a reliable criterion for this effect.
It appears then that the distinction between a true inversion twin and a non-centrosymmetrically refined centrosymmetric structure is best made in real space by a symmetry-detection procedure like ADDSYM (Spek, 2003), followed by a detailed inspection of the weak reflections after refinement in the suggested centrosymmetric space group.

Recommended procedure
Current versions of refinement programs cannot use the value of y to take absolute structure into account. We therefore recommend to refine the structure including the Flack x parameter (e.g. use the TWIN/BASF instructions in SHELXL). The value of y can then be determined separately using a utility that explicitly calculates structure factors for the Bijvoet pairs (e.g. the Bijvoet Pairs option in PLATON). 4 This procedure will account for any correlation between the structural parameters and the absolute structure.

Conclusions
A new probabilistic procedure was introduced that can be used to establish the absolute structure. The procedure is especially suitable for biologically active compounds, which often do not contain atoms with a larger resonant scattering signal than that of oxygen.
The only special requirement for the data collection procedure imposed by the new probabilistic calculation is that Bijvoet pairs should be present in the data set. In contrast, the determination of the Flack x parameter also works for data sets that have a Bijvoet pair coverage of 0%, although this is not recommended practice.
One of the results of the procedure is a value y, which can be directly compared with the value of the Flack x parameter. We observe for our test data sets that the standard uncertainty in the value of y is roughly half of the standard uncertainty in the value of x. The observed deviation from 0.0 is consistent with the standard uncertainty. These observations are comparable with the results obtained using invarioms but without the significant efforts associated with the calculation of invarioms.
The calculations also give explicit probabilities for the absolute structure assignment, without referring to the value of y and without the assumption that the distribution of y is Gaussian. The explicit probability of an absolute structure assignment error makes our procedure suitable to regulate the probability of erroneous assignments in pharmaceuticals. The probability calculations can be based either on a model with two hypotheses for the two absolute structures or optionally take the chance of a racemate into account as a third hypothesis.
The procedure was tested on a number of light-atom structures (no atoms with a stronger resonant scattering signal than that of oxygen). For those data sets collected using Mo K radiation, a mixed result was obtained: some structures could receive quite a good absolute structure assignment; most structures show at least a clear direction. For all data sets measured using Cu K radiation (resulting in a roughly five times larger resonant scattering signal) an excellent absolute structure discrimination was obtained with the chance of error for most structures below 10 À100 . Of course, for most of the Cu K structures the standard uncertainty in the Flack x parameter is small enough for an unambiguous assignment.
The current method offers an alternative method to look at the same experimental data as addressed by the Flack x approach.