Method for Modelling the Electron Density of a Crystal

This specification includes microfiche Appendix having 3 microfiche sheets with 255 frames. Additionally, there is a hard copy of the microfiche appendix (not printed in the patent) in the application file.

A portion of the disclosure of this patent document contains material which is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure as it appears in the Patent and Trademark Office patent file or records, but otherwise reserves all copyright rights whatsoever.

The present invention relates to crystallographic methods and apparatus for determining the three-dimensional structure of macro-molecules by crystallography or electron microscopy.

Under special conditions, molecules condense from solution into a highly-ordered, crystalline lattice, which is defined by a unit cell, the smallest repeating volume of the crystalline array. The contents of such a cell can interact with and diffract certain electromagnetic and particle waves (e.g., X-rays, neutron beams, electron beams etc.). Due to the symmetry of the lattice, the diffracted waves interact to create a diffraction pattern. By measuring the diffraction pattern, crystallographers attempt to reconstruct the three dimensional structure of the atoms in the crystal.

A crystal lattice is defined by the symmetry of its unit cell and any structural motifs the unit cell contains. For example, there are 230 possible symmetry groups for an arbitrary crystal lattice, and each symmetry group may have an arbitrary dimension that depends on the molecules making up the lattice. Biological macromolecules, however, have asymmetric centers and are limited to 65 of the 230 symmetry groups. See Cantor et al., Biophysical Chemistry, Vol. III, W. H. Freeman & Company (1980), which is incorporated herein by reference for all purposes.

A crystal lattice interacts with electromagnetic or particle waves such as X-ray or electron beams, respectively, that have a wavelength with the same order of magnitude as the spacing between atoms in the unit cell. The diffracted waves are measured as an array of spots on a detection surface positioned adjacent the crystal. Each spot has a three-dimensional position, hkl, and an intensity, I(hkl), that both are used to reconstruct the three-dimensional electron density of the crystal with the so-called Electron Density Equation: ##EQU1## where .rho.(xyz) is the electron density at the position (xyz) in the unit cell of the crystal, V is the volume of the unit cell, and F(hkl) is the structure factor of detected spot located at point (hkl) on the detector surface. As expressed above, the Electron Density Equation states that the three-dimensional electron density of unit cell is the Fourier transform of the structure factors. Thus in theory, if the structure factors are known for a sufficient number of spots in the detection space, then the three-dimensional electron density of the unit cell could be calculated using the Electron Density Equation.

A number of problems exist, however, in actual parctice. The Electron Density Equation requires knowledge of the structure factors, F(hkl), which are generally complex numbers that consist of both an amplitude and a phase. The amplitude of a structure factor, .vertline.F(hkl).vertline., is simply the square root of the experimentally measured intensity, I(hkl). The phase of each structure factor, on the other hand, is not known and cannot be measured directly in a diffraction experiment. Nor can it be derived directly for macromolecules. Without the phase of each structure factor, determination of the three-dimensional structure of most large structures by the use of the Electron Density Equation is impossible, except for special cases.

Theoretical methods are exemplified by the Direct Method and the Patterson Method or their extensions, such as the maximum entropy method or use of simulated annealing in both reciprocal and Patterson space. These methods calculate the phases directly from the measured intensities of the diffracted waves and allow routine computer solutions for molecules having typically less than approximately 100 non-hydrogen atoms. (As is known in the art of crystallography, hydrogen atoms contribute little to the diffraction process.) For structures having more than 100 non-hydrogen atoms, such as proteins, peptides, DNA, RNA, virus particles, etc., such direct methods become impractical and, in most cases, impossible. Fortunately, experimental methods, such as Multiple Isomorphous Replacement and Anomalous Scattering, exist to aid in the determination of these phases.

Multiple Isomorphous Replacement is based on the observation that the absolute position and, therefore, the phase of the structure factor of a heavy atom incorporated into an otherwise unmodified crystal lattice can be determined. With this knowledge, the phase of each structure factor in the derivative is determined relative to that of the heavy atom. Except for crystals having centrosymmetric symmetry, at least two heavy metal derivatives are required to unambiguously determine the phase of a structure factor. Furthermore, Multiple Isomorphous Replacement requires that each heavy metal derivative not otherwise change the structure of the molecule or distort the unit cell of the crystal.

Other experimental techniques, used in conjunction with Multiple Isomorphous Replacement allow the crystallographer to forego analysis of some heavy metal derivatives. One such technique, Anomalous Scattering, is based on the observation that particular heavy atoms scatter radiation of different wavelengths significantly differently. With this technique, one heavy metal derivative studied at two wavelengths yields data equivalent to two heavy atom derivatives studied at one wavelength.

Other techniques completely circumvent the preparation and study of heavy metal derivatives.

Molecular replacement, as the name suggests, uses a molecule having a known structure as a starting point to model the structure of the unknown crystalline sample. This technique is based on the principle that two molecules that have similar structures and the same orientation and position in the unit cell diffract similarly. Effective use of this technique requires that the structures of the known and unknown molecules be highly homologous.

Molecular replacement involves the complicated task of positioning the known structure in the unit cell in same location and orientation as the unknown structure. Difficulty in using this technique arises because the result is critically dependent on the exact position of the known structure. Slight variations in either the location or orientation of the known structure often results in complete failure. Once positioned, the atoms of the known structure in the unit cell are used in the so-called Structure Factor Equation to calculate the structure factors that would result from a hypothetical diffraction experiment. The Structure Factor Equation takes the form: ##EQU2## where F(hkl) is the structure factor of the molecule at the point (hkl) on the defector surface, f.sub.j is the atomic structure factor (that is, it represents the scattering properties of the individual atom), N is the number of non-hydrogen atoms, and x.sub.j, y.sub.j, z.sub.j are the factional coordinates of atom j in the unit cell. The structure factor calculated is generally a complex number containing both the amplitude and phase data for the molecular replacement model at each point (hkl) on the detector surface. These calculated phases are used, in turn, with the experimental amplitudes measured for the unknown structure to calculate an approximate electron distribution. By refinement techniques, this approximate structure can be fine-tuned to yield a more accurate and often higher resolution structure.

The molecular replacement technique requires knowledge of the number of molecules in the unit cell, and the orientation and position of each molecule within the unit cell. Initially the electron density calculated from the phases from molecular replacement model and experimental amplitudes closely resembles the electron density of the model. Only after refinement of the initial structure will the success or failure of the method be apparent. For instance, failure occurs if the initial structure fails to converge (as represented by a correlation value) or if the refined structure diverges from the structure of the model during the refinement process. In cases where the unknown structure is a substrate or an intermediate bound to an protein, molecular replacement's success is evident when the result is a structure whose only difference is added electron density that represents the protein-bound molecule. The determination of such structures is important in the area of pharmaceutical drug testing where the structure of protein-bound drugs and intermediates yield important information about binding and mechanism. Similarly, new mutants of a protein or variations of protein-bound inhibitors are well suited for molecular replacement, as are structures of the same molecule that have crystallized in different symmetry groups.

Molecular Replacement is not always effective, however. Determination of the number of copies of the model in the asymmetric unit and the correct location and orientation of each copy is critical and time consuming, since ideally one samples all rotational and translational degrees of freedom in the asymmetric unit to determine the correct set of parameters.

Multiple Isomorphous Replacement, Molecular Replacement, and their related techniques, do not work for all cases, however, and there exists needs for simplified, efficient methods to determine the structure of crystalline molecules. The present invention fulfills these and other needs.

The present invention produces a model of the electron density distribution of a macromolecule in a defined asymmetric unit of a crystal lattice in a multi-step methods. These methods are simple and rapid ways of modelling the electron density of a crystal, without the need to determine the phases of the reflection data. According to one aspect of the invention, data collected from an X-ray diffraction experiment of a crystal lattice are inputted into a computer and are converted into normalized amplitudes. The crystal lattice has a defined asymmetric unit having locations of uniformly diffracting electron density. An initial distribution of scattering bodies is produced within a asymmetric unit that has the same dimensions as the defined asymmetric unit of the crystal lattice. The scattering amplitudes and phases of the initial distribution are then calculated, and the correlation between the calculated scattering amplitudes and the normalized amplitudes is calculated to determine the fit between the two data sets. At least one of the scattering bodies within the unit cell is moved to create a modified distribution, the scattering amplitudes of this modified distribution is calculated, and the correlation between the calculated amplitudes and the normalized values is recalculated. A final distribution of scattering bodies is produced by repeating the steps of moving at least one of the scattering bodies to create a modified distribution and determining the correlation between the calculated amplitudes and the normalized values, until the correlation between the calculated scattering amplitudes and the normalized amplitudes is effectively maximized. This final distribution of scattering bodies defines the electron density of the crystal.

In a preferred embodiment of the present invention, a scattering body is moved a predetermined distance. In another preferred embodiment, the method includes additional steps. The final distribution is refined by at least one refining step that reduces the predetermined distance that a scattering body is moved. At least one of the scattering bodies is moved by the reduced distance within the asymmetric unit to modify the distribution, the scattering amplitudes of this modified distribution are calculated, and the correlation between the calculated amplitudes and the normalized amplitudes is determined. Finally, the refining step produces a final distribution of scattering bodies by repeating steps of moving at least one scattering body and calculating the correlation between the calculated amplitudes of the distribution and the normalized values, until the correlation is effectively maximized. The final distribution of scattering bodies define a refined electron density of the crystal.

Other preferred embodiments include one or more of the following features. The refining step is repeated with decreasing move distances, until the move distance is reduced to a predetermined final value. The scattering bodies are translated in a random translation direction, and the random translation direction is randomly selected from predefined translation directions that are parallel to the axes defined by the crystal lattice unit cell.

In another preferred embodiment, the steps of determining the fit between the calculated amplitudes of scattering bodies and the normalized amplitude consist of calculating a correlation coefficient between the calculated amplitudes and the normalized values.

In another aspect of the invention, a model of the electron density distribution of a macromolecule in a defined asymmetric unit of a crystal lattice having locations of uniformly diffracting electron density is produced by inputing data collected from an X-ray diffraction experiment into the computer. A plurality of scattering bodies is randomly distributed in another asymmetric unit having substantially the same dimensions as the defined asymmetric unit. The plurality of scattering bodies is then moved into a final distribution, whereby the scattering amplitudes of the final distribution have a maximum fit with amplitudes from the data. This final distribution defines the electron density distribution of the macromolecule in the defined unit cell.

FIG. 1 a is a schematic drawing of an X-ray diffractometer;

FIG. 1b is a schematic representation of the detection plate of FIG. 1;

FIG. 1c is a block diagram illustrating computer hardware to which the invention may be applied;

FIG. 2 is a schematic flowchart showing the steps in a microcycle;

FIG. 3 is a stereoscopic view of an initial random distribution of scattering bodies;

FIG. 4 is a stereoscopic view of the final distribution of scattering bodies after the application of the condensing protocol of the present invention;

FIG. 5 is a stereoscopic view the alpha carbon backbone of Elastase superimposed on the final distribution of scattering bodies of FIG. 4;

FIG. 6a and 6b is a graph showing the behavior of the Pearson correlation coefficient, r, during the condensing protocols;

FIG. 7 is a graph showing the behavior of the spatial distribution, j, during the condensing protocol;

FIG. 8 is a stereoscopic view of an initial distribution of scattering bodies used in modeling the electron density of the Rlt69 the bacteriophage 434 repressor protein

FIG. 9 is a stereoscopic view of the final distribution of scattering bodies after application of the condensing protocol;

FIG. 10 is a stereoscopic view of the alpha carbon backbone of the protein superimposed on the final distribution of scattering bodies of FIG. 9;

FIG. 11 is a graph showing the behavior of the Pearson correlation coefficient, r, during the condensing protocol; and

FIG. 12 is a graph showing behavior of the spatial distribution, j, of during the condensing protocol of the present invention.

The present invention will be described by providing details of each step involved in the modelling.

When macromolecules crystallize, solvent typically occupies large portions of the crystal lattice including space on the interior as well as space exterior to the macromolecular envelope. The bulk of the incorporated solvent is fluid and consists of randomly oriented molecules that do not diffract X-rays, electron beams, neutron beams, and the like. The portion of the unit cell occupied with the macromolecule is termed the "positive image," while the portion occupied with solvent is termed the "negative image." The term "macromolecular envelope" denotes the surface of the macromolecule and may be determined at different resolutions. Thus, a high resolution macromolecular envelope includes many detailed features of the macromolecular surface, such as location of side chains, clefts, etc. Conversely, a low resolution envelope includes few detailed features and provides details about the general shape of the macromolecule.

The term "macromolecule" includes, but is not limited to, the following: biological macromolecules such as proteins, peptides, RNA, DNA, complexes of peptides and nucleic acids, virus particles, organelles, and the like; organic molecules such as organic polymers, plastics; inorganic molecules such as zeolites; and other large molecular structures. Although the term "protein" is used below in conjunction with the description of illustrative embodiments, the method is fully suited for structure determination of other macromolecules that crystallize into a unit cell having solvent space.

In a preferred embodiment, an initial distribution of scattering bodies is rearranged into the final distribution that represents the negative image of the crystal. That is, the scattering bodies in the final distribution are located in the portion of the crystal occupied with solvent rather than the portion occupied with the macromolecule. The electron density distribution of the macromolecule is modelled as the portion of unit cell that is substantially devoid of scattering bodies.

The invention methods include different steps: collecting diffraction data; inputing experimental crystallographic data into a computer which is used to perform the method; distributing scattering bodies in a corresponding asymmetric unit; calculating scattering data from this distribution; determining the correlation between the experimental and calculated scattering data; moving at least one of the scatterers; calculating the new scattering data from this new distribution; determining the correlation between the new scattering data and the experimental data; and producing a final distribution by repeating the moving, calculating and correlating steps until the correlation between the experimental and calculated data is effectively maximized.

Collection of diffraction data from scattered waves is well known in the art of crystallography. Referring to FIG. 1, a diffractometer 10 (also known as an X-ray Set) for use with the present invention includes a source of X-rays 12, a sample holder 14, and a detection apparatus 16. X-ray source 12 produces a collimated beam 18 of X-rays having a relatively narrow cross section and is typically a mercury flash tube, copper cathode or rotating anode that produces X-rays having a narrow and well-defined wavelength spectrum.

In other preferred embodiments, alternate forms of radiation and radiation sources are used. For example, an alternate source of X-rays, such as a tunable X-ray source (e.g., radiation from a synchrotron or other source that emits X-rays of different wavelengths) is preferred for use with techniques such as Anomalous Scattering or Multiple Wavelength Scattering. Alternate forms radiation include electron beams (e.g., those typically used in electron microscopy), neutron beams (e.g., those typically used in neutron beam diffraction), and the like.

Sample holder 14 consists of a capillary tube 20 having a crystallized sample 22 located within its lumen. Capillary tube 20 and crystallized sample 22 are positioned in the path of collimated X-ray beam 18. X-rays 24 diffracted by the crystal impinge a detection apparatus 16 that is positioned generally opposite X-ray source 12 and that consists of a detector surface 26 mounted on an arm 28. Detection surface 26 takes on many shapes, such as a two-dimensional disk, a three-dimensional cylindrical surface, and the like, and is adapted to record the position and intensity of diffracted X-rays 24. Examples of suitable detection surfaces for use in the present invention include photographic film, gas chamber or multiwire area detectors, CCD channel plates, image plates, diffractometers including precession cameras, and the like. In many apparatus, e.g., a precession camera, the rotation of sample holder 14 and detector surface 26 are coupled during the diffraction experiment to accumulate data corresponding to an entire plane in reciprocal space.

Referring now to FIG. 1a, detector surface 26 after exposure to diffracted X-rays consists of an array of spots that each have a position, (hkl), and an intensity, I(hkl). The data form an array having a circular boundary 30 that represents the high resolution limit of the diffraction experiment. Data representing low resolution features of the crystal are located near the center of the circular array while data representing the high resolution features are near the outer edge. For example, data lying within the circle defined by line 28 represent lower resolution features of the crystal (e.g., down to 10 .ANG.), while data lying between circles 29 and 30 represent features of higher resolution (e.g., between 2 .ANG. and 5 .ANG.).

During a diffraction experiment, data collection typically involves accumulation of a large number of data points, often over 10,000. After accumulation on an appropriate detection surface, the position and intensity of each data point is measured, as is known in the art of crystallography, and the data are put into a computer for storage and further processing. In a preferred embodiment, the computer is a digital computer, such as a VAX 8550, produced by Digital Equipment Corporation of Maynard, Mass. Other computers of varying computational power are also suitable: supercomputers, multiprocessor computers, mainframe computers, work stations, personal computers, and the like. Exemplary computers for use with the present invention include computers produced by Cray Research, Digital Equipment Corporation, Thinking Machines, Data General, International Business Machines, Apple Computer, Sun Computers, and Silicon Graphics. In other embodiments special duty computers are used. For example, an appropriate digital computer incorporated into a diffractometer is suitable for performing the invention method.

Referring now to FIG. 1c, a computer 32 used in conjunction with the present invention include an interface 33 to receive data from a diffraction apparatus 34, memory 35 (e.g., RAM), file storage 36 (e.g., magnetic disk or tape) to store the data, and a CPU 37 to process the data. In preferred embodiments, the computer further includes an output device 38, such as a printer, plotter, or graphics display that allows the resulting electron density of the crystal to be displayed graphically. Typical graphic displays are produced by Evans & Sutherland and Silicon Graphics.

As previously mentioned, the crystallographic symmetry and dimensions of the unit cell and asymmetric unit are determined directly from the data (See Blundell and Johnson "Protein Crystallography" Academic Press, NY 1976, which is incorporated herein by reference for all purposes). The unit cell is the smallest portion of the crystal lattice that repeats upon operation of a translation. Thus, the crystal is composed of a repeating array of unit cells in three dimensions, and determining the electron density of the unit cell is equivalent to determining the electron density of the crystal lattice. In most space-groups, the unit cell has multiple copies of the crystallized macromolecule, and may have internal symmetries, such as a plane of reflection, an n-fold rotation axis, etc., termed "crystallographic" symmetries. Subspaces of the unit cell related by such crystallographic symmetries, are termed "asymmetric units." Determining the electron density of the asymmetric unit is also equivalent to determining the electron density of the crystal lattice, by first applying the appropriate symmetry operations to reconstruct the unit cell and then translating the unit cell. As used herein, the term "asymmetric unit" refers to portion of the crystal lattice that can be repeated by translations, rotations, or combinations of thereof, to reconstruct the crystal lattice. Thus, for some crystals, the asymmetric unit is the same as the unit cell (when there is no known crystallographic symmetry). In other crystals for example, the asymmetric unit may be half of the unit cell. Once the symmetry and dimensions of asymmetric unit are determined, as is known in the art, the experimental data are converted from unscaled F-values into a form convenient for use in the method of the present invention.

In a preferred embodiment, a portion of the experimental data are converted into normalized structure-factor magnitudes (i.e., E-values) that are conventionally used in Direct Methods (See Karle, Acta Crystal. 1989, vol. A45, pp. 765-781, which is hereby incorporated by reference for all purposes). As used herein, "portion" is used to indicate a subset or the whole of the experimental data, since some cases require conversion of the entire data set, while others require conversion of a subset. While F-values represent scattering by atoms that have a finite electron distribution, normalized E-values represent scattering by bodies that have no spatial distribution and have a simple scattering cross-section. Thus, this conversion models the crystal as an array of point scatterers rather than atoms. Alternatively, the accumulated data is converted into properly scaled F-values, as described in Blundell and Johnson. The use of properly scaled F-values or normalized E-values depends on the details of the calculations. In a preferred embodiment, E-values are used in conjunction with high resolution data. With data having resolutions lower than approximately 8 .ANG., however, either F-values or normalized E-values are used (with a corresponding change in the Electron Density and Structure Factor Equation, as is well known in the art) with no substantial difference in the results. In ensuing discussions, therefore, properly scaled F-values and normalized E-values are used interchangeably at lower resolutions unless otherwise specified, and are collectively referred to as "experimental data."

Once diffraction data are collected and converted to properly scaled F-values or normalized E-values, and the crystallographic symmetry and dimensions of the asymmetric unit have been determined, the electron density of the crystal is modelled by a distribution of point-scatterers.

To model the electron density of the crystal, an initial distribution of scattering bodies is created and allowed to condense (that is, to rearrange) into a final distribution, but two constraints. The first constraint is based on a "physical" packing of the scattering bodies, and the second is based on correlating the scattering resulting from the distribution of scattering bodies to the experimental data. The initial distribution of scattering bodies is created by placing a plurality of scatterers into an asymmetric unit having the same symmetry and dimensions as the asymmetric unit of the crystal lattice. The scattering bodies are hypothetical objects used to model the electron density of the crystal and have physical characteristics such as a radius and scattering cross-section. As described below, the number and properties of the these scatterers depend on the properties of the asymmetric unit of the crystal.

Like physical objects, no part of more than one scattering body can occupy the same space in the asymmetric unit; in other words, two scattering bodies can approach each other until their surfaces touch. Without this physical limitation, the scattering bodies may condense into the same small region in the asymmetric unit. The scattering bodies have an outer surface, and can be of many shapes and sizes. In a preferred embodiment, a scattering body is a sphere, ellipsoid, cube, tetrahedron, etc. In more preferred embodiment, the scattering bodies are spherical in shape and have a predefined radius, r, although the radius of the scattering body may vary, as described more fully below. In a most preferred embodiment, each scatterer has a radius of 1.5.

Although each scattering body is preferably spherical and has a radius, each is treated as a point scatterer having a scattering factor of unity. Such a choice is made for the convenience of computation, since application of the Structure Factor Equation is easier to compute for point scatterers. In other preferred embodiments, however, the scatterers have scattering profiles that approximate the spatial distribution of a real atom, such as a Gaussian profile or a Normal profile. The Structure Factor Equation becomes more complex, however, and the computation time increases upon increasing the complexity of these scattering profiles.

Point scattering bodies cannot provide a phase description of the macromolecule to a resolution better than their inter-sphere collision distance (e.g., 3 .ANG. of spheres having a 1.5 .ANG. radii). Application of the present invention to solving structures at resolutions higher than 3 .ANG. requires scatterers having smaller radii, as discussed below.

The number of scatterers distributed within the asymmetric unit depends on many factors, such as the radii of the scatterers, the number of non-hydrogen atoms in the asymmetric unit, the solvent fraction of the asymmetric unit (or, equivalently, the expected packing fraction of the macromolecule), and the resolution and number of reflections experimentally collected.

The number and physical characteristics of the scattering bodies that optimally satisfy the expected protein packing fraction are chosen. Optimal packing of the scatterers allows the distribution of scatterers to have mobility within the solvent portion of the asymmetric unit during the condensation process. That is, the scatterers have enough free space in which to move. The number of spherical scattering bodies of a radius, r, required to maximally fill the volume occupied by the solvent incorporated into the asymmetric unit in a particular crystal, N.sub.max, is determined. Once this maximum number is known, allowance is made for mobility within the distribution of scatterers.

Biological macromolecules generally crystallize with a large amount of solvent incorporated into the asymmetric unit, and it is possible to estimate the percentage of solvent content with reasonable accuracy, as is known in the art. In the case of proteins, enzymes, polymeric nucleic acids etc., the primary sequence may be used to estimate the volume of the protein. Typical algorithms for calculating the volume of amino acids and polypeptides are well known in the art. In a preferred embodiment, an average value the solvent fraction of about 0.4 (i.e., 40% of the unit cell volume) is used as an estimate for the solvent fraction of the asymmetric unit.

For a plurality of spheres having the same radius, the theoretical value for the maximal random packing fraction is about 0.6. That is, spheres randomly packed to maximize their density will fill 0.6 (i.e., 60% of the volume). This value, combined with the average solvent content of the asymmetric unit (i.e., 40%), provides the number of spherical bodies needed to maximally fill the solvent space as:

where r is the radius of each sphere. This value is unrealistic for the invention method, however, since it requires that the spheres lie in an extremely tightly, packed lattice. Even if there were no Fourier constraints to be satisfied, the requirement that the distribution of scattering bodies be mobile and fluid requires that the optimum number of scattering bodies be significantly smaller than the maximum number. Additionally, there are limitations imposed by the Fourier constraints that further decrease the number of scatterers, as described below.

In a preferred embodiment, the optimal number of scattering bodies, N.sub.hs, is between (0.01) V.sub.asymmetric unit /V.sub.sphere and (0.15) V.sub.asymmetric unit /V.sub.sphere. In a more preferred embodiment, the number of scattering bodies is between (0.04) V.sub.asymmetric unit /V.sub.sphere and (0.08) V.sub.asymmetric unit /V.sub.sphere.

In the case of protein or peptide crystallography, when the packing fraction of the peptide is in the range P=0.5 to 0.6, the number of alpha carbons (C.sub..alpha.) of a typical protein, N.sub.o.alpha., is approximately 1.5 times greater than the value for N.sub.hs. Thus, in a preferred embodiment, N.sub.hs =0.67 N.sub.o.alpha..

Once the number of scatterers is determined, as described above, the asymmetric unit is randomly filled with the appropriate number (e.g., N.sub.hs =0.67 N.sub.o.alpha. for typical proteins) of scattering bodies of a predetermined radius, for example 1.5 .ANG.. As described above, the centers of any two 1.5 .ANG. radius spheres cannot approach closer than 3 .ANG..

In a preferred embodiment, the initial distribution of the spheres is produced by positioning a first scattering body within a volume having the same dimensions and symmetry as the asymmetric unit of the crystal. Subsequently, a second scattering body is positioned in the asymmetric unit under the constraint that its center is at least 3 .ANG. away from the center of the other body. Placement of scattering bodies continues until the asymmetric unit is filled with the desired number of scattering bodies.

Many methods are suitable for the positioning of each scattering body. In a preferred embodiment, each scattering body is placed at a random position, subject to the constraint that the centers of no two spheres are closer than 3 .ANG.. In another preferred embodiment, each scattering body is placed at a regular position in the asymmetric unit, such as a grid point of a rectilinear array, subject to the same constraint as described above.

Before further steps are performed, the unit cell is reconstructed from the asymmetric unit, if necessary by the appropriate symmetry operations, as described above.

Once the initial distribution of the scattering bodies is determined, the Fourier amplitudes of this distribution are calculated by a trigonometric summation using the Structure Factor Equation. In a preferred embodiment, the following equation is used: ##EQU3##

In other preferred embodiments, other methods for calculating the Fourier amplitudes are used. Suitable methods include, for example, Fast Fourier Transfer methods (See Press et al. Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 1988, incorporated herein by reference for all purposes.)

The calculated Fourier amplitudes and the experimental data are then correlated to determine the fit between the positions of the scattering bodies in the initial distribution and the positions of the atoms in the crystal. Many methods for determining the correlation between two sets of data exist. In a preferred embodiment, the Pearson correlation coefficient, r, is used. The Pearson coefficient takes the form:

where E is the experimental data, E.sub.c is calculated amplitude, and the summations are taken over all experimental data points. To maximize the fit between data, the value of the Pearson coefficient is maximized.

In another preferred embodiment the crystallographic correlation, R, is used to correlate the data sets (See Blundell and Johnson). This correlation exhibits behavior opposite of the Pearson Coefficient: a closer correlation between the experimental data and the calculated amplitudes results in a smaller value of R. Thus, to maximize- the fit between the two data sets, R is minimized. Other methods of correlating the two data sets are well known and will be apparent to those skilled in the art.

The choice of resolution, K, at which the two data sets are correlated depends primarily on two conditions. First, the resolution cannot exceed the inter-collision distance of the scattering bodies 3 .ANG.. Second, there should be sufficient over-determinacy; that is, the number of experimental data points should be larger than the number of scattering bodies. One can use more of the experimental data, but for a low resolution image such data makes little difference. For instance, when the crystallized molecule has a large unit cell, data having a resolution of approximately 7 to 10 .ANG. generally provides sufficient over-determinacy. For smaller unit cells, however, sufficient over-determinacy is attained by using higher resolution data, for instance down to approximately 4 or 5 .ANG. resolution.

After the parameter K has been chosen, the correlation coefficient between the Fourier amplitudes is calculated from the initial distribution of scattering bodies and the experimental data.

After the initial correlation has been determined, the distribution of scattering bodies is modified to increase the correlation between the calculated amplitudes and the experimental data.

The general strategy used to modify the distribution is to randomly select one of the N.sub.hs scattering bodies and randomly move it. The distance and direction of the movement are constrained, however. In a preferred embodiment, the scattering body is moved an initial predetermined distance in a random direction. The predetermined distance is chosen according to the other parameters, such as the asymmetric unit dimensions, the correlation coefficient, number of scattering bodies, etc., and is chosen to allow the scattering bodies to "explore" the asymmetric unit. In other embodiments, the scattering bodies each are moved a random distance within a predetermined range. For example, a scattering body may be moved a random distance between zero and one-third the average dimension of the asymmetric unit. In another preferred embodiment, each scattering body is constrained to move parallel to one of the six directions defined by the unit cell edges.

The random movements of the scattering bodies are constrained: some movements are not allowed. For instance, after a scattering body is moved, it cannot occupy the same space as another scattering body (i.e., the surfaces of any two scattering bodies cannot intersect). If this constraint is violated, the move is rejected and another scattering body is moved randomly, subject to the same constraint. If the scattering body is moved to a location that is outside of the asymmetric units it is repositioned back into the unit by use of the appropriate space-group dependent symmetry operator (the scattering body is "folded" back into the asymmetric unit). Upon completion of an allowable move, the scattering amplitudes for the new distribution is calculated, and the correlation coefficient is reassessed for this new distribution of scattering bodies. If the correlation coefficient is more favorable (indicating a closer fit), then the move is accepted. Otherwise the move is rejected, and the sphere is returned to its original position. In this way only moves that result in a closer fit are allowed.

This process of moving a scattering body is defined as a "microcycle." A microcycle is "attempted" if the movement of the body is allowed (that is, the "physical" constraint is satisfied) and the correlation coefficient is calculated. Otherwise the microcycle is "rejected." If the correlation coefficient calculated for the new distribution indicates a closer fit between the experimental data and calculated amplitudes, the move is "accepted," and a new microcycle is started. If, however, the correlation coefficient indicates a worse fit, the move is "rejected."

Referring now to FIG. 2, the steps involved in each microcycle are shown schematically in a flowchart. After the initial distribution of scatterers is selected 40 and the correlation coefficient, r, is calculated 42, one of the scatterers is randomly selected and moved 44 under the constraints described above. If the new position is out the unit cell 46, then the appropriate symmetry operator is applied 47 and the scatter moved back into the asymmetric unit. The new position of the scatterer is determined and if it is further than 3 .ANG. away from another scatterer (steps 46 and 48, respectively), the correlation coefficient is calculated for this new distribution 50. If the movement criterion 48 is violated or r increases 51 (indicating a worse fit), the scattering body is returned to its original position 52, and the microcycle begins again.

After the first microcycle is accepted, another scattering body is randomly selected and moved randomly according to the previously-described constraints. In a preferred embodiment, the scattering bodies are moved by the same predetermined distance, or step-size, x, during each microcycle until the distribution has condensed to a stable state (i.e., when the correlation coefficient has converged to a stable environment that indicates an effective maximum between the correlation between the calculated amplitudes of the scattering bodies and the experimental data). The number of microcycles that are attempted as well as the number of microcycles that result are accepted and rejected are tabulated throughout the condensation process. When 10 accepted moves occur before a total of 100 attempted microcycles have occurred, the collection of attempted moves is collectively defined as a "condensing microcycle," and indicates that the distribution of scattering bodies, as a whole, is converging to a closer fit with the experimental data. If 10 moves have not been accepted before 100 attempted moves, the set of 100 attempted moves is defined as a "condensed macrocycle," which suggests that the distribution us in a stable environment. At this stage, the distribution of the scattering bodies is effectively maximized at the current step-size. That is, further movement of scatterers will probably not increase the correlation between the data. A set of approximately 200 consecutive macrocycles together constitute a "supercycle." In a preferred all embodiment microcycles within a supercycle have the same step-size (i.e., each scattering body is moved the same distance within a supercycle). A supercycle can have less than 200 macrocycles if 40 condensed macrocycles occur consecutively, which indicates that the distribution has converged to a close fit with the experimental data using the current step-size that the scattering bodies are moved.

Before the next supercycle is started, the step-size, x, is decreased. Thus, consecutive supercycles have a value of x that typically decreases from the initial step-size, x.sub.i, to a final step-size, x.sub.f. In a preferred embodiment the initial step-size, x.sub.i, is between one and one-eighth of the average dimensions of the asymmetric unit, and the step-size is decreased by integral angstrom units. In a more preferred, the initial step-size is between one-quarter and one-half of the average dimensions of the asymmetric unit. The final step-size, x.sub.f, is typically not much smaller than the resolution, K, of the supplied data.

This method of varying the step-sizes is designed to first perform large random movements in early supercycles that span the full size of the asymmetric unit, while subsequent smaller moves in later supercycles sample the asymmetric unit and experimental Fourier data more finely. Not unexpectedly, step-sizes smaller than about a third of the data resolution contribute little to the final outcome.

This procedure--from the first supercycle to the last one--is called the "condensing protocol." When the last supercycle is terminated, the calculated data resulting from the final distribution of scattering bodies will be highly correlated with the experimental data. For example, the value of the Pearson Coefficient, r, is typically in the range of 0.6-0.8, as compared with the typical values of r=0 for random starting configurations, and r=1 for an exact fit.

Direct visualization of the final distribution of scattering bodies is optionally performed on a display device using suitable molecular graphics computer software. Suitable display devices include cathode ray tubes or printed output. Methods for displaying the final distribution are well known in the art and include the use of computer graphics software such as FRODO, HYDRA, McIMDAD, MIDAS, MOGLI. Visualization of the final distribution of scattering bodies outlines the shape of the macromolecular envelope and shows low-resolution features. For proteins, structural motifs such as inter-domain clefts and other prominent surface indentations, are typically observed. The computing time for the condensing protocol is modest--many macromolecular problems requires roughly on the order of an hour of computer time on a mainframe computer such as a VAX 8550, produced by Digital Equipment Corporation of Maynard, Mass. Other computers are suitable for practicing the invention, the choice of which will be apparent to one skilled in the art, as described above.

Many of the parameters used with the present invention may be varied. For example, in other embodiments of the invention, the numbers of individual microcycles in a macrocycle, or the number of macrocycles in a supercycle are modified. The constraint that must be observed is that the number of microcycles must be sufficient to sample a sufficient number of the allowable moves during the condensing procedure. In other embodiments, the radii of the scattering bodies are modified concurrently with modification of the scheduling of step-sizes and subject to the constraints imposed by the experimental data.

As the condensing protocol proceeds and the scattering bodies condense into stable positions, the refined distribution of homogeneous and featureless scatterers may occupy the macromolecular volume or the solvent void. Using the invention method described herein, scattering bodies preferentially condense into regions of the asymmetric unit occupied by solvent and away from the regions occupied by protein thus modelling the negative image of the electron density. However, the scattering bodies occasionally condense into the regions of the unit cell occupied by protein. Not wishing to be bound by speculative theory, I believe that the condensation preference is determined by the inherent properties of the solvent, relative to that of the macromolecule.

Except at extremely low resolutions, the solvent void and macromolecular volume have an important difference: the macromolecule is not featureless, in contrast to the solvent void. That is, at most resolutions a macromolecule and, therefore, the associated positive image have internal variations in electron density that are due to the local electron density variation within atoms or event due to small regions of solvent. In contrast, the negative image, which corresponds to the regions of the asymmetric unit occupied by solvent, is essentially featureless since the solvent randomly and uniformly fills these regions. Atoms that lie within the molecular envelope of the protein diffract to give the positive image internal features, while atoms within the solvent space, due to their generally random orientations and locations, do not give a diffraction pattern. Thus, the distribution of scatterers used in this method, which are featureless and homogeneous, models the negative image more readily.

On occasions, however, the scatterers condense to model the positive image. Such solutions occur infrequently, and occur especially when a poor set of parameters is chosen. Thus, an incorrect number of scatterers, a very small initial step size, etc., may contribute to modelling of the positive image. This possibility of mistakenly interpreting a positive image for the negative one can be avoided.

Whether a given final distribution of scattering bodies (i.e., image) is positive or negative is assessed by different methods, such as a suitable density-based refinement. This type of refinement compares the final image with its inverse as a function of changing parameters. For example, the final distribution of condensed scatterers are assigned variable X-ray scattering parameters, such as the scattering cross-section, which are used to calculate a new set of amplitudes for both the image and its inverse. By first calculating the correlation of each set of calculated amplitudes with the experimental data, and comparing the correlations, there is a basis for discrimination, particularly at higher resolutions. The negative image generally has a higher correlation.

Another method for assessing the final distribution incorporates higher resolution experimental data. After the final distribution of scatterers is determined, additional higher resolution data is incorporated into the calculation of the correlation coefficient. Thus, the correlation coefficients between the final distribution and the experimental Fourier coefficients incorporating higher resolution data, and the correlation coefficient between the inverse of final distribution and the same data are calculated. As more of the high resolution data are incorporated into these calculations, one expects that the correlation coefficient of the positive image would indicate a worse fit at a faster rate than the negative image, due to the aforementioned increase in internal structure of the positive image obtained from incorporation of higher resolution data. Modelling of the solvent space is not as sensitive to higher resolution data.

In preferred embodiments, the results from the low resolution technique aids in high resolution structure determination in the following cases: (1) Traditional molecular replacement with full or partial models; (2) Phase extension based on three- or more fold non-crystallographic symmetry; (3) Available high-resolution SIR data; (4) two dimensional electron microscopy, where Fourier data are typically of higher resolution than the available direct phases and (5) verification of suitable heavy atom candidate data sets.

The following examples illustrate, but in no way limit the invention.

The method of the invention was applied to model the electron density of the protein, Elastase from Pseudomonas aeruginosa. The experimental data collected from previous diffraction studies indicated that 70% of the reflections were collected to a resolution of 2.0 .ANG. and to an accuracy of R.sub.sym of 3.5%. Examination of the diffraction data indicated that the protein crystallized in the P.sub.212121 space group having unit cell dimensions of a=124.4 .ANG., b=51.5 .ANG. and c=44.5 .ANG.. The packing fraction of the protein was estimated to be P=0.6, and because of the symmetry of the P.sub.212121 space group, each unit cell had four symmetry-related solvent voids that together accounted for an estimated 40% of the unit cell volume.

Using the experimental data set for the Elastase protein, the condensing protocol was utilized to model the electron density of the unit cell and determine the shape of the molecular envelope of the protein. The following parameters were used:

Number of alpha carbons, N.sub.ca =298;

Number of hard sphere scatterers, N.sub.hs =199;

N.sub.ref (number of reflections)=470;

Resolution, K=7 .ANG.;

Initial move size, x.sub.i =12 .ANG.; and

Final move size, x.sub.f =8 .ANG..

In an initial step, the 199 hard sphere scattering bodies, each having a radius of 1.5 .ANG., were placed in an initial random distribution by placing a first body in a random location in the unit cell, followed by placing a second body in a random location, and so on, until all 199 bodies had been assigned positions. The only constraint on this process, other than the requirement that each sphere lie fully within the boundaries of the unit cell, was the physical limitation that the outer surfaces of any two spheres not intersect. FIG. 3 shows the initial distribution of scattering bodies. Once the initial distribution was determined, the Fourier amplitudes and phases of hypothetically scattered X-rays were calculated based on the assumption that each sphere, although having a 1.5 .ANG. radius for the purposes of distribution and movement within the unit cell, scattered as it were a point scatterer. As is known in the art, this calculation is obtained by a trigonometric summation over the positions of the scattering bodies. The Pearson correlation coefficient was calculated from this initial distribution and had an initial value of r= -0.06, indicating a poor fit.

The condensing protocol, as described above, was utilized to refine the positions of the scattering bodies. After about 40 minutes of computation time on a VAX 8450 computer, the condensing protocol produced a final distribution of scattering bodies having a Pearson correlation of r=0.85. The final distribution of the unit cell, which represents the effectively maximized correlation, is illustrated in FIG. 4.

The structure of Elastase has previously been solved to a resolution of 2 .ANG.. Comparison of the results obtained from the condensing protocol with the previously solved structure verified the validity of the model. The alpha-carbon backbone of the solved structure was superimposed on the final distribution of scattering bodies, and is presented in FIG. 5, which clearly shows that the scatterers had preferentially condensed into the solvent void and defined the molecular envelope of the protein, even on the order of 10 .ANG..

In order to quantitatively assess the progress of the fit during the condensation procedure, a 10 .ANG.-resolution molecular envelope of the known structure of Elastase was made on a 2 .ANG. grid. The envelope was chosen such that 50% of the grid-points corresponding to the highest electron density were located within the envelope. The accuracy of the model resulting from the condensing protocol was assessed in terms of a spatial distribution ratio, j, which is defined as the number of scattering bodies exterior to this 10 .ANG. resolution envelope, divided by the number interior. For a random distribution of scattering bodies, j is approximately unity, and as condensing protocol proceeds and scatterers move from the interior of the molecular envelope to the exterior, j is increases. For analysis purposes, the value of j, as well as the Pearson correlation coefficient, r, were tabulated and plotted as shown in FIGS. 6a-b, for each macrocycle of the condensation protocol.

FIG. 6b shows that j increased to 1.93, while the correlation coefficient concurrently increased to 0.85 during the protocol. While this simple measure, j, approximately discriminates between solvent and macromolecule, it does not measure the spatial distribution of scattering bodies within either the solvent void or, more importantly, within the macromolecular envelope. In particular, FIG. 5 shows that the scattering bodies lying within the macromolecular envelope are preferentially distributed toward the macromolecular surface rather than clustered about the molecular centroid.

Finally, to determine whether the final distribution was a relatively unstable local maximum of the correlation function, a different random starting configuration and a different resolution range was used. The following parameters were used:

Number of alpha carbons, N.sub.ca =298;

Number of hard sphere scatterers, N.sub.hs =199;

N.sub.ref (number of reflections)=2195;

Resolution, K=4 .ANG.;

Initial move size, x.sub.i =11 .ANG.; and

Final move size, x.sub.f =2 .ANG..

The results are shown in FIG. 7.

Since more reflection data were used (as required by the resolution range used in this case), the computer processing time increased to 3.5 hours. The molecular envelope obtained with these parameters was sufficiently detailed to place the alpha carbon model of the protein within a few angstroms of its true center of mass. A limited translational and rotational R-factor search, as is known in the art of molecular replacement, identified the exact location of the molecule within the unit cell.

The experimental diffraction data, collected to 1.9 .ANG. with an R.sub.sym of 5.4%, indicated that the protein crystallized in the P.sub.212121 space group with unit cell dimensions of a=32 .ANG., b=37.5 .ANG., and c=44.6 .ANG.. The packing fraction of the protein was assumed to be P=0.6.

The crystal structure had been previously solved with a refined r-factor of 19%.

The parameters used in this example were:

Number of alpha carbons, N.sub.ca =63;

Number of hard sphere scatterers, N.sub.hs =43;

N.sub.ref (number of reflections)=291;

Resolution, K=5 .ANG.; and

Initial move size, x.sub.i =8 .ANG.; and

Final move size, x.sub.f =3 .ANG..

Because of the small size of the unit cell, higher resolution data was required to have sufficient over-determinacy. Furthermore, since the scattering bodies were not as numerous as in the previous examples and since the unit cell is small, all allowable moves could be sampled using fewer macrocycles. The minimum number of macrocycles per supercycle, therefore, was decreased from 40 to 20, which required less computer processing time. The correlation coefficient for the initial distribution was r=-0.1, and the initial value of j was 0.87. The index of spatial distribution, j, was calculated as described above but using a 6 .ANG. envelope of the Rlt69 model, made on a 2 .ANG. grid.

After about 10 minutes of computer processing time, the condensing protocol produced a final distribution of scattering bodies having an effectively maximized correlation of r=0.65, and a j=2.6. FIGS. 8-10 show the initial distribution of scattering bodies, the final distribution of scattering bodies, and the alpha-carbon backbone of the protein superimposed on the final distribution, respectively. FIG. 11 and FIG. 12 show the behavior of r and j during the condensation protocol.

Other embodiments are embraced within the present invention. For instance, substantially any crystalline molecule having uniformly diffracting voids within the unit cell can be modeled by the methods described herein. Examples include the modelling of the electron density of polymeric nucleic acids, such as DNA or RNA fragments, as well as nucleic acid-peptide complexes, virus particles and the like.

Other methods for maximizing the correlation between the experimental data and the calculated amplitudes include the downhill simplex method, direction-set methods (such as Powell's method), conjugate gradient methods (such as the Fletcher-Reeves and Polak-Ribiere algorithms), variable metric methods (such as the Davidon-Fletcher-Powell algorithm), simulated annealing, random-walks, and the like. Other methods that maximizes a correlation between the calculated amplitudes and the experimental data can be used.

Moreover, beyond the phaseless Fourier inversion problem of the present invention, other problems can be reformulated in a manner analogous to that presented above. Such as the determination of three dimensional structures using electron microcopy. In another preferred embodiment, the selection of the point scatterers to be moved is not random. For example, the selected point scatterer is the one having the largest increase in the correlation between the changed distribution and the experimental data. In another preferred embodiment, the selected point scatterer is not moved in a random direction, but in a direction corresponding to a predetermined algorithm. For example, the direction is selected to effectively maximize the increase of the correlation between the distribution and the experimental data.

Moreover, in other embodiments, a plurality of point scatterers are simultaneously moved in one cycle. The number of point scatterers that are moved simultaneously is between one and all of the point scatterers present. In other embodiments, the distance that a point scatterer is moved in a cycle is variable. That is, the distance moved is not a fixed, pre-determined distance. In a preferred embodiment, the move distance is randomly chosen, but constrained in a predetermined range. Alternatively, the correlation coefficient versus the move distance is calculated, and the scatterer is moved to maximize the correlation. Rather it is varied dynamically in a manner that is consistant with the other aforementioned standard maximization methods. However, the move distances cannot be far different from the appropriate pre-determined distances that are characterized.

The present invention provides new methods for modelling the electron density of a macromolecule in a crystal lattice. It is to be understood that the above description is intended to be illustrative and not restrictive. Many variations of the invention will become apparent to those of skill in the art upon review of this disclosure. The scope of the invention should, therefore, be determined not with reference to the above description, but instead should be determined with reference to the appended claims along with their full scope of equivalents.