JP2003524831A - System and method for exploring combinatorial space - Google Patents

Abstract

(57) [Summary] A method and system for exploring every corner of a combination space without causing a combination explosion. The search is performed according to at least one desired characteristic of the combination, which can be converted, for various combinations of the basic elements, into a quantitative measure of the success of the search. Since the number of variables, and thus the number of combinations, can be very large, preferably a sample of the combinations is examined. Elements of the combination that have consistent maximization and / or enhancement of quantitative measurements are retained, and other elements are discarded. This process may be repeated until a certain minimum number of combinations is found, which may optionally be further evaluated according to similar parameters and / or some other parameters or characteristics.

Description

Detailed Description of the Invention

FIELD OF THE INVENTION [0001] The present invention discloses systems and methods for searching combinatorial spaces in depth, and in particular, the position of one or more combinations of primitives having desired properties. And systems and methods that can be identified within. The desired property must have a numerical basis or be at least convertible to some form of numerical measure and / or equivalent. The present invention, according to the above properties, allows a quick and efficient search of the combinatorial space without causing combinatorial explosion. The present invention accomplishes the above task by examining each value of a primitive element at least once, and preferably multiple times during the search process. Therefore, it can be said that each value can be searched for in an exhaustive search without having to exhaustively search all combinations of values for the basic elements. Thus, the present invention combines the effects of a non-exhaustive stochastic search process and the effects of an exhaustive search.

BACKGROUND OF THE INVENTION This section has been divided into a number of subsections for the sake of brevity. Put simply,
The first item describes the general problem of combinatorial spaces and the search within that space. The next section describes the solutions that have been attempted so far to solve many different biological problems, which are background techniques for dealing with combinatorial search spaces in relation to biological problems. It is also an example showing the inadequacy of the measures. These items include the placement of polar protons for biological molecules such as proteins, the placement of side chains for amino acids in proteins, and the prediction of loop structures in proteins.

Combinatorial Space A combinatorial space is defined as having multiple combinations of basic elements. These combinations may differ depending on the format of the element values, the resulting structure of the element combinations, or may be generated as a result of both factors. At a more basic level, each combination can be considered to consist of variables, each of which can have one or more values. In the present specification, each variable preferably takes on one value in a set of disjoint values, but instead each variable may, for example, be in a range of continuous values or in a function. May take one value of

Many basic biomaterials are made in biology because they are themselves made up of relatively basic building block combinations, even though the resulting structures and / or functions are very complex. Often results in combinatorial spaces. Examples of combinatorial space include, but are not limited to proteins. Proteins are produced by the combination of amino acids as building blocks and ultimately fold into a conformation known as the "tertiary structure", which is a set of values in the combinatorial space. Searching every corner of such a combinatorial space to find this single structure can also be called a “combinatorial search”. However, folded proteins in the biological milieu may exist in many equilibrium conformational metastable states without being locked into one "tertiary structure". Therefore, a combinatorial space must be searched everywhere, and preferably more than one solution must be found.

Even though these various types of building blocks are obviously simple, the exhaustive search of the combination space is not exhaustive because of the enormous number and complexity of combinations available. Since it exceeds the processing range of, it is a difficult problem to search every corner of the combination space. For example, even simple small proteins with relatively short amino acid sequences have a huge number of different potential structures, even if only one or a few are actually stable structures. In terms of protein structure, this problem includes the placement of amino acid side chains and polar protons (to determine hydrogen bonds within the protein and optionally between the protein and other molecules), and loops. It consists of several subproblems, such as the location and type of the larger structures in the protein. Therefore, it has been found that exploring combinatorial spaces of the above type in depth, especially due to biological problems, typically also withstands modeling and prediction by various computational approaches.

Many attempts have been made to address the problem of searching combinatorial spaces in every corner. Generally, these attempts include transforming a large search space into a smaller search space by subdividing the combinatorial space, changing the display of the space, and reducing the determinism of the criteria for finding potential solutions. And so on, including a directed search in which the search strategy was guided towards the desired solution. However, it has been found that none of the above solutions are suitable for biological problems such as prediction of protein structure. Within the protein structure there are many "dead ends (
Directed search is not useful for problems such as the one above because of the possibility of "dead-end)" and the lack of a set of decisive rules for protein folding. Similarly, the search space cannot be reduced at this time because the laws for protein structure are unknown. Also, because the protein structure has basic fixed building blocks of amino acids, the display of the combination space cannot be changed. Finally, identifying less deterministic solutions proves to be less useful for the overall prediction of protein structure. This is because the definite "rule" itself that can be guided to fold a protein is specific to that protein at the moment and has not been generalized. There is no effective search strategy that can find a population of appropriate value sets in the combinatorial space.

First Biological Problem: Polar Protons The generally attempted solutions above can be further described in terms of specific solutions to particular biological problems. For example, the positions of polar protons (hydrogens) are important for determining hydrogen bond relationships and specificities within biologically important molecules such as proteins and DNA molecules, as well as for these and other molecules. It is also important for determining hydrogen bonds such as between and. Inclusion of all hydrogen atoms in protein and nucleic acid models is necessary for energy minimization, for a more accurate representation of biological systems during molecular dynamics simulations, and for understanding molecular recognition (Jones, et al., J. Mol. Bio
1, 24, 43-53, 1995). Polar hydrogen plays an important role in determining secondary structure and protein packing, and the correct placement and reliable formation of hydrogen bonds is very important for energy evaluation. One misplaced polar hydrogen in the active site has been found to dramatically change the conformation of the substrate during molecular dynamics simulations (bass (B
ass) et al., Proteins, 12: 266-277, 1992).

At present, the main source for obtaining high resolution data of biomolecules is X-ray crystallography. However, X-ray crystallography is useful for locating heavy atoms when protons have not been localized. Although neutron diffraction tests can locate protons, at the moment only a few X-ray / neutron diffraction combined tests have been deposited with the Protein Data Bank (PDB).

Several computer-based methods have been proposed for locating polar hydrogens in proteins. First, common to most molecular modeling software packages, hydrogens can be placed in a non-specific manner, after which the structure can be optimized by an energy minimization algorithm subject to the multiple minimum problem. Energy minimization algorithms do not take into account alternative positions for flexible polar hydrogens, many of which may form hydrogen bonds.

[0010] Brunger and Karplus (Prota)
Ins, Vol. 4, pp. 148-156, 1998), a second method utilizes the search for a local minimum conformation at each polar proton turn by twist angle rotation. Then, this iterative process is continued until the convergence point is reached. Since this method does not take into account the effects of adjacent rotatable hydrogens, it may be said to be accurate only for systems where there is no such intimate hydrogen-hydrogen contact.

Bass et al. (Proteins, 12: 266-277, 19)
1992) proposed a third method based on partitioning the system into a network of interacting hydrogen-bond donors and acceptors. The algorithm tries to maximize the number of hydrogen bonds that can be formed in each network and minimize the total distance between donor and acceptor. Since each network can probe the most likely set of hydrogen bonds, the number of comparisons (between options) is commensurate with the factorial of the elements in the network, which limits the computation to small networks. End (Bass et al., Proteins, Volume 12, 2)
66-277, 1992). No energy rating has been used to select the best structure. As a result, the output may contain high energy interactions between the localized hydrogens and their environment.

Richardson et al. (J. Mol. Biol., 2nd
85, 1711-1733, 1999) and Word et al. (J. M.
ol. Biol. , 285, pp. 1735-1747, 1999), by recently considering the Asn / Gln "flip", the protonation state of the histidine ring, and a simple model of water interactions. It extends the approach. For rotatable protons, optimize the set of local H bonds with respect to distance and van der Waals overlap. Unfortunately, none of the above solutions to the problem of the placement of polar protons in biomolecules such as proteins
Even for a class of such molecules, it cannot be generalized to be accurate within that class and effective for practice.

Second biological problem: Arrangement of amino acid side chains As another example of the problem in searching the combinatorial space,
The arrangement of amino acid side chains may be mentioned. This problem itself is solved by combinatorial space, but this is only part of the general problem of protein structure prediction. However, this problem has been found to be difficult to handle with currently available methods that attempt to predict the position of the side chains.

Accurate placement of protein side chains is essential for both theoretical and experimental purposes. From a theoretical perspective, this is a subproblem in de novo protein structure prediction. Also, structure-based drug design (Defay and Cohen, Proteins, Volume 23, 4).
31-445, 1995), reverse folding and threading (th
reading algorithm (Bahar and Jerigan)
nigan), J. Mol. Biol. 266, 195-214, 19
1997) Understanding of folding process and structural stability (Zhukov (Zhukov)
) Et al., Protein Sci. , Vol. 9, pp. 273-279, 2000),
Ab-initio prediction of protein tertiary structure (Huang (H
Uang) et al., Proteins, 33, 204-217, 1998).
, And homology-based modeling (Brundell et al.,
Nature, 326, 347-352, 1987), which is an unavoidable problem. From the viewpoint of an X-ray crystallist, it is considered that the position determination of the side chain using the electron density map of the main chain, which is performed prior to the refinement calculation, can be speeded up. The main limitation is the large number of possible conformations that each side chain can have (Lee and Subbiah, J. Mol.
Biol. 217, 373-388, 1991). An exhaustive search for all possible protein conformations is beyond the state of the art computer.

X-ray crystallography usually provides a single structure characterized by the R factor. The crystal structure describes the biomolecules in a highly ordered crystal lattice, as opposed to the solution environment of the more physiologically relevant NMR structure. X-ray crystal structure may not be found in the ensemble of conformations in solution,
It may be biased to a special conformational sub-state in the crystal (Brunger, Nat. Struc. Biol., 4 supl.
, 862-865, 1997). Observing alternative rotors exceeds the detection limits of conventional X-ray crystallography techniques, unless at very high resolution. At least 10% of all side chains in a protein adopt multiple disjoint conformations in a carefully refined crystal structure (Smith et al., B.
iochemistry, 25, 5018-5027, 1986). MacArthur and Thornton (Ac
ta. Cryst. D Biol. Cryst. , D55, pp. 994-1004, 1999) found a significant and unexpected correlation between χ ₁ mean and resolution, mainly for small flexible side chains. All data support the hypothesis that this observation reflects local conformational flexibility and disorder that would be interpreted as a single distorted conformer at low resolution.
The results of all these studies have focused on a diagram of protein structure, which is dynamic rather than static, and the need to extract this dynamic information from the NMR population to gain a more detailed understanding of protein function. (Philippoul)
and Lim, Proteins, Vol. 36, pp. 87-110, 1.
999). Protein function and molecular recognition depend on structural plasticity (Garcia et al., Science. 279: 1166-1172).
P., 1998, conformational flexibility of receptor proteins is one of the major factors affecting ligand docking (Desme.
t) et al., FASEB J .; 11: 164-172, 1997). However, due to the large number of minimal energy conformers on the potential energy surface, computer-accurate determination of protein side-chain positions is a complex task, even with a rigid backbone. Conventional methods for side chain addition usually result in a single protein structure, which is compared to the X structure when available. The conformational space is ignored. Recently, many protein NMR studies have been carried out, suggesting many conformations for each protein (Schneider).
) Et al., J. Mol. Biol. 285, 727-740, 1997).
However, whether such conformations represent an alternative solution to the distance limitation emerging from two-dimensional and three-dimensional coupling maps, or whether those conformations can contribute to the total population at equilibrium It is also unclear whether it is a conformation of. The classical molecular dynamics (MD) approach is the technology of choice for simulating biomolecules. With current technology, it is becoming commonplace to perform MD simulations of systems consisting of tens of thousands of atoms in a few nanoseconds (Sagui and Darden,
Ann. Rev. Biophys. Biomol. Struct. , Volume 28,
155-179, 1999). However, the range of time scales associated with biomolecular function can range from nanoseconds to more than seconds. By MD, the time required to reach equilibrium between the different conformers of the protein is extraordinary for such simulations, and we can only glance at the behavior of the protein in its environment. Seem.

Current strategies for side chain addition fall into three different categories.
The first category is the conformational space of each side chain. Continuous space method (Eisenmenger, J. Mol. Biol., No. 231
Vol. 849-860, 1993; Reitberg and Elber, J. Am. Chem. Phys. Volume 95, 9277-9287
Page, 1991), the torsion angles of any side chain can be sampled. The discrete space method is based on the assumption that the side chains exist in an energetically favorable conformation called rotamers, which rotamers are
Local minimum conformers collected by statistical analysis of known structures (Chandracekaran and Ramachandran).
machandran), Int. J. Protein Res, Volume 2, 22
3-233, 1970; Sasisekharan and Ponnuswamy, Biopolymers, vol. 9, 12
49-1256, 1970; Sasisekharan.
And Ponnaswami, Biopolymers, Vol. 10, pp. 583-592, 1971; Ponder and Richards (R).
Richards), J. Mol. Bio1. , 193, 775-79l,
1987; Gelin and Karplus, Bioc.
Chemistry, Vol. 18, pp. 1256-1268, 1979; Dunblack and Karplus, Nat. Struc
t. Biol. Vol. 1, pp. 334-340, 1994). Discrete space methods cannot predict conformations that do not exist in rotamer databases. However, large rotamer databases containing very rare conformations do not always yield better predictive results than smaller databases (Holm and Sander, Proteins).
, 14: 213-223, 1992; Laughton,
J. Mol. Biol. 235, 1088-1097, 1994; Tanimura et al., Protein Sci. , Volume 3, 2358
~ 2365, 1994; Vasquez, Biopolym.
ers, 36, 53-70, 1995). Databases can also be classified into backbone-dependent and backbone-independent. The former is based on the relationship between the side chain conformation and the local backbone conformation, while the latter is not.

The second category is the cost function for evaluating the solution. Energy-based methods include non-bonded terms (Laughton, J. Mol. Biol,
Volume 235, pp. 1088-1097, 1994; Baskets (Vasquez
), Biopolymers, 36, 53-70, 1995; Wilson et al., J. Am. Mol. Biol. , Volume 229, 996-100
6 1993; Vasquez, Curr. Opin. St
ruct. Biol. , Vol. 6, pp. 217-221, 1996). We assume that the lower the energy, the higher the prediction accuracy.

The following knowledge-based methods have also been proposed. Sutcli
ffe) et al. (Protein Eng. Vol. 1, 385-392, 1987) from side chains at topologically equivalent positions (provided such a correlation is observed). It suggests how to assemble side chains using spatial information and the most likely side chain conformation in each secondary structure type. Sali and Brundel 1 (J. Mol. Bio
l. , 234, 779-815, 1993), describes comparative protein modeling methods designed to align with related structures and find the most likely structure for a sequence. There is. Bower et al. (J. Mol. Biol., Vol. 267, pp. 1268-1282, 1997).
Identified the location of residues within the most desirable backbone-dependent rotamers and systematically analyzed collisions caused by the structure.

The third category is search strategies. There are many examples of search strategies used. Metropolis Monte Carlo method (Holm and Sander, Proteins, 14: 213-223, 1
992), Gibbs Sampling Monte Carlo (Vasquez)
, Biopolymers, 36, 53-70, 1995), neural networks (Hwang and Liao, Protei.
n Eng. , Vol. 8, pp. 363-370, 1995), Genetic algorithm (Tuffery et al., J. Biomol. Struct. Dy.
nam. Volume 8, pp 1267-1289, 1991; Tuffee
ry) et al. Comput. Chem, 14: 790-798, 199.
3 years), simulated annealing (Lee) and Subia (Subbi)
ah), J. Mol. Biol. , 217, 273-288, 1991), mean field optimization (Koehl and Delarue, J.
． Mol. Biol. , 239, 249-275, 1994) and L.
ES (Locally Enhanced Sampling) method (Roitberg and Elber, J. Chem. Phys.
95, 9277-9287, 1991).

Combination search (Dunblack and Karplus)
lus), J. Mol. Biol. 230, pp. 543-574, 1993.
Year; Tuffery et al. Biomol. Struct. Dy
nam. 8: 1267-1289, 1991; Wilson
on) et al. Mol. Biol. 229, pp. 996-1006, 199
3 years) may be used for the individual conformers, followed by successive minimizations at the final stage of refinement. There is no guarantee that any of the above will converge to an effective solution. Another widely used method is the dead end exclusion method (
DEE: Dead End Elimination). This method is based on identifying rotamers that are not completely compatible with the global minimum energy conformation and excluding rotamers that cannot contribute to a local energy minimum above a certain order. A conformation consisting of such rotamers can be considered dead-end (Desmet (Desmet)).
et al., Nature, 356, 539-542, 1992; Desmet et al., FASEB J. et al. , Vol. 11, pp. 162-172,
1997; Lasters and Desmet, P.
protein Eng. , Vol. 6, pp. 717-722, 1993). A global minimum can be found if sufficient rotamers are excluded by repeated application (Goldstein, Biop.
hys J.H. 66, 1335-1340, 1994). However, DEE is unable to find a population of low energy solutions.

The A ^* algorithm finds an optimal path from a root node of a search tree to a goal node by using a cost function shown by f ^* (Leach and Lemon, Proteins, Vol. 33, Pp. 227-239, 19
1998). Each node has a unique f ^* value consisting of the cost of searching the node from the start node and the estimated cost of reaching the goal node. Optimize f ^* iteratively, expand the node for which f ^* has a minimum value, and calculate a new value of f ^* for the successor node.

[0022] for identifying the low energy side-chain conformations of the proteins, the best way known at present is a method that combines the DEE and A ^* algorithm, A ^* algorithm, build a distribution function Has been used to The A ^* algorithm approach can find the best N solution, but is limited to relatively small proteins. Currently, the largest protein analyzed by this algorithm contains 68 amino acids, which consist of about 10 ⁴³ combinations, depending on the complexity of the rotamer library, although more combinations are possible. The proteins that they have are also common. As a "stand-alone" algorithm (without the DEE pre-processing stage), the A ^* algorithm reaches up to 10 ²¹ combinations. For an efficient search by the A ^* algorithm, one must have a good cost estimate to reach the goal node. This is problematic due to the interaction between residues that have not yet been assigned. These limitations require new and powerful algorithms to find global minimums and lowest energy conformations in larger systems. Unfortunately, no such algorithm is currently provided.

Third Biological Problem: Prediction of Loop Structure As mentioned above , prediction of protein structure requires search in combinatorial space, and there is no suitable solution at this time. The prediction of protein structure itself is divided into many smaller problems, and even those small problems require a search in the combinatorial space. One example of such a problem is the prediction of very complex loop structures.

Structural genomics projects have been undertaken to provide experimental structures or good models for newly discovered sequences from various genomic projects. Brenner and Levitt (Pro
tein Sci, vol. 9, 197-200, 2000), the number of new folds has been gradually reduced based on the analysis of sequence similarity databases.
It suggests that the most common folds are readily found. Therefore,
Another protein in which one protein has a known tertiary structure (“template”)
Homology modeling may be a powerful tool for predicting the three-dimensional structure of a protein sequence if it shows a reasonably high sequence similarity with. In that case, the elements of secondary structure are passed from the template to the target protein. However, the "loop" or "coil" elongation is not determined and must be predicted. Homology modeling has long been used and has been successful. Rous sarcoma virus aspartic protease (Weber, Sc
H.I., 243, 928-931, 1989).
Prediction of one protease and amyloid precursor protease inhibitor domain from bovine pancreatic trypsin inhibitor (Struthers)
Et al., Proteins, Vol. 9, pages 1-11, 1991), and lessons from breakthrough predictions have been useful for subsequent sequence homology modeling for many structures. Most of the predicted loops are variable in length and sequence within the protein family, and their modeling requires a three-dimensional insertion and deletion process. Even when such sequence and length discrepancies are localized to short stretches of protein, significant rearrangement of the backbone conformation can occur. The other regions of the two proteins that are highly similar in sequence and length are called "framework regions." An aperiodic structure that connects two consecutive secondary structures is called a "loop" and is described as having an "irregular conformation" or a "random coil" (Oliver (Ol
iva) et al. Mol. Biol. 266, 814-830, 199.
7 years). Finding appropriate coordinates for variable regions in protein homology construction by loop construction, ie insertion and deletion, is an important task in globular proteins (Alwyn and Thirup).
up), EMBO J. 5: 819-822, 1986; Brooks et al., J. Am. Comput. Chem. , Volume 4, 187-217
Page, 1983; Bruccoleri et al., Nature, 335, pp. 564-568, 1988; Bruccoleri.
i) and Karplus, Biopo1ymers, Vol. 26, 1
37-168, 1987; Geer, Proteins, Volume 7,
317-334, 1990; Palmer and Shelaga (Sch.
eraga), J. et al. Comp. Chem. , Vol. 12, pp. 505-526, 19
1991; Summers and Karplus, J. et al. M
ol. Biol. 216, pp. 991-1016, 1990). The construction of loops is also a major issue in other fields of structural biology. This is loop 2
It is a tremendously complex combinatorial problem that involves finding a properly insertable fragment between two endpoints and evaluating its energy.

Extensive structural studies have been carried out on immunoglobulins having almost the same sequence and structure but having markedly different binding specificities (Bruccoli (Brucco).
Leri) et al., Nature, 335, 564-568, 1988; Chothia et al., Science, 233, 755-758.
1986; Chothia and Lesk, J. et al. Mol. B
iol. 196, 901-917, 1987; Fine.
Et al., Proteins, Vol. 1, 342-362, 1986; Tramentano and Lesk, Proteins, Vol. 13, 231-245, 1992). The specificity is due to six hypervariable loops called complementarity determining regions (CDRs). Understanding the structure of these loops teaches us much about antigen binding, catalytic antibodies, and molecular recognition.

As another example, transmembrane helices within a G protein-coupled receptor (GPCR) are linked to form a latent ligand binding domain on the extracellular side, or a G protein-coupled on the intracellular side of the membrane. Examples include intracellular loops that form domains (Kazmi et al., Biochemistry, Vol. 39, 3734-).
3744, 2000; Iiemann et al., J. Am. Struc
t. Bio1. 128, pp. 243-249, 1999). Prediction of the three-dimensional conformation of the intracellular loop is important for understanding the mechanism (Gather et al., Nature, 362, 345-34).
8 p. 1993; Kyle et al. Med. Chem. 37: 1347-1352, 1994; Cypes et al., J. Am. Bi
ol. Chem. 274, 19455-19464, 1999; Zhang et al., Protein Sci. , Vol. 3, pp. 493-506,
1999; Lee et al., J. Am. Biol. Chem. , Vol. 275, 92
84-9289, 2000), and is subsequently important in the design of GPCR-related drugs (Mukherjee et al., J. Biol. Chem.
． 274, 12984-12989, 1999).

Similar problems arise in modeling chemical rings, but the constraint in this case arises from the need to close the rings with chemically reasonable bond lengths and angles (go (
Go) and Shuraga, Macromolecules, 3rd
Vol. 178-187, 1970; Bruccoleri and Karplus, Macromolecules, Vol. 18, 2
767, 1985; Shenkin et al., Biopolyme.
rs, 26, 2053-2085, 1987; Palmer.
) And Sheraga, J .; Comp. Chem. , Volume 12, 5
05-526, 1991).

Many current strategies divide the above problem into two sub-problems. First, find a geometrically acceptable conformation for inserting a polypeptide backbone fragment of appropriate length between two endpoints of a loop within a framework of known protein structure. Must Go (Go and Scheraga, Macromolecules, Volume 3, 178-18
P. 7, 1970; Weiner et al., J. Am. Amer. Chem. S
oc. , 106, 765-784, 1984; Bruccoli.
coleri) and Karplus, Macromolecule
S., Vol. 18, pp. 2767, 1985; Shenkin et al., B.
iopolymers, 26, 2053-2085, 1987). This step usually results in multiple solutions. Second, the appropriate polypeptide is selected from the solutions suggested in the first step, usually by energy criteria.

A method by lattice search for most of the backbone dihedral of the loop (Bruccoleri and Karplus, Biopoly)
mers, 26, 137-168, 1987; Mault.
James, Proteins, Volume 1, 146-163,
1986), in which the first grid point is selected from the allowed region of the Ramachandran map. The allowed regions are those determined from the distribution of backbone twist angles observed in the set of protein structures determined by high resolution crystallography. Such a method is limited to relatively short loops because the number of degrees of freedom of the system increases exponentially in the lattice search.

The database search was first done by Jones and Thirlap.
Rup) (EMBO J., Vol. 5, pages 819-822, 1986), proposed by Summers and Karplus (J.
Mol. Biol. 216, pp. 991-1016, 1990). A high resolution X-ray structure database is consulted to select amino acid segments with similar geometric descriptors and sizes as the desired loop. Dock selected segments into proteins and examine their feasibility using geometric and energy criteria. Koehl and Delarue
) (Nat.Struct.Biol., Vol. 2, pp. 163-170, 1995).
Years) to add side chains, self-consistent
nt) Database search combined with mean field theory was used. It has been shown that the loop length poses a limitation for reliable database searching. That is, Summers and Karplus (J. Mol.
Biol. , 216, 991-1016), this method is limited to 6 residues in length. Furthermore, this method needs to cover most of the test cases, which requires a fairly large database of well-parsed structures. Deane and Brundel
(Proteins, 40, 135-144, 2000) recently adopted an exhaustive appetition search on computer generated fragments to generate loops of up to 8 residues.

Tree Search Method (Brooks et al., J. Comput. Chem.
． , Vol. 4, pp. 187-217, 1983; Bruccolei.
ri) et al., Nature, 335, pp. 564-568, 1988), the search strategy is node-based, which may be expanded during the conformational search. Yields are low for medium-sized structures and prohibitively low for large structures (Sherkin et al., B.
iopolymers, 26, 2053-2085, 1987).

Random tweak method (Sherkin
Kin) et al., Biopolymers, 26, 2053-2085, 19
1987), a structure that is not constrained is generated, in which all dihedral angles are set to random values. Subsequently, the loop constraint is geometrically strengthened by collectively tweaking all dihedral angles in the iterative process. Fine et al. (Proteins, Volume 1, 342)
~ 362, 1986) also generate a number of random conformations to the backbone of the desired loop, followed by minimization and / or molecular dynamics with the remaining molecule immobilized. Kyle et al. (J. Med.
Chem. , 37, 1347-l352, 1994), combining homology modeling of known transmembrane structures of bacteriorhodopsin with a two-step conformational search for energy minimization, molecular dynamics and docking simulations. I am using. These techniques search the conformational space near the starting point for the local energy minimum.

Bond scaling-relaxation
procedure meets both geometric and energy requirements at the same time (
Zheng et al., J. Am. Comp. Chem. , Volume 14, 556-56
P. 5, 1993). Generate a random initial conformation with standard bond lengths and angles. The bond length for each initial conformation is increased or decreased to fit the distance constrained by the loop, relaxing the system to a local energy minimum. This method will be described later in terms of the multiple copy sampling method.
Enhanced by combination with a copy sampling method (Zheng et al., Protein Sci., Vol. 2, 124.
2-1248, 1993; Zheng et al., Protein Sc.
i. , Vol. 3, pp. 493-506, 1994). The improved method was used to process loops of up to 12 residues.

A significant problem with all of the above methods is the large number of different closures covering a large conformational space that are needed to maximize the likelihood that the correct (very close) structure may be included. It is difficult to give (Zh
eng) et al., J. Comp. Chem. 14: 556-565, 199
3 years). Current approaches are either limited to small loops or only a section of conformational space. Therefore, potentially good solutions may have been overlooked. There is a need for a more effective search strategy that considers the entire conformational space and finds all possible conformations that follow the geometrical criterion of loop closure. These solutions can later be evaluated on a finer scale.

Other Biological Issues In addition, there are many other biological issues that may require combinatorial search, such as those associated with theoretical drug design. The basic assumption for theoretical drug design is that the drug action is molecularly bound to the pocket of another molecule where one molecule (ligand) is usually a larger molecule (receptor, generally a protein). Is obtained. In their active or binding conformations, molecules exhibit geometric and chemical complementation, both of which are essential for successful drug action.
By binding to such macromolecules, drugs can modulate signaling pathways, for example, by altering their sensitivity to hormonal action or by altering metabolism, for example by interfering with the catalytic activity of enzymes. Most commonly, this is done by interfering with the access of the natural substrate by binding within the specific cavity (active site) of the enzyme that catalyzes the reaction. In other cases, such as a transmembrane protein, an "antagonist" may be designed to block the binding of an "agonist" (a natural molecule that activates signal transduction), or
More potent binding agonists may be required as drugs in reducing the biological response.

Modeling molecular structure is a complex task, as most molecules are flexible and can assume many different conformations of similar or close energy states. Modeling the binding process is also a complex task, because the properties of the receptors, ligands, and the solvents in which they are found must be taken into account. Chemists are struggling to get as accurate a model as possible, but in practice some approximations have to be made. Clearly, the more accurate models used, the more opportunities chemists face in predicting molecular interactions. Nevertheless, a large number of predictions made using the fitted model have been confirmed by experimental observations. Recently, some drugs have been designed by computer theory. Encouraged by this, researchers have built tools that use approximate models and have investigated how useful these tools are. Such approximate models pose difficult algorithmic problems. More accurate molecular modeling, obtained through deeper theoretical understanding or increased computational power, can only improve the technology developed by simpler models.

Depending on whether the chemical and geometric structure of the receptor is known, the problems that arise can be divided into two broad categories. If the receptor is known,
Chemists are interested in finding out if it is possible to place the ligand within the binding pocket of the receptor in a conformational low energy conformation. This problem is called the docking problem. There are several variations to this, which may require an accurate explanation of binding interactions, or searching for an approximate estimate of which ligand fits within the receptor from a large database of ligands. May be issued.

Very often the binding pocket is unknown. In fact, X-ray crystallography or N
The number of three-dimensional structures determined by MR technology is relatively small, although the number is rapidly increasing. In this case, an indirect approach with many ligands that interact with specific receptors must be taken. These ligands have been discovered primarily by experimentation. Using the geometrical and chemical properties of these molecules, chemists try to infer information about receptors. Of particular interest to chemists is the identification of pharmacophores that are present within these ligands. The pharmacophore is the set of features in a particular three-dimensional arrangement that is contained within all the active conformations of the molecule considered. The predominant hypothesis is that the pharmacophore is the part or part of the molecule responsible for drug activity, the rest of the molecule being the underlying structure for the features of the pharmacophore. Once the pharmacophore is determined, by examining the different activities, relative shapes and chemical structures of the starting molecules, chemists can use it to design more potent drugs.

Techniques that have been used so far for computer-based drug design include robotics (kinematics and planning), graphic algorithms (visualization of molecules), geometric calculations (surface calculation), numerical methods (energy). Minimization), graph theory methods (identification of invariants), randomization algorithms (conformation search), computer vision methods (docking), and various other techniques such as genetic algorithms and simulated annealing. Numerous tools are now available to perform complex geometry and energy calculations, and the success of these computer-aided methods is being evaluated.

The other general problem in drug design is that which emerges from biomolecular targets. Advances in genomics, proteomics, and bioinformatics are rapidly creating new therapeutic targets for drug discovery research. Due to the myriad of compounds that may be tested for activity against these targets, biopharmaceutical research is rapidly becoming dependent on synergistic approaches to accelerate drug discovery. .. Novel computational methods are required to improve our ability to process much of the information emerging from the sequences of new proteins in order to convert sequences into structures. Moreover, most X-ray studies yield single structures lacking important information about protein and solvent proton positions, but even detailed protein structural studies have at best the limited set of static consonants. Since it only creates formations, it is limited to information content. Even if the complete structure of the targets is known, designing drug candidates that may interact with these targets is a daunting task. Thus, in the field of structure-based drug design, there is a great need for improved methods that facilitate the above tasks.

The most severe limitation of most of these problems is the high complexity due to the large number of variables. The search for any "real" molecular structure requires "optimization" of this large number of variables. Such a complex potential energy surface (PES) has many minima, one of which is the "global minimum" and is probably associated with the native structure of the molecule.

Despite recent significant advances, the general problem of global optimization remains unsolved (Wales and Scheraga, Sc).
ience, 285, 1368, 1999). Traditional minimization techniques are time consuming and tend to converge to a local minimum. Sampling of the "phase space" of biomolecules (Bern and Straub, Curr)
． Opin. Struct. Biol. , Vol. 7, p. 181, 1997),
It may help to search for the smallest areas and reduce the time required to reach such areas. Some of the main global minimization algorithms
The order of applicability to biological applications using computers is shown.

Since the number of conformations increases exponentially with the number of residues, protein structure prediction can be shown to be a NP-hard problem. The native conformation of proteins accounts for only a very small sub-population of these, thus requiring powerful search algorithms for exploration.

Simulated Annealing (SA) is a generalization of the Monte Carlo method for investigating equations of state and frozen states of many systems (Metropolis (Metr).
opolis), J. Chem. Phys. 21: 1087, 1953
Year). In the annealing process, the initially high temperature unaligned melt is slowly cooled. As the cooling progresses, the system becomes more aligned and approaches the "frozen" ground state at T = 0. The initial structure is perturbed and the energy change dE is calculated. If the change in energy is negative, new structures are allowed. If the energy change is positive, this is the Boltzmann coefficient exp- (dE / kT
) Is allowed based on This process is repeated a sufficient number of times to give good sampling statistics for the current temperature, then the temperature is reduced and the whole process is repeated until a frozen state at T = 0 is achieved. SA is particularly suited to large-scale optimization problems where the required global minimum is hidden in many smaller trivial local minimums (Holm and Thunder (
Sander), Proteins, Vol. 14, pp. 213, 1992; Lee (
Lee) and Subbiah, J .; Mol. Biol. , Volume 217,
373, 1991; Hwang and Liao, Prot.
ein Eng. , Volume 8, p. 363, 1995; Presse (Presse
t) et al., "Numerical Recipes", p. 326, Cambridge University Press, New York, NY, 1986).

Genetic algorithms (GA) have been applied to many optimization problems, some of which have been successful (Tuffery et al., J. Compu.
t. Chem. , Vol. 14, p. 790, 1993). GA is inspired by Darwin's theory of evolution, in which the most suitable one in nature is selected and survives (Forrest S (1993); Sci.
ence, Vol. 261, p. 872). Each time the GA is repeated, a competitive choice is made to eliminate weak solutions. High "fitness" by exchanging part of the solution for another
Solutions that have are recombined with other solutions. The solution is also "mutated" by making small changes to a single element of the solution. GA is simple, does not easily get stuck in the local minimum, and can often find the optimal solution globally. There is no need to perform any derivative or other problem specific calculations. However, there is no guarantee that it will converge to a valid solution, and it will require a large number of iterations to reach the convergence condition.

Tabu Search (TBS) (Glover, Computers
and Operations Research, Vol. 5, pp. 533, 19
1986) outperforms SA in both the time required to obtain the solution and the quality of the solution (Cvijovic and Klinowski).
ki), Science, 267, 664, 1995). At the time of initialization, the purpose is to roughly examine the space of the solution, and as candidate positions are specified,
The focus of the search is on obtaining the local optimal solution. TBS is problem-independent,
It can be applied to a wide range of tasks. TBS is easy to implement and occupies only a few lines of code throughout the procedure. And, conceptually, it is much simpler than SA or GA. However, TBS cannot guarantee that it will solve the problem of multiple minimums in a limited number of steps, which may increase the computation time.

H. The group at Scheraga is actively researching to devise methods for global optimization. Latent functional transformation and smoothing methods (Piela et al., J. Phys. Chem., 93, 3339, 1989: Pillardy and Piela, J. et al.
Phys. Chem. Vol. 99, p. 11805, 1993; Philady (P
illardy) et al. Phys. Chem. , 103, 7353, 1
999) transforms the energy "landscape" of biomolecules, making it possible to study parts of the PES that are more relevant to finding a global minimum. However, the deformed surface has too many "cathement".
, and even smoothing by the diffusion equation method does not guarantee isolating the lowest energy minimum in a multidimensional problem. Conformation space annealing (Lee et al., J. Comput. Chem.,
Vol. 18, p. 1222, 1997) is also limited to a small number of variables. The conformational annealing narrows the scope of the search from the entire conformation space to the low energy region, starting with a "pool" of minimized conformations. This minimized conformation is later modified by randomly picking variables from the "pool".

The Dead End Exclusion (DEE) method is based on identifying solutions that are totally incompatible with the global minimum (Desmet et al., Nature).
356, 539, 1992; Lasters et al., J. Am.
Prot. Chem. 16: 449, 1997). Solutions that cannot contribute to the local energy minimum above a certain order are excluded. The energy (cost) function must be described as the sum of terms, each of which is a function consisting of at most two variables. The value for the i-th variable x _i should not match the globally optimal solution if another value x ′ _i for the same variable is found.

[0049]

[Equation 2] By repeating the process, the solution can be sufficiently excluded to find the global minimum (Goldstein, Biophys. J. et al.
66, 1335, 1994). Another method is to combine individual search strategies with continuous minimization in the final stages of refinement (Dunblack and Karplus, Mol. Bi.
ol. 230, p. 543, 1993; Vasquez,
Biopolymers, 36, 53, 1995). The best-known application for this algorithm is the determination of the side chain conformation of proteins. If the DEE method cannot reach a unique structure, additional steps such as a forceful combinatorial search or an agglomerative approach are required for the remaining conformations (Becker, Protein).
s, vol. 27, p. 213, 1997). DEE faces significant practical problems. That is, while the above conditions require minimization for all possible variables, DEE can only minimize one per criterion. A further drawback is
The inability to find a group of low energy solutions.

Statistical methods (SM) use a model of the objective function to bias the selection of new sample points. These methods are justified by Bace's reasoning, which proposes that the particular objective function to be optimized is derived from a class of functions modeled by a particular probability function (Mockus, J. . Global
Optim. , Vol. 4, 347, 1994). Information from previous samples of the objective function can be used to estimate the parameters, and this refined model can be used to bias the selection of points in the search domain. The problem with using the statistical SM is whether the statistical model is appropriate for the problem. Moreover, mathematical complexity makes it difficult to write computer code for high-dimensional optimization problems. In many cases, SM relies on partitioning the search area into partitions, which limits these methods to problems of moderate dimensionality.

Unfortunately, none of the above-mentioned attempted solutions provide an adequate answer to the above-mentioned particular problem among the major problems of protein structure prediction, rather than to the search in combinatorial space. It only provides a general solution.

Thus, a solution to combinatorial in-space searching that is efficient, fast, easy to implement, and useful for various types of biological problems, such as those in the big problem of predicting protein structure. Is needed and it is useful to have such a solution.

SUMMARY OF THE INVENTION The present invention discloses a system and method for exploring combinatorial spaces without combinatorial explosion. At a more basic level, each combination is considered to consist of variables, each of which is at least 1
Can take one value. According to the invention, each variable is preferably one of a set of discrete values.
Takes one value, but each variable is, for example, a range of consecutive values or 1 in a function
It may take one value. These variables interact with each other in a manner known for individual interactions. Preferably individual interactions can be described for variable pairs such that the interactions are pair interactions. The search is preferably performed by sampling the value of one of each variable to obtain the combination. This process is then typically repeated many times. Each combination is evaluated by a quantitative measurement. The quantitative measure is preferably a cost function, and the desired result is typically maximized over the cost function during the process of determining which combination best meets the cost function. Or at least increased. For example, if the cost function is an energy minimization function, then a combination is selected that preferably has a low energy cost or value.

The present invention then provides for which factors do not contribute to the combination giving at least some minimum desired values for the quantitative measurements, and / or which are certain cutoffs or thresholds for the desired values. Try to determine which contributes to the combination giving a value for the quantitative measurement below. In other words, these elements are
It does not contribute to the "best" or most satisfactory combination for the system. These elements are then preferably removed or at least separated from the rest of the potential elements for making a combination.

The process of expelling the value of the variable is preferably repeated until only a predetermined number of combinations of elements that have not been removed and / or separated remain. At this point, an exhaustive search is best performed according to quantitative measurements and / or according to some other measure parameter.

The present invention accomplishes the above task by examining each value of the primitive at least once, and preferably multiple times, during the search process. Therefore, it can be said that each value is searched by the exhaustive search without the need to exhaustively search all combinations of values for the basic element. Therefore, the present invention has both the effect of non-exhaustive stochastic search processing and the effect of exhaustive search.

According to the present invention there is provided a method of searching through a combinatorial space characterized by a plurality of combinations, each combination consisting of at least one element, the steps of the method being carried out by a data processor, The method comprises: (a) providing a measurable quantitative parameter for each combination to determine whether the combinatorial space search result is successful, and (b) at least 1 each.
Dividing the combinations in the combination space into a population characterized by one combination; (c) calculating a value for the quantitative parameter for at least one combination of each population; and (d) the quantitative parameter. Determining the effect of each element on said value of, and (e) holding at least one combination according to said effect and providing the result of searching every corner of the combination space.

Hereinafter, the term “amino acid” shall refer to both natural and synthetic molecules capable of forming peptide bonds with other like molecules. The present invention will now be described by way of example with reference to the accompanying drawings. DESCRIPTION OF THE PREFERRED EMBODIMENTS The present invention discloses a system and method for exploring combinatorial spaces without causing combinatorial explosion. This search is carried out for various combinations of building blocks according to at least one desired property of the combination, which can be converted into a quantitative measure of the success of the search. The present invention then provides for quantification which factors do not contribute to the combination giving at least some minimum desired values for the quantitative measurements, and / or which are below a certain cutoff or threshold for the desired values. Attempts to determine which contributes to the combination giving the value for the static measurement. Preferably, elements are selected that contribute only to combinations that do not meet the minimum threshold for the desired value. In other words, these factors do not contribute to the "best" or most satisfactory combination for the system. These elements are then preferably removed or at least separated from the rest of the possible elements to make a combination.

The process of sorting the elements is preferably repeated until only a predetermined number of combinations of elements that have not been removed and / or separated remain. Such a predetermined number is optionally, or more preferably, an actual numerical value for the total number of combinations, or a minimum desired for quantitative measurements that a combination must meet in order to be included in the remaining combinations. It may be a threshold for the value. At this point, an exhaustive search is best performed according to quantitative measurements and / or according to some other measure parameter.

Due to the huge number of combinations, a sample of the combinations is preferably examined for the effect of the elements on the quantitative measurements for evaluating the combinations. Retain the elements of the combination with consistent maximization and / or promotion of quantitative measurements and exclude the other elements. The process may be repeated until a certain maximum number of combinations is found, after which these combinations may be evaluated according to similar parameters and / or some other parameter or characteristic. Such a group of combinations may optionally be observed as a population of combinations having a particular minimum for the quantitative measurement.

At a more basic level, each combination is considered to consist of variables, each of which has at least one value. According to the invention, each variable preferably takes on one value in a set of discrete values, but each variable may take on one value in a range of continuous values or functions, for example. These variables interact with each other in a manner known for individual interactions. Preferably individual interactions can be described for variable pairs such that the interactions are pair interactions. The quantitative measure of the combination of variables is preferably a cost function, and during the process of determining which combination best meets the cost function, the desired result is the overall maximum for the cost function. Or at least increased. For example, if the cost function is an energy minimization function, then a combination is selected that preferably has a low energy cost or value.

Various properties can optionally be used to create a cost function for performing the search in combinatorial space. For example, to search various protein structures according to amino acid sequences, the cost function can optionally be a combinatorial energy minimization, such that the selected structure represents at or near the energy minimum. Such energy cost functions are also useful for more specific or "subordinate" problems within the major problem of protein structure prediction. For example, minimizing the predicted positions of polar protons and side chains relative to amino acids provides quantitative parameters useful for these types of combinatorial exploration. However, it should be noted that in this case the maximization of the desired quantitative parameter is in fact achieved through minimization of the values of the energy calculation for the combination.

However, virtually any cost function can optionally be used with the method of the present invention. The cost function may not necessarily be associated with a biological problem, but may be associated with other types of problems, such as optimization of the cost function with respect to monetary value (literally economic “cost”). unknown.

The present invention accomplishes the above task by examining each value of the primitive at least once, and preferably multiple times during the search process. Therefore, it can be said that each value is searched by the exhaustive search without the need to exhaustively search all combinations of values for the basic element. Therefore, the present invention has both the effect of non-exhaustive stochastic search processing and the effect of exhaustive search.

As noted above, additional exhaustive searches may optionally be performed after practicing the present invention, for example, to identify absolute minimums as well as multiple local minimums. Such an additional exhaustive search is particularly preferred if the initial search process according to the invention comprises a stochastic search and / or comparison component, which is a preferred embodiment of the invention. The present invention is in many ways clearly distinguished from the background art search methods. First, the present invention is not based on any method known in the art, nor is it modified. Second, unlike other probabilistic search methods that cannot guarantee that every individual value in the combinatorial search space will be proved, the individual values of all variables in the combinatorial search space must be removed from the search space. You have to be certified. Third, the invention can also optionally and preferably obtain a population of local minimums in addition to the global minimum. The present invention is described in detail below, while providing the effectiveness of an exhaustive search, as evidenced below by comparison of the results of the present invention with the results of using a complete exhaustive search alone. , A probabilistic search can achieve the above objectives.

The principles and operation of the present invention may be better understood with reference to the drawings and the description set forth in several sections. The first part of the description (in this section) focuses on an exemplary general method according to the invention and a basic exemplary system for implementing the method. The following sections describe specific biological issues, and each section is labeled with its respective problem type. These sections describe examples for the proper practice and application of the invention and are not intended to give any limitation.

Referring now to the drawings, FIG. 1 is a flow chart of an exemplary but preferred general method for searching an entire combinatorial space in accordance with the present invention. As shown, in step 1, a combinatorial space is provided. Such combinatorial spaces feature multiple combinations of basic elements. The combinatorial space is optionally created, for example, by creating multiple structures with basic elements according to some pattern, plan and / or scheme. Alternatively, the combination space may optionally be pre-defined. For example, for biological problems, the combinatorial space may already be defined according to the type of biological structure to be analyzed.

In step 2, according to a preferred embodiment of the present invention, each combination is optionally and preferably constructed from variables each of which can take at least one value. According to the invention, each variable preferably takes on one value in a set of discrete values, but each variable may take on one value in a range of continuous values or functions, for example. These variables interact with each other in a manner known for individual interactions. Preferably individual interactions can be described for variable pairs such that the interactions are pair interactions.

In step 3, quantitative parameters are determined and according to the quantitative parameters it is measured whether the search was successful. Quantitative parameters must be measurable for each combination in the combination space. Quantitative parameters are typically calculated according to the basic elements of each combination, and the effect of structural features and / or interactions on this measure may optionally be further taken into account. For biological problems such as protein structure prediction, the types of quantitative parameters for investigating a particular problem may be known. For example, for the prediction of the position of polar protons within a protein, the best quantitative parameter is preferably the energy for the combination, determined according to equations well known in the art and described in more detail below with respect to Section 1. Is the minimization of.

The quantitative measure of the combination of variables is preferably a cost function, and in the process of determining which combination best meets the cost function, the desired result is generally the maximum for the cost function. Or at least increased. For example, if the cost function is an energy minimization function, the combination with the lower energy cost or value is preferably selected.

In step 4, the contribution of each element or variable is evaluated to determine the effect of the value of that particular element or variable of each combination on the quantitative parameter or cost function. Such effects are preferably determined both from the values of the variables and the interactions between the variables, as assessed by the cost function.

The preferred effect is for consistent maximization of the cost function. Consistent maximization is optionally measured according to the distribution of cost function values for all variable set combinations or large groups of "configurations". According to a preferred embodiment of the present invention, the effect of these different values is preferably determined by probabilistic analysis, especially with a large number of variables. This is because the efficiency of exhaustive analysis is considered to be orders of magnitude worse and time-consuming. The probabilistic analysis is preferably performed by randomly choosing values for each variable to create a combination, more preferably a plurality of different combinations. Most preferably, a predetermined number of such combinations are formed as part of the sampling process. The result or value of the cost function for each combination is then calculated according to both the value of the variables and the interaction between the variables.

In step 5, optionally and preferably, the values of the elements or variables that do not contribute to the consistent maximization of the desired result of the cost function are removed as described above. More preferably, the values of variables that contribute only less desirable results of the cost function, or results below a certain minimum threshold, but not above a certain threshold for the desired result, are eliminated. For example, for a cost function that includes energy minimization, a combination where the energy cost is found only in a combination above a certain threshold (less desirable result) and the energy cost falls below another low energy threshold (more desirable result). Values for variables not found in) are preferably removed. Alternatively, these values may optionally be "labeled" and / or separated for further analysis rather than being removed.

In step 6, if the total number of combinations has reached a certain minimum value,
These combinations are optionally and more preferably further analyzed according to a cost function and / or some other parameter to determine the results of the combinatorial search.
An exhaustive search may optionally be performed among the minimum number of combinations for the combination of interest, evaluated according to a cost function and / or some other parameter of interest. Such a set of combinations can optionally be observed as a population of combinations that have a particular minimum for the desired outcome of the quantitative measurement.

If not, preferably step 4 until the number of combinations reaches the above minimum number.
And 5 are repeated. “Minimum number” is an optional absolute number of combinations, and / or
Or, it refers to the minimum value for the cost function that such a combination must satisfy as a "cutoff" threshold.

FIG. 2 shows an exemplary system according to the invention for implementing the method of FIG. As shown, the system 10 features a computing device 12. In the present embodiment, the computing device 12 computes a number of functional modules so that the method of FIG. These functional modules are optionally and preferably implemented as software modules, but may alternatively be implemented as hardware, firmware or a combination thereof.

As shown, one such module is the combination storage module 14, which holds the combinations currently being considered in each population. The quantitative parameter calculation module 16 then calculates the value of the quantitative parameter for at least one combination in each population from the combination storage module 14. The evaluation module 18 creates a sample of multiple combinations from the elements of the combination, preferably such that one element is retained as consistently contributing to the maximized value for the quantitative parameter for that combination. , Evaluate the effect of each factor on the value of the quantitative parameter for that combination. These modules preferably interact until a certain minimum number of combinations is retained in the combination storage module 14, which represents the result of the search in the combination space.

The above description is a general description of the method and system of the present invention. The following sections describe specific model systems that are addressed by the present invention as specific problems, to which the present invention can provide solutions. These sections search the combinatorial space to identify polar proton positions (Section 1), identify amino acid side chains in proteins (Section 2), and identify loop structures in proteins. It includes predicting (section 3) and an explanation of other miscellaneous biological problems solved by the present invention (section 4).

Section 1: Polar Proton Localization The present invention is useful in solving the problem of accurately locating polar protons within biological molecules (eg, protein molecules or DNA). The localization of such polar protons then determines the position of hydrogen bonds within the biological molecule itself or between the biological molecule and another molecule. Therefore, this particular implementation of the present invention solves an important scientific problem.

The particular implementations of the invention described in “ Methods ” of this section were also tested against other methods known in the art as described in “ Results ”. It should be noted that these methods and results are presented for illustrative purposes only and are not intended to be limiting in any way. The interpretation of these results is then described in " Discussion ".

Method For testing purposes, the method of the present invention was implemented as a computer software program written in C ++. This program operates as shown in the flowchart of FIG. In step 1 as indicated, the program optionally reads the Protein Data Bank coordinate file format (PDB file) or receives input information from another source. The program uses an auxiliary ASCII file that acts as a database for parameterizing the atoms of the system. These files include connectivity of all atoms, atom charges, A and B parameters of the Lenard-Jones function, and the bond length between hydrogen and heavy atoms. The user can easily edit these files to add, delete, and modify residue types. These values are read from a file or input from another source to parameterize the atoms in step 2.

In step 3, the hydrogen to be added and the lone electron pair are divided into three categories: (1) trivial hydrogen-
Hydrogens that can be located using heavy atom coordinates and bond hybridization (eg, aliphatic and aromatic hydrogens), (2) nontrivial hydrogens (non trivia).
l hydrogen-polar hydrogens with rotatable degrees of freedom (eg, hydroxyls of serine, threonine, and tyrosine), (3) non-trivial lone pair of electrons-
Non-trivial hydrogen with the same geometric properties.

In step 4, trivial hydrogen is first added. These coordinates are calculated using heavy atom coordinates from the database, bond lengths and angles, and standard dihedral angles.

In step 5, the non-trivial hydrogen and the non-trivial lone pair of electrons (ens)
emble) and these coordinates are not calculated yet. An ensemble is defined as a collection of nontrivial hydrogen or nontrivial lone pairs of electrons that interact between themselves. The collective cutoff is user defined. The user can assign a large population cutoff value so that the system behaves as one large population. Collective cutoffs are measured from the coordinates of the heavy atoms bound to the non-trivial atom, since the position of the non-trivial atom has not yet been identified. The group is “segment (se
gment) ”. Each segment contains a rotation around a bond connecting two heavy atoms, one of which is bonded to a polar proton. Each segment can use various positions in the space that satisfy the H-bond condition.

In addition to the population cutoff, optionally and preferably two other cutoff conditions are used. The energy cutoff (the energy cutoff in the normal sense used in non-bonded energy calculations, ie the default is no cutoff) is used. Another cutoff is used to locate hydrogen-bonding partners around rotatable segments (see below) (this is
It may be larger or smaller than the "collective cutoff", but always> 3 to avoid the risk of losing the solution of the segment, including all adjacent partners for H-bonds.
Å). When this cutoff exceeds 4.5Å, many unrealistic non-essential partners are created, increasing the time to search for a solution. Collective cutoffs are used to create a collection of related heavy atoms (hydroxyl oxygen, water oxygen, NH ₃ ⁺ , amines, etc.) that must be solved for all members. Therefore, if this cutoff is 4Å, the distance of each pair of atom A and atom B or atom A and atom C may be less than 4Å, but R _{B, C} may be> 4Å, At the same time, all three atoms are part of the same population.

Preferably, each population is treated separately. In order to calculate the coordinates of the non-trivial hydrogen and the non-trivial lone pair of electrons, a two-dimensional matrix is formed in step 6. This two-dimensional matrix is a list of all hydrogen bonds that can be formed between the donor and acceptor. The larger the H-bond cutoff, the more
There are more choices of hydrogen bond linkages that are formed, and alternative interactions (alternat
The 2D matrix of ive interaction becomes larger.

As an example, the population shown in FIG. 4 contains only two carbonyls (1,2), one amide, and one hydroxyl, which together form one population. Hydroxyl donates one non-trivial hydrogen (3) and two non-trivial lone pair electrons (4,5), and amide donates one trivial hydrogen (6). A segment is defined as a collection of nontrivial hydrogen and nontrivial lone pairs of electrons attached to one heavy atom. For example, atom 3 and lone electron pairs 4 and 5 are one segment because they are linked to the same oxygen. The hydroxyl hydrogen (3) can form a hydrogen bond with any of the carbonyls, and the hydroxyl lone pair (4,5) is N
Assuming that hydrogen bonds can be formed with -H, the complete 2D matrix has the shape shown in Figure 5A. The two lone electron pairs are degenerate to form a hydrogen bond to the amide, so one of these can be omitted to form the first alternative combination of the 2D matrix ( 4-> 6 or 5-> 6). The omitted lone pair is automatically added after hydrogen and the first lone pair are located. Therefore, the initial 2D matrix has the shape shown in FIG. 5B. Module (mod
ul) refines the 2D matrix. That is, high energy value (
The positions that produce a salient value) are deleted. The energy threshold is user defined and a non-bonded energy equation is used.

A 3D matrix is formed in step 7 using the refined 2D matrix. Here, every combination in the population is uniquely defined (ie, in any combination, there is only one choice for any nontrivial (rotatable) hydrogen and nontrivial lone pair). . In this example, the 3D matrix has the shape shown in Figure 5C. Each pair of lines constitutes one contribution. Each combination is evaluated and the best combination is the population result. For more than one population, this process is repeated for each population.

The energy criterion used to assess the quality of each combination is the pairwise “unbonded” energy function:

[0090]

[Equation 3] Is. _Where A _{i, j} is a repulsion parameter of two (i, j) atoms, B _{i , j} is a polarizability attraction parameter of two (i, j) atoms, and q _i is a partial charge. , R _ij are distances between atoms. ε is the dielectric constant chosen to be 4 in testing the algorithm. This code is mutable and the force field can be easily modified in any way desired.

Energy calculations extend to the “boundaries” of each population. The cutoff distance for calculating E (r _{i, j} ) is user defined, but it is recommended to avoid cutoff so that long range electrostatic interactions can be explained. The main problem with this collective approach is to compute the interactions between non-trivial atoms in one population and non-trivial atoms in another. In this case, it is assumed that the coordinates of the non-trivial hydrogen in the second population whose configuration has not yet been determined match the coordinates of the heavy atom to which the non-trivial hydrogen is bound. This is a "unified atom" approximation that is tolerated by the relatively long distance between the known position of one atom and the as yet undetermined position of another atom. However, by forcing the program to treat the system as one huge ensemble (where all non-trivial hydrogen and non-trivial lone pairs are added simultaneously with exact positions), the user This approximation can be avoided.

It is clear that for large biological systems that make up one population, very large combinatorial problems arise. For example, RNase A (5RSA) has 1.76 ^* 10 ⁵⁹ alternative combinations for all rotatable hydrogens. Attempts to create 3D matrices from 2D matrices are beyond the capabilities of computers. A unique probabilistic approach was developed to reduce the size of this problem. This algorithm switches the exhaustive computation of the population to probabilistic computation if the number of combinations exceeds a user-defined threshold. The position in the d ₀ segment in this population is unknown. Usually, each non-trivial hydrogen or non-trivial lone pair of electrons has more than one position, but only one position exhibits the minimum energy. Nontrivial hydrogen and nontrivial lone pairs of electrons influence each other. That is, when the position of a lone pair of electrons is specified, the position determines the position of hydrogen bonded to the same heavy atom, and vice versa.

Let X = (X ₁ , X ₂ ... X _d0 ) be the configuration of the d ₀ segment in a population. For each configuration X, the energy E = E (X) is calculated according to the above energy function. The goal is to find the configuration that minimizes E. Since it is not possible to evaluate all alternative configurations for a large number of combinations, the following combinations are evaluated as an example for performing step 8 after the above steps.

1. N configurations from a large population of combinations, X ₁ = (X ₁₁ , X ₁₂ ... X _1d0 ). ． ． , X _n = (X _n1 , X _n2 ... X _nd0 ) are randomly sampled. Where X ₁₁ is the first randomly selected conformation of the first segment and X _n1 is the nth randomly selected conformation of this segment. FIG. 6A shows the first three configurations and the nth configuration sampled from the 2D matrix. Corresponding energy value: In case of configuration X ₁

[0095]

[Equation 4] In case of configuration X _n

[0096]

[Equation 5] To calculate.   2. Distribution Fⁿ _E(N = about 10^Three) Is created. Fⁿ _EIs for the whole protein
Is a set of energies corresponding to the n sampled configurations. Fⁿ _E Define cutoff points H and L at. H is E_i >Fⁿ _E(1-α)
Whereas all configurations that satisfyⁿ _E(Α) is Fⁿ _EThe first
α is the percentile), L is E_i <Fⁿ _EIncludes all configurations that satisfy (α).
The number of configurations in both H and L is n₀= N^*It is α. Maximum energy solid
N = 1000 configurations and α = for the configurations and minimum energy configurations
N if 1%₀= Α^*n = 0.01^*1000 = 10, so L
= 10 and H = 10. In other words, H is the maximum energy system of 10
In contrast (FIG. 6B), L represents the 10 systems with the lowest energies.

3. Create a vector h for a position in the configuration that matches the energy in H. The vector h is the element-wise intersection of all configurations in H as follows. That is, if all the configurations in H share the same value, then if the component j (corresponding to X _{nj of the} configuration X _n ) is 5-> 1, then h _j = 5-> 1, and otherwise If h _j = 0 (segment j in all high energy configurations
No common position about). For example, in FIG. 6B, all configurations of H share the same values 5-> 1 and 23-> 34, so these configurations are vector h = (5->1,23-> 34. ., 0). This is the vector created for the n ^* α high energy configuration of the d ₀ segment,
Value 5-> 1 in segment 1 and value 23-> 34 in segment 2
Appears in all high energy configurations (Fig. 6C). No common position was found for the last segment d ₀ in the high energy region.

4. Create a vector l for the position in the configuration that matches the energy in L. Vector l is the union of all configurations in L shown in FIG. 6D. Unlike vector h, more than one configuration may appear in each segment of l.

5. Compare h and l. vector component j with which both h _j and l _j are similar
, Which also contributes to the low energy value, is left as a viable configuration for this segment. However, if h _j ≠ l _j , the corresponding segment component h _j is excluded from later iterations. It should be noted that in a segment with a size equal to 1, h _j is the only available solution and is not excluded from later iterations. FIG. 6E shows that the value 5-> 1 is excluded from the further calculations because it exists only in the high energy vector h, whereas the value 23-> 34 in segment 2 is also excluded in the vector l. Is shown.
The new 2D matrix does not include the pair 5-> 1 as shown in Figure 6F. To avoid configurations with very high energies that can skew the results,
The number of configurations n and the percentile value α were chosen according to statistical equations that specifically deal with the probabilities of correct and illegitimate removal of configurations from a large set of combinations. Minimization of falsely excluded cases can be achieved by increasing α and n. However, the expected number of correctly excluded cases is also reduced, but the slope is smaller. Values of n = 500 and α = 0.008 were chosen as a good compromise (step 8 in FIG. 3).

6. Repeat steps 1-4 for the reduced position-space until the number of possible configurations is below a user-defined threshold (step 9 in Figure 3).
. 7. To find the best configuration, calculate E for all remaining configurations (exhaustive search; step 10 in FIG. 3).

Results The algorithm consists of five high-resolution crystal structures: Brookhaven Pro.
tein Data Bank (Bernstein et al., J. Mol. Biol
． 1997; 112: 535-542) file: Bovine pancreatic trypsin inhibitor (5PTI), RNAse-A (5RSA), trypsin (1NTP), and carbon monoxide-bound myoglobin (2MB5) (these neutron diffraction coordinates are related to the proton position). Available), as well as phosphate binding proteins (1IXH
) (This ultra-high resolution result is reported by X-ray).
The algorithm was run to remove all hydrogen atoms from the PDB file and reconstruct the positions of the hydrogen atoms, assuming that the hydrogen atoms were at optimal positions in the crystal.

Each system was treated with the following two variations of the method. 1. Combination "group-stochastic approach": Each system is divided into groups. Each population is treated separately. All possible combinations in the population are evaluated, the combination with the lowest energy is the result. In a population with very large combination requirements, a "probabilistic approach" was activated to reduce the number of combinations to a number that could be exhaustively evaluated. The advantage of this approach is that the CPU time required for the calculation is short. As an example, the calculation for 3INS and its aquifer by this method is Silicon Graphi
It is interactive on the cs R10000 machine and takes about 4 minutes. However, this approach requires an approximation of the distance between the non-trivial hydrogen and the non-trivial lone pair of electrons in different populations, as mentioned above, with some loss of accuracy.

2. Pure "stochastic approach": This program is forced to treat the system as one huge ensemble. This algorithm reduces the number of combinations to a number that can be exhaustively evaluated. All hydrogens in the protein are added simultaneously in each combination, so no approximation is applied during energy estimation. This is important if the energy difference between the few combinations is small and the accumulation of many long-range electrostatic interactions can add a significant contribution to the final result. This approach has a larger CPU demand. That is, the calculation for the same system takes about 15 minutes on a Silicon Graphics R10000 machine.

A test was devised to reveal whether a pure “stochastic approach” could find a global energy minimum from a large number of possible combinations, given the system and the energy function. To overcome the time limitation of this kind of calculation, a hypothetical protein was constructed. This protein has 1186 amino acids, as shown in Figure 7, of which 13 are serines (shown as a CPK model) (13 segments) and 1173 are glycines (0 segments). It has a spherical shape with a size of 64Å ^* 64Å ^* 61Å. The oxygen of the serine hydroxyl was placed at least 10Å apart. In this case, the interactions between the hydroxyls can be ignored and each segment can be treated as a separate population. All possible combinations in this population can be evaluated to obtain a global minimum of the system. Thus, the pure “stochastic approach” was compared to the collective approach (which in this particular case is roughly equivalent to a full exhaustive evaluation). The probabilistic search started with a total of 5.02 ^* 10 ¹⁰ combinations, reaching 2.7 ^* 10 ³ combinations after 204 iterations, and then 2.7 ^* 10 ³ combinations were exhaustively evaluated. The group method is 1 + 4 + 15 +
Only 12 + 10 + 11 + 2 + 10 + 6 + 12 + 5 + 8 + 11 = 107 calculations (sum of positions of all segments) were needed. The two methods gave the exact same results in terms of energy and proton position.

Five protein systems with high resolution coordinates from a combination of X-ray and neutron diffraction were analyzed. Only 5PTI, 5RSA, and 2MB5 have many water molecules (including proton positions) within the solvation shell.

Bovine Pancreatic Trypsin Inhibitor (5PTI, 1.8Å Resolution) The structure of trypsin inhibitor was determined by cooperative X-ray (1.0Å resolution) and neutron diffraction (1.8Å resolution) (Wlodawer et al., J. Mol. B
iol. 1987; 193: 145-156). This PDB file contains the coordinates of 58 amino acid residues and 63 water molecules. 2.5 containing 54 water molecules
A water layer of Å was included in this calculation. Potassium ions from PDB and PO ₄ ³⁻ ions were also included in the calculation. Atoms in the side chains of residues GLU7 and MET52 were found to occupy two major sites. The ^* A ^* format was selected for the calculation. We defined a group of rotatable atoms at a distance shorter than 4.5 Å as a group. There were a total of 21 populations and 256 possible positions.

A “combinatorial population-stochastic approach” was used. In Table I, the number of possible combinations in the population is the product of the number of combinations in each segment. The total number of combinations in group 9 is the result of multiplying the number of positions in each segment,
Therefore, 6,403,320 is reached. This group solved probabilistically, while other groups solved it exhaustively. "Total energy" is the sum of all populations.
The minimum energy of all combinations in separate populations is -121.0 Kcal /
The maximum energy was 2.9E + 16 Kcal / mo, while it was mole.
It was le. This maximum energy value is the result of a "protrusion" between rotatable hydrogens that could not be eliminated at this stage because only one proton protrusion was tested in the pretreatment stage.

The behavior of the system in the pure “stochastic approach” is shown in FIGS. 8 and 9a. FIG. 8 shows In (total number of possible combinations) vs. number of iterations. The initial number of combinations is 1.19 ^* 10 ³⁰ , of which only 2690 remain for exhaustive computation after 443 iterations.

FIG. 9a shows the energy distribution in the first and fourth iterations. The x-axis does not hold the same energy value for all iterations. That is, the average energy of the extracted sample decreases in progressive iterations. Therefore, the specimen is divided between 30 columns. That is, the minimum energy sample is in column 1 and the maximum energy sample is in column 30. The number of samples taken in all iterations is constant. It can be seen that the algorithm eliminates spikes in energy values during the iterative process.
Therefore, the energy distribution becomes more and more bell-shaped.

The structure of RNase-A (5RSA, 2.0Å resolving power) ribonuclease A was determined by cooperative X-ray and neutron diffraction (2.0Å resolving power) (Wlodawer et al., Acta. Crystallogr. B.
1986; 42: 379-387). This PDB file contains the coordinates of 124 amino acid residues, PO ₄ ^3- ions, and 128 water molecules. A 2.5 liter water layer containing 90 water molecules was incorporated into this calculation. 4 in 5 RSA during calculation
The four histidine residues were retained in the protonated form as seen in the PDB file.

A collection of rotatable atoms at a distance shorter than 4.5 Å was defined as one ensemble. A total of 37 populations and 485 possible positions (Table II) were received. "
The "combined population-stochastic approach" was used. Populations 2, 7, 10, 29 containing many combinations solved probabilistically, while the other populations solved exhaustively. The minimum energy, which is the sum of all minimum energy combinations in the separate populations, was -60.8 Kcal / mole, while the maximum energy was 261.3 Kcal / mole. Combinations with high energy values were excluded at an early stage by the pretreatment "protrusion value" calculation.

The behavior of the system in the pure “stochastic approach” is shown in FIGS. 8 and 9b. The initial number of combinations is 1.76 ^* 10 ⁵⁹ , of which only 2772 is 668.
It remains for exhaustive computation after iterations.

FIG. 9b shows the energy distribution at the 1 st and 4 th iterations. The energy distribution remains bell-shaped during minimization due to the lack of energy spikes. The structure of myoglobin (2MB5) myoglobin was determined by neutron diffraction (1.8Å resolution) (C
Heng and Schoenborn, Acta Crystallogr. B
1990; 46: 195-208). This PDB file contains the coordinates of 153 amino acid residues and 89 water molecules (including their protons). This contains protoporphyrin having Fe, ammonium ion, and sulfate ion. All water, ions, and protoporphyrin moieties were included in the calculation. HEM CO atoms are disordered. The ^* A ^* format was selected for the calculation.

As shown in Table III, a “combinational population-stochastic approach” was used. We defined a group of rotatable atoms at a distance shorter than 4.5 Å as a group. A total of 43 populations were obtained.

The behavior of the system in the pure “stochastic approach” is shown in FIGS. 8 and 9c. The initial number of combinations is 4.98 ^* 10 ⁵² , of which only 2400 have 552
It remains for exhaustive computation after iterations. FIG. 9c shows the energy distribution at the 1st and 4th iterations.

Trypsin (1NTP, 1.8Å resolution) The structure of trypsin was determined by neutron diffraction (1.8Å resolution) (Ko
ssiakoff, Basic Life Sci 1984; 27: 281-.
304). This enzyme was inhibited by the monoisopropylphosphoryl derivative, which was taken into account in the calculation. Calcium ion with 2+ charge PDB
Added according to file indications and placed near GLU70, ASN72, VAL75, and GLU80. This structure contains no water of crystallization. A group of rotatable atoms at a distance shorter than 4.5 Å is defined as a group.

Once again, the "combinatorial population-stochastic approach" was used. Table IV is 4
A total of 33 populations with a minimum energy of 83.9 Kcal / mole are shown. The behavior of the system in the pure "stochastic approach" is shown in Figures 6 and 7d. The initial number of combinations is 9.63 ^* 10 ¹⁰ , of which only 1152 remain for exhaustive computation after 14 iterations.

Phosphate-binding protein (1IXH, 0.98 Å resolution) The structure of the phosphate-binding protein was determined by X-ray diffraction (Wang et al., N.
at. Struct. Biol. 1997; 4: 519-522). The PDB file contains 321 amino acid residues. The coordinates of the water molecule have not been reported.
This protein forms a complex with PO ₄ phosphate, which has a −3 charge.
This ion was included in the calculation. This entry contains 6 irregular residues: Glu1.
, Ser3, Thr162, Pro216, Ser234, Lys245. The ^* A ^* format was chosen for all of these. As shown in Table V, the "combinational population-stochastic approach" was used. We defined a group of rotatable atoms at a distance shorter than 4.5 Å as a group. A total of 45 populations were obtained.

The behavior of the system in the pure “stochastic approach” is shown in FIG. The initial number of combinations is 1.18 ^* 10 ²¹ , of which only 2400 remain for exhaustive calculation after 51 iterations.

Discussion The five types of systems should be divided into two categories. The first category is systems that lack experimental data for water molecule coordinates. Such a system has trypsin (1
NTP) and phosphate binding protein (1IXH). Figure 10 shows a ribbon representation of 1NTP and its polar residues. Many polar hydrogens should form hydrogen bonds with water molecules. However, this PDB entry does not include water coordinates. In this case, the method of the invention lacks the data essential for the correct localization of the polar protons of the residues on the protein surface.

5RSA, 5PTI, and 2MB5 are systems with lots of experimental data on water location. These are the three most important systems in this study, and good algorithms are expected to yield accurate proton predictions for these systems.

The results of methods for locating protons in biomolecular structures should be evaluated by a few criteria. First, the quality of the results is related to the ultimate goal of achieving a negligible RMS of theoretical proton coordinates as compared to the previously described method as well as to experimentally obtained proton coordinates. Should be evaluated.

As shown in Table VI, the results of the “combinational population-stochastic approach” and the pure “stochastic approach” were compared to the experimental results, MSI Discov.
CVFF minimization using er / InsightII software package,
Compared to the method of Brunger and Karplus, and the method of Bass et al. CVFF minimization used a "steepest descent" algorithm in the first 100 iterations, followed by a conjugate gradient method until convergence with a maximum derivative of less than 0.001 Kcal / Å was achieved.

A standard program (for example, Discover / CVFF (BIOSYM /
Molecular Simulations. Discover 2.9.7
Force field simulations user guide
1995; Part 1; BIOSYM / Molecular Simulati
ons. InsightII 95.0 Molecular Modelin
Insight (BIOSYM / Molecular Simulatio) with further optimization according to g System User Guide; 1995.
ns. Discover 2.9.7 Force field simula
Tions user guide 1995; Part 1; BIOSYM / M
olecular Simulations. InsightII 95.0
Molecular Modeling System User Guide
1995)) compared to proton localization according to the invention, the improvement by the methods of the invention (population-stochastic and "pure stochastic" methods) was clearly demonstrated.
Brunger and Karplus (Proteins 1988; 4:14).
Self-consistency algorithms, such as the algorithm of 8-156), usually yield better results than non-specific methods. However, these are not as accurate as the method of the invention. The method of the invention yields better results than Bass et al. In the more experimentally accurate 5RSA and 5PTI (Ser, try, and water R in Table VI).
See MS values), giving similar results in the less accurate System 1 NTP.

The present invention has two additional improvements over Bass et al. First, 9
Since a system with as many as two segments (5RSA) can be treated as one huge ensemble ("pure stochastic" approach), there is no restriction on ensemble size, whereas Bass et al. (Proteins). 1992; 12: 266-.
277) is limited to a very small size. From Tables IV, close distances between rotatable protons in some regions of the protein should be considered when considering together with the number of options for localization of each proton. It is clear that some cases require extremely long calculations. The probabilistic method is better equipped to handle fairly large molecules, since special attention must be paid to the need to locate protons in large molecules in a relatively short time. Second, the method of the present invention is energy-based, whereas the method of Bass et al. (Proteins 1992; 12: 266-277) is not energy-based.

No consistency was found with regard to the improved prediction of one residue type compared to the other. This also applies to other methods of locating protons. However, in some cases, we find a correlation between the order of our RMS results and the order of our Bass RMS results for residue types. This may be related to the spatial distribution of these residue types in each protein (ie, some residues are closer to the protein core and some are closer to the protein surface) and May not be accurate due to lack of information about.

It is expected that the pure stochastic approach will drop some low energy solutions during the iterative calculation that excludes the solution, and thus yield inaccurate results. It is noteworthy to find out how well a pure stochastic approach works, as compared to ensemble-stochastic methods that solve most of the system exhaustively. The "virtual protein" mentioned in the results section was used to compare the "pure" stochastic approach with the exhaustive search. Both approaches yield the same minimum from 5.02 ^* 10 ¹⁰ possible combinations. This is a complementary suggestion of the certainty of a "pure" stochastic approach as a tool to find a global minimum.

Information on the characteristics of the system can be obtained by examining the energy distribution table in detail (FIG. 9). The bell-shaped distribution in the first iteration shows no protrusion between the rotatable hydrogens. The "ordered" bell shape of the energy distribution for the positions of rotatable protons obtained after a few iterations could be an indication of protein density near these protons. That is, "crowded"
The protein should increase the barrier to rotation. Therefore, its energy should be biased to the upper end of the energy spectrum. The bell shape may be evidence of the relatively "free rotation" of these protons in the less dense surroundings.

Section 2: Localization of Amino Acid Side Chains The present invention is also particularly useful in solving the problem of correctly locating amino acid side chains within proteins. This particular implementation of the invention solves the difficult problem by allowing such a position determination with considerable accuracy without undue assumptions and combinatorial explosions.

Certain implementations of the invention described in “ Methods ” of this section were also tested against other methods known in the art as described in “ Results ”. It should be noted that these methods and results are presented for illustrative purposes only and are not intended to be limiting in any way. The interpretation of these results is then described in " Discussion ".

Method Exploration Approach The code is based on a backbone-dependent rotamer library (backbone).
Dependent Rotor Library (Bowe)
r et al. Mol. Biol. 1997; 267: 1268-1282; Dun.
black and Karplus, Nat. Struct. Biol. 1994
1: 334-340; Dunblack and Karplus, J .; Mol.
Biol. 1993; 230: 543-574). For testing purposes only, and not in any way limiting, the August 1997 updated version of the Dunblock and Karplus rotamer library was used in the tests described below. Use the fused atom model (Weiner et al., J Amer. Chem. Soc.
1984; 106: 765-784). Energy is AMBER unbound 12-
6 Calculated according to Eq. 1 using the Lennard-Jones and electrostatic energy terms. _Where A _{i, j} is the repulsion parameter of two (i, j) atoms,
B _{i, j} is a polarizability attraction parameter of two (i, j) atoms, q _i is a partial charge, r _ij is a distance between atoms, and ε is a dielectric constant. A distance dependent permittivity ε = r is used. V _n is the twist potential barrier height of the twist angle φ, n is the multiplicity, and γ is the phase factor.
The potential of V _n is chosen from the AMBER force field parameters. The non-bonding energies for the interactions of the backbone and other residues with rotamers are calculated. Torsional energy terms are calculated for all dihedral angles of the rotamer of each residue. For a particular atom pair, the non-bonding energy term is 10 Kcal /
If it exceeds the value of mole, it is rounded down to 10 Kcal / mole.

[0132]

[Equation 6] Bower et al. (J. Mol. Biol. 1997; 267: 1268-128.
As suggested by 2) and implemented in the SCWRL algorithm, every rotamer is given a local energy based on its probability in a backbone-dependent rotamer library. _Energy, as _{-ln (p rotamer / p 0)} , obtained from the probability of the backbone-dependent rotamer library. Wherein, _{p 0} is the probability of the most probable _{rotamers, p rotame r} is (assuming kT = 1) is a certain probability of a particular rotamer. This search strategy involves several steps:

(I) Steric collision elimination step and preliminary rotamer localization : The input for calculation is the backbone of a protein with a known structure (N, C _α , C, O).
Coordinates. These coordinates and the φ and ψ angles of the backbone are used to create the initial configuration of the possible rotamers of each residue. Possible disulfide bonds between cysteine residues are calculated by the distance between the sulfur atoms. All rotamers that collide with the backbone are excluded. When all rotamers of a residue collide with the backbone, the rotamer with the lowest "collision energy" is left behind. This algorithm treats multiple single rotamers as part of the backbone (ie, other rotamers that conflict with these residues are also excluded).
. The algorithm also uses rotamer i of amino acid j and rotamer k of amino acid l.
Search for all side chain collisions with. The algorithm excludes such pairs from the solution, so such pairs are not sampled at the stochastic stage (see below).

[0134]   (II) Probabilistic stage: Very large for large biological systems such as proteins
It is clear that there will be a problem of proper combination. For example, hydrolase (1
arb) (Tsunasakawa et al., J. Biol. Chem. 1989; 264.
: 3832-3839), after step I, 2.29.^*10¹⁰⁵Alternative to
There are placement options. In order to reduce the size of the problem, a new stochastic algorithm
Is used. In this protein, d₀Side chain rotation difference in amino acids
The sex is unknown. Usually more than one rotamer for each amino acid
However, only one rotamer shows the minimum energy. X_j= (X _j1 , X_j2． ． ． X_jd0) Is d in the protein₀Amino acids randomly
It is the conformation of the protein containing the selected rotamer. Each
Conformation X_jAbout energy E_j= E (X_j) Is the energy
It can be calculated according to the Ruggie function. The goal is a conformation that minimizes E.
To find a solution. All alternative conformations for numerous combinations
Since it is not possible to evaluate the formation, the following steps are taken.

1. From the population of large combinations, n conformations, X ₁ = (X _{1 1} , X ₁₂ ... X _1d0 ) ,. ． ． , X _n = (X _n1 , X _n2 ... X _nd0 ) are randomly sampled. Where X ₁₁ is a randomly selected rotamer of the 1 st amino acid in the 1 st conformation and X _n1 is a randomly selected rotation of the same amino acid in the n th conformation. Is an isomer. We used n = 1000 to create many sufficient numbers of protein conformations, and the corresponding energy values: E ₁ = E (X ₁ ) ˜.
Calculate E _n = E (X _n ).

2. Create a distribution F ⁿ _E (n = approximately 10 ³ ). F ⁿ _E is the set of n samples extracted conformational all the energy for the entire protein. Defining a cut-off point H and L in F ⁿ _E. H is E _i > F ⁿ _E (1
While all variable values that satisfy −α) are included (wherein F ⁿ _E (α) is the ^αth _percentile of F ⁿ _E ), L is E _i < F ⁿ _E ( Includes all variable values that satisfy α). The number of conformations in both H and L is n ₀ = n ^* α. If n = 1000 conformation and α = 0.01 (1%) for maximum and minimum energy conformations, then n ₀ = α ^* n = 0.01 ^* 1000 = 10, so L = 10 and H = 10
Is. In other words, H represents the 10 highest energy conformations, while L represents the 10 highest energy conformations.

3. Create a vector h for all rotamer variables that match the conformation in H. Vector h is the elementwise intersection of all rotamer states in H, as follows: That is, if all rotamer states in H share the same rotamer in component j (corresponding to X _{nj of} conformation X _n ), then h _j = rotamer Number, otherwise h _j = 0 (no common rotamer for j in all high energy conformations).

4. Create a vector l for rotamer variables that match the conformation in L. Unlike the vector h, for each amino acid j more than one rotamer may appear up to the maximum of the n ₀ values in l _j . It is the union of all rotamers of component j that appear in the low energy conformation of L.

5. Compare h and l. If both h _j and l _j have similar rotamers, this also contributes to the low energy value and is therefore left as a viable rotamer state. However, if h _j does not match any element of l _j , then
The corresponding rotamer h _j is excluded from later iterations. If an amino acid has only one rotamer, this rotamer is the only remaining solution and is not excluded from later iterations.

6. Repeat steps 1-4 for a reduced set of variable values until the number of possible combinations of all variables is less than the user defined "probabilistic step criterion limit".

The value of α used to determine n ₀ should be chosen carefully. If α is too large, rotamers are not excluded. If α is too small, unreasonable exclusion of rotamers can occur. At best, α should be adjusted by the number of possible rotamers of each amino acid to give the same probability for elimination of rotamers. To explain the determination of α, it is assumed that no rotamer is affected by its interaction with any other amino acid in its environment.
2 to 2 of one residue that correctly excludes rotamers with> 99.983% certainty
The α values for the 9 possible rotamers are shown in FIG. These values were calculated as follows. Given a residue with three rotamers, if we want to exclude one rotamer with greater than 99.99% certainty ( _Pcorrect ), the probability of _error ( _Perror ) is 0. It must be less than 01% (0.0001). In order to remove the rotamer by mistake, this rotamer must first appear in all high energy conformations. In this case, the probability is (1 /
3) ^α . Furthermore, this rotamer must not appear in any low energy energy conformation. In this case, the probability is (2/3) ^α .
All error probabilities are P _error = (1/3) ^α (2/3) ^α . Therefore,
By using the general formula P _error = (1 / m) ^α (m−1 / m) ^α ,
The calculation can be adjusted to almost 100% confidence. In the formula, m is the number of variable values (rotomers). If m = 1 (one rotamer), then P _error = 0. By assigning a value of P _error = 0.0001 and solving the equation, α = 6
． A value of 12 is obtained. When α is very large, P _error = 0, but the probability of excluding any variable value is very low. Therefore, P _correct = 99.983
% Values that allow the removal of variable values between% -99.99988% are preferably in accordance with FIG.
Used from.

(III) End of search : When the number of remaining combinations is less than M (M is about 1
0 ⁵ ), perform an exhaustive search to obtain the N minimum energy conformations of the protein.

[0143]result   Stochastic algorithms cover a wide range of protein fold families
Chosen to vary in size (46-263 residues) and complexity (back
1.04 after elimination of rotamers that collide with Kuborne^*10¹⁴~ 2.29^*10 ¹⁰⁵ (Possible combinations of) and 10 proteins. These one
6 out of 0 proteins (46-48 residues) are DEE / A^*Algorithm
And Lemon (Proteins 1988; 33:
227-239), the stochastic algorithm and DEE / A^*A
Useful for comparing rugorism. These proteins are:
Crambin (PDB entry 1crn) (Teeter et al.,
J Mol Biol. 1993; 230: 292-311), ribosomal tan.
Quality (1ctf) (Leijonmarck and Liljas, J Mol)
  Biol. 1987; 195: 555-579), complement regulatory protein (co
implement control protein) (1hcc) (Norm
an et al., J Mol Biol. 1991; 219: 717-725), Obom.
Coid third domain (2ovo) (Empie and Laskowski, Bi
Chemistry 1982; 21: 2274-2284), erabutoxy.
B (3ebx) (Smith et al., Acta Crystallogr A.1).
988; 44: 357-368), and rubredoxin (5rxn) (Wat
empauugh et al., J Mol Biol. 1980; 138: 615-633.
). The remaining proteins selected were larger (residues 129-263) and
Resolvable X-ray structure (resolution <1.5Å, R factor <0.17): lysozyme
(2ihl), ribosomal protein (1whi) (Davies et al., Stru.
Culture 1996; 4: 55-66), endonuclease (2end) (
Morikawa et al., Science 1992; 256: 523-526),
And hydrolase (1 arb) (Tsunasakawa et al., J. Biol. Ch.
em. 1989; 264: 3832-3839). Table VII shows a stochastic algorithm
The results of applying the rhythm to 10 kinds of proteins are summarized. For each protein
About the combination (after the initial elimination of rotamers that collide with the backbone)
Number and some values of multiple single conformations (global for each protein)
Luminum) and 1000 low energy conformations
The average value of "population" is shown. The best possible RMS for each protein is shown. Finally,
Shows the average energy gap (without weighting) of 1000 conformations
. Global energy minimum conformation compared to X-ray conformation
Side chain atom (C_βRMS) was calculated. Global minimum
The RMS range for is 1.32 to 2.60. 1000 low energy
The average RMS value for the conformation is somewhat higher than the global minimum, but
Energy quality is higher than the global minimum and higher than the minimum.
A conformation with a small RMS has been found. 1000 minimum
The energy value range of the energy conformation is 5.52 Kc from the global minimum.
al / mole There is a large value. 1000 minimums from the global minimum
The average energy gap of the energy conformation is always small (for all proteins
2.20 Kcal / mole).

Testing the Validity of the Search Method To test the effectiveness of the probabilistic search, and in view of the values reported in Table VII, a number of issues were raised. The first question is whether, given a particular rotamer library, the results obtained by the exhaustive search can be obtained by the stochastic search. The second question is whether such a search will identify the crystal structure of the protein if the rotamer library contains the original X-ray rotamers.

The first problem requires the testing of relatively small proteins that can carry out such an exhaustive search. Given the constraints of the energy function and rotamer library, our probabilistic algorithm finds the lowest energy combinations in the test proteins and compares the lowest energy combinations with the results of the exhaustive search. Was imposed on. The selected proteins are clambin (PL format), which is a high-quality X-ray structure (1.05Å resolution, R factor = 0.105), Brookhaven Protein Data Bank (Bern).
Stein et al. Mol. Biol. 1997; 112: 535-542) file 1 cnr (Teeter et al., J Mol Biol. 1993; 230:
292-311). Crambin is large enough to be a reliable test case, but not very large and not large enough to require long operations in an exhaustive search. . The entry contains 46 amino acid residues (see Figure 11) and the coordinates of the ethanol molecule. Eight irregular residues (Thr1, Thr2,
Ile7, Val8, Arg10, Asn12, Ile34, Thr39). To evaluate this protein at the appropriate time, Arg10 (A form for this irregular residue), Arg17, Glu23, Ile33, and Ile35 were fixed in the original position. The initial number of combinations (after the step of eliminating steric collisions) was 6.79 ^* 10 ⁸ . In FIG. 12, the results of the probabilistic search and the exhaustive search for a wide range of N low energy conformations are compared. The energy values for each of the 10,000 conformations and the% difference between the two searches are shown. The 485 minimal energy conformations of this protein were found to be exactly the same by stochastic and exhaustive searches. The difference is small for many conformations. The energy of conformation number 10,000 is 4. From the global minimum in the exhaustive search.
71 Kcal / mole large, 4.80 Kcal / mole in stochastic search
It turns out to be big. Therefore, the difference between these searches for this conformation is 0.56
There is only%. This test consists of three different genus (10000
0, 200000, and 300000) with similar results.

To test the ability to reproduce X-ray coordinates, a stochastic algorithm was used with an expanded rotamer library (introduced 1 cnr of crystalline rotamers). No residues were fixed during this search. The energies were calculated by Eq. 1 without the probability term (probability term not available for crystal coordinates). The following residues: 4
1 Gly (no side chain), 5 Ala (only one possible rotamer),
And 6 Cys, which are rotamers-all of which form S—S bonds, were not included. Therefore, from amino acid 46 to 31 in the sequence were left for this comparison. The energy of the protein in this crystal structure coordinate is 3.4 from the global minimum found by the stochastic algorithm.
It was 1 Kcal / mole higher. 20% of the residues (6 amino acids: Ser6,
In Val8, Thr21, Tyr29, Ile33, Ile34), the search results were perfectly superimposed on the X-ray results. High quality overlays were found in 58% of the residues (18 amino acids). That is, it has been found that the absolute angular deviation of all twist angles is less than 40 °. Thus, the expanded rotamer library (this RMS should be 0.0) pinpointed approximately 80% of the side chains. For example, the atomic CGs of Leu 18 are 0.18 apart,
The CG atoms of Asn14 were separated by 0.23Å. RMS values between the global minimum structure and crystal structure of the side chain atoms (excluding C _beta) was 1.16.

To test the limits of the original rotamer library (no crystallographic rotamers), each rotamer was placed as close as possible to the relevant side chain in the crystal structure. The RMS value obtained was 1.15. The RMS value between the global energy minimum and the crystal structure in the stochastic search is 1.97.
Was found.

Comparison of algorithm with results from X-ray, NMR and MD Results for 1000 lowest energy conformations can be taken for each side chain under different conditions: X-ray crystallography, NMR and MD. To compare with the experimental and theoretical methods to detect conformation, E. The conformational space of E. coli ribonuclease HI was probed using the method of the invention.

A set of 8 NMR structures (PDB entry 1rch) was reported based on distance constraint from the experiment. Phi
lippopoulos and Lim (Proteins 1999; 36: 8).
7-110) provides an expanded set of NMR results with high resolution (2rn2,1.4).
8Å) Crystal structure (Katayanagi et al., J Mol Biol. 1992;
223: 1029-1052), low resolution (1 rnh, 2.05Å) crystal structure (
Yang et al., Science 1990; 249: 1398-1401), and these MD simulations. NMR and MD simulations yielded very few twist angle results for each, and the resulting conformations were classified as a population. Each population has an average dihedral angle and an order parameter S (Hyberts et al., Protein Sci. 19).
92; 1: 736-751) (representing the deviation of each dihedral angle from its mean value). The S-parameter of each dihedral angle at each residue is then calculated over the population. Order parameter S of angle α _i of residue i
(Α _i ) (where α = φ, ψ, χ ₁ , χ _2, etc.) is

[0150]

[Equation 7] Is defined as Where N is the total number of structures in the population and α ₁ ^j (j = 1
,. ． ． , N) is a 2D unit vector with a phase equal to the dihedral angle α _i , where i
Represents the residue number, and j represents the population number. If the angles are the same in all structures, S has a value of 1, but values of S much smaller than 1 indicate irregular regions of the structure. Philippopoulos and Lim classify these categories as 0
． Limited to S-values greater than 8.

[0151]   The stochastic algorithm is a 2 rn2 backbone with higher X-ray accuracy.
Used against. Calculated as 1.61^*10⁸⁷Started with the possible combinations of
. The conformational space is reduced to 1.3 by removing the high energy conformation.^*10 ³³ Algorithm is used to refine the best conformation of
And remains in a conformation sufficient to assess the conformational flexibility of the protein.
You

Table VIII contains a comparison of X-ray crystallography, NMR, and MD results with stochastic algorithms. This table focuses on the highly likely conformational residues according to the following assumptions. In some cases, the twisted angle has a single conformation in the MD population, and in the NMR population, a plurality of conformations, and in some cases, the opposite is obtained. The present inventors show that the rotamer obtained from the experiment appears in the following rules: (1) high resolution crystal structure (2rn2). (
2) The following three: low resolution crystal structure (1 rnh), NMR model, and MD
According to one or more of the found in at least two of the simulations,
Assume a high probability of rotamers obtained from the experiment. "Hit" is "correct"
It was considered to be any result of a stochastic algorithm with a variation from conformation to ± 30 °. Each such hit is marked with a "+" in the table.
In some cases, angles such as χ1 of M47 are indicated in the table by one rotamer and are marked with “(+)”. Such angles have additional values that do not follow the above two rules. Other angles are considered to have low probabilities and are listed in Table VII.
Not shown in I. Seven of the 115 dihedral angles in Table VIII are missing from the rotamer library (see FIG. 13A) and the other two are offset by about 40 °, thus a “hit”. Was not included in the evaluation by the inventors. Therefore, we can obtain up to 106 compared to X-ray, NMR and MD.
Can be expected to be a "hit". The probabilistic algorithm correctly predicts 87 angles (82%) (see Figure 13B).

[0153]Comparison between algorithm and DEE / A ^* algorithm   Leach and Lemon (Proteins 1998; 33: 227-.
239) selected to cover a wide range of protein fold families
DEE / A for each set of 8 proteins^*Algorithm
I searched the space of formation. Then, six of these proteins (1c
rn, 1ctf, 1hcc, 2ovo, 3ebx, 5rxn)
Method was used. Snake Venom Neurotoxin (1nxb) (Tsernoglou et al., M
ol Pharmacol. 1978; 14: 710-716) is an unknown residue
Excluded due to ip (residue 59). Bovine pancreatic trypsin inhibitor (5pt
i) (Wloderer et al., J Mol. Biol. 1987; 193: 145.
-156) is the occupancy of the two major sites by residues Glu7 and Met52.
Was excluded because of. Leach and Lemon also refer to "standard" and "
The effects of "fused" and "whole atom" models with "decreased" electrostatic display
I searched. Unfortunately, they do not report the RMS value for each system separately,
Only the average value of all eight systems of Table IX shows the stochastic method and DEE / A ^* Includes comparison with. Of combinations solved by a stochastic algorithm
The maximum number is 2.29^*10¹⁰⁵While DEE / A^*Is 2.48^*1
0³⁴I only reached the combination of. Solved by a stochastic algorithm
The largest protein system is 263 amino acids, while DEE / A^*Is the maximum
Solved 68 residues. The predicted consonants are then used to assess model accuracy.
Side-chain atoms (C) in conformation and X-ray conformation_β(Except)
The average RMS was calculated. Best possible RMS of current rotamer library
Indicates. DEE / A^*In the same system also calculated by, the stochastic algorithm
The RMS value range of the rhythm is 1.32 to 2.48, and the average is 2.07.
Was found. The best possible RMS for our rotamer library is
, 1.18 for all proteins in our test case. Lea
ch and Lemon 1 depending on atom model and rotamer library
． 77-1.92 average RMS value and 0.75-0.83 rotamer library
Reported Lee's best possible RMS. DEE / A for combinatorial explosion ^* In a larger system that could not be used in a stochastic algorithm
Found an average RMS in the range of 2.22 to 2.60 and an average of 2.40. Different rotation
The best possible RMS for the sex library was 1.23.

Discussion The above discussion relates to the application of a novel probabilistic search approach to search the conformational space of protein side chains. This is a development and improvement of the previous example of the previous section for searching the position of polar protons in proteins. This algorithm can successfully search a large number of conformational spaces of a protein and process many combinations after eliminating rotamers that collide with the backbone.

The certainty of probabilistic algorithms in the processing of complex combinatorial searches is shown in Table VI.
Proven in I and IX. Comparison of probabilistic search and exhaustive search (
Figure 12) proves the reliability of the probabilistic algorithm in finding the increased amount of the minimum energy conformation. No difference was found between the stochastic and exhaustive searches for the 485 low energy conformation of this protein. When approaching the limit of 10,000 minimum energy conformations, a slight deviation of 0.56% was detected. This population contains a major contributor to the molecular properties according to the Boltzmann distribution, as this large number of conformations reach an energy gap that is up to 4.71 Kcal / mole greater than the global minimum. These are the major contributions to the molecular partition function that can be used to compute the conformational entropy.

There is no difference between the probabilistic search and the exhaustive search for low energy conformation populations, and another issue is the comparison between our search results and experiments.
This is detailed in Table VII by comparing our global minimum and crystallographic results for 10 proteins and detailed with our low energy population for one protein. It is shown in Table VIII by comparing the X-ray, NMR and MD results. Global minimum RMS values for 10 proteins are strongly influenced by rotamer libraries, but not completely by rotamer libraries (ie,
The energy formula is limited to the reproduction of the structure). The inclusion of the original crystallographic rotamer in this library proved this point. Even in this case, only about 80% of the residues were calculated with high accuracy. The energy of the 1 cnr crystal coordinate was found to be 3.41 Kcal / mole greater than the global minimum found by our probabilistic algorithm. However, the RMS value of this structure is only 1.16. 1cn without crystallographic rotamers
A stochastic search at r yields an RMS of 1.97. The limits imposed by the rotamer library are shown in Table VII for each protein by the column showing the best possible RMS of this library. These values result in 50-75% of the errors in the RMS value of the global energy conformation. From Table VII it can be seen that the global minimum energy conformation does not necessarily have the minimum RMS value. There is a higher energy conformation and its structure is closer to the X-ray result.

Table VIII shows the E. coli expected to be detected by comparison with X-ray, NMR, or MD. The E. coli ribonuclease HI contains 106 angles. The algorithm correctly detected 87 angles, which is 82% of the total. Some of the deviation from 100% accuracy may be due to the quality of the rotamer library, but most is due to the energy function. Mendes et al. (Proteins 19
99; 37: 530-543), a continuous population of conformations in which one rotamer is grouped around a classical fixed rotamer.
Ensemble). Such an approach can enhance the effectiveness of rotamer libraries. These results (RMS of 1.16 in the presence of the crystallographic rotamer, 1.97 in the absence of the crystallographic rotamer)
Supporting the claim that the larger rotamer library does not guarantee a dramatic improvement in RMS values (Proteins 1992; 14: 213-223; J
． Mol. Biol. 1994; 235: 1088-1097; Tanimur.
a et al., Protein Sci. 1994; 3: 2358-2365; Vasq.
Uez, Biopolymers 1995; 36: 53-70).

Currently, there are four ways to study the conformational space of a particular protein.
There are two main methods: methods based on X-ray crystallography, NMR, MD, and rotamer libraries. X-ray crystallography usually suggests a single structure that may be biased to a particular conformational metastable state in the crystal (Brunger).
, Nat. Struct. Biol. 1997; 4 Addendum: 862-865). Observation of different conformations may only be possible with the highest resolution. The advantage of our algorithm is simple and unambiguous: it extends from a single conformation to a family of viable conformations.

Unlike X-ray crystallography, NMR suggests alternative conformations by decoding the 2D and 3D coupling maps. From NMR, the shape of the minimum energy value on the potential energy surface is unknown. Protein N
MR is a long and tedious experiment limited by the timescale of conformational changes, especially for large proteins. In this case, the method of the present invention may be an additional tool for suggesting alternative conformations. NMR
Allows the determination of the amount of conformational energy, if a structure is available,
Therefore, the method of the invention can be used to expand this information by allowing the assessment of the contribution of the energy content to the total population at equilibrium.

MD simulations require a large CPU timescale for biomolecules, which prevents a perfect search of the conformational space. MD is NM
It suggests a conformation that cannot be detected by either R or X-ray crystallography. The time scale and barrier-passing ability of MD are not yet reliable enough to detect global minimum or minimum energy conformation populations in large biomolecules. The reliability of our probabilistic algorithm in both of these findings is proved in this document. However, M
The D trajectory shows a conformational interconversion mechanism, whereas the stochastic approach focuses on products rather than pathways.

Dill and Chan (Nature Struct. Biol. 1997.
4: 10-19; Chan and Dill, Proteins 1998,3;
0, 2-33), the native state of a particular protein is consistent with the global minimum in free energy (not necessarily the global minimum potential energy). Therefore, the algorithm that adds side chains should yield most of the minimum energy conformation to allow entropy estimation. Currently, the method of the present invention meets this need.
In the “mean field” approximation, each side chain is “affected” by the average of all possible conformations of the side chains next to it. The conformational entropy is then evaluated from side chain probabilities at certain possible positions (Vesquez, Biopo.
lymers 1995; 36: 53-70; Koehl and Delarue.
, Nat. Struct. Biol. 1995; 2: 163-170; Koeh.
1 and Delarue, Curr. Opin. Struct. Biol. 19
96: 222-226). The probabilistic search of the present invention is unique in that, in addition to finding the global minimum, it provides the N closest best solutions for rotamers in large proteins without any mean field approximation . Therefore,
The probabilistic search of the present invention can be used to study the thermodynamic properties of complex molecular systems. The probabilistic algorithm can handle over 250 residues (maximum value at this stage is 2.29 ¹⁰⁵ combinations). The DEE / A ^* algorithm processed up to 68 residues and the maximum number of combinations (before backbone collision exclusion) was 10 ⁴⁴ . After applying the DEE algorithm, the size of the remaining space searched by the A ^* algorithm can be reduced up to 10 ²¹ .

The quality of the method of the invention (using the energy scheme of the method of the invention and the backbone-dependent rotamer library) was compared to the combination of DEE / A (using different energy schemes and two different libraries). ^* Algorithm (Leach and L
Emon, Proteins 1998; 33: 227-239). Since RMS is affected by rotamer libraries, comparisons between RMS approaches and experiments are limited. That is, an RMS value of 2.0 with a rotamer library where the protein has a minimum RMS value of 1.9 has an RMS value of 1.5 obtained from the library with an optimal RMS of 0.1. It shows a better search approach. The RMS value should be compared to the optimum RMS value that can be achieved within the framework of the rotamer library. It can be seen in Table IX that there is a correlation between the best possible RMS value of the library and the RMS of the global energy conformation. This could explain the difference between these results and the results of Leach and Lemon. Another advantage is that the probabilistic algorithm can be used in an "independent" form without any pre-processing algorithm (eg DEE in the case of the A ^* algorithm). The A ^* algorithm requires a good estimate of cost dependence to reach the goal node. This can be difficult to achieve. Because interactions between residues that have not yet been assigned to any position cannot be easily calculated. It should also be noted that the number of combinations shown in Tables VII and IX for the stochastic algorithm refers to the number of possible combinations left after removing rotamers that collide with the backbone. Therefore, the true number of possible combinations is much higher.

Section 3: Prediction of Loop Structures in Proteins The present invention is also particularly useful in solving the problem of correctly predicting loop structures in proteins. This particular implementation of the invention solves the difficult problem by allowing the determination of such predictions with considerable accuracy, without undue assumptions and combinatorial explosions.

The particular implementations of the invention described in “ Methods ” of this section were also tested against other methods known in the art as described in “ Results ”. It should be noted that these methods and results are presented for illustrative purposes only and are not intended to be limiting in any way. The interpretation of these results is then described in " Discussion ".

Construction of method loops can be accomplished by several strategies. Most of these use standard bonds and bond angles, but only change the dihedral angles.
This particular implementation of the method of the invention follows this general policy, but deviates from it in some steps.

Geometric Assumptions Figure 15 shows an example of 6 residues (0-5). Residue 0 and residue 5 are in the constant part of the protein. A search is performed for the conformation of residues 1-4. Loops were constructed simultaneously from both N- and C-termini (Mault and James, Proteins 1986; 1: 146-163) and loop closure is tested between residues 2 and 3. Such a construction strategy is such that when accumulating errors (ie, constructing a loop by a dihedral angle from one end to the other, an error in the first residue causes an increased amount of deviation in additional residues). Cause).

FIG. 16 shows dihedral angle definitions for certain residues. That is, the ψ of the residue n in the construction strategy is the ψ of the previous residue on the N-terminal side. The idea behind such a definition is that both φ _n and ψ _n determine the positions of the N and C atoms at residue n. When constructing a loop from the N-terminus (starting at residue 1 in Figure 15), the position of the nitrogen in residue 1, the first expected, should be specified according to the angle ψ of the previous residue. . The exemplary method of the invention described below assumes a Cα-C-N-Cα trans (180 °) structure. Therefore,
At residue 1, the position of Cα is specified according to this premise. The position of the carbonyl carbon of residue 1 is specified according to φ ₁ retrieved from the search (see below). The nitrogen position of residue 2 is specified according to ψ ₂ (which is regularly defined as ψ ₁ ) (and so on). When constructing the loop from the C-terminus, the position of the carbonyl carbon of residue 4 is specified by φ ₅ . Residue Cα is 180 with residue 5 Cα
It is in the ° position. The N position of residue 4 is specified according to ψ ₅ . Similarly, residue 3
Position of is specified based on the phi ₄ and [psi ₄ as defined in Figure 16. Therefore, the value of phi ₃ and [psi ₃ is not required.

Assignment of (φ; ψ) Possible Angles of Residues It is extremely unlikely to find all peptide sequences with a length of more than 6 residues in the structural database (Oliva et al., J. Mol. Biol). .199
7; 266: 814-830). Therefore, the method of the present invention is based on SWISS-PRO.
A search for a segment consisting of three overlapping residues of each loop in T is used (Bairoch and Apweiler, Nucleic Acids Re.
s. 2000; 28: 45-48). An array. ． ． ACGDEIL. ． ． Given a protein having: where “A” is residue 0 from FIG.
CGDE is a loop), the method of the present invention uses ACG, CGD, GDE, DEI.
, And search for the EIL segment. Broo for all (φ; ψ) angles of related residues in the segment detected by SWISS-PROT search.
khaven Protein Data Bank (Bernstein et al.,
J. Mol. Biol. 1997; 112: 535-542). The search is performed only on the 2nd and 3rd residues in each triplet so that any φ; ψ combination found should be associated with the order of the loop structure. Such a search results in multiple allowed conformations (
(Including rare conformations) can be obtained, and hundreds of pairs of φ; ψ angles can be obtained for a particular residue. These are stored as a database for later processing. If both φ and ψ angles differ from another pair of the same residue by less than 2 °, they are discarded from the database. To avoid any bias, the dihedral angle pair values of any protein searched for were excluded from the database for testing this particular protein.

[0169]Search of conformational space by stochastic algorithm   Large combination problems for medium and large loops
That is clear. For example, 1.15Å resolution bacteriorhodopsin complex
(Lücke et al., J. Mol. Biol. 1999; 291: 899-911.
) (Brookhaven file 1c3w.pdb) in the second loop
The number of combinations is 5.5 * 10²⁸Is. These database conformations
Only part of the loop may close loops that follow geometric criteria. Big problem
The method of the invention is used to reduce the texture. d₀Has an unknown angle pair of
Given a loop that does
involved in the possible structure of the function (see below). X_j= (X _j1 , X_j2． ． ． x_jd0) Was randomly selected ((φ_j1; Ψ_j1),
(Φ_j2; Ψ_j2) ,. ． ． , (Φ_jd0; Ψ_jd0)) A loop containing an angle
Information. Each conformation X_jAbout the cost
Function C_j= C (X_j) Can be calculated. The goal is to minimize all C
To find a conformation. This method is detailed in the two sections above
As mentioned above. Briefly, follow the steps below.

1. Randomly choose a value for each angle pair. That is, the sum constitutes the conformation of the entire loop. 2. We use a "cost function" to calculate the value of this conformation.

3. Continue to calculate the values of n such conformations (each conformation has all its variable values chosen randomly). 4. Create a histogram of the distribution of values for all sampled energies (n = 1000).

5. Compare all variable components that contribute to the part α (α = number of conformations) of the complete histogram in the upper region. 6. All components that contribute to the maximum conformation are removed, but none that contribute to the minimum conformation are removed .

7. The process is iteratively repeated until the remaining combinations can be exhaustively evaluated. At the end of this phase, there are many conformations that follow the criteria of geometric loop closure, but are not necessarily collision free. So in the next step,
Add side chains and evaluate loops by energy criteria. In various tests of this algorithm, at the end of this step, 10 ² -10 ⁵ conformations were left for further processing.

Score function in the stochastic phase The purpose of the stochastic phase is to create a population of loops that can potentially be closed. By using Equation 2, the loop that remains open is eliminated. The method of the present invention uses the cost function in Equation 2 to search the conformational space.

[0175]

[Equation 8] Where the distance d _i is shown in FIG. These are calculated after placing the last linking residue from the N- and C-termini. The value of r _i ^experimental is, N-C bond length (Shenkin et al, Biopolymers 1987;
26: 2053-85).

If there are less than M combinations of the loop scoring M (M is approximately 10 ² to 10 ⁵ ), then an exhaustive search is run to obtain the N minimum energy conformations of the loop. A backbone-dependent rotamer library (Dunblock and Karplus, J. Mol. Biol. 1993) was used to score the energy of the remaining loops.
230: 543-574; Dunblock and Karplus, Nat.
Struct. Biol. 1994; 1: 334-340), with the addition of the side chains of the loop. For atoms, a fused atom model was used (Weiner et al., J Amer. Chem. Soc. 1984; 10
6: 765-784). The backbone NH and the polar hydrogens of the side chains are clearly indicated. AMBER (Weiner and Kollman, J. Comp
． Chem. 1981; 2: 287-303; Weiner et al., J Amer.
Chem. Soc. 1984; 106: 765-784) binding and non-bonding energy terms were used with a distance dependent permittivity of ε = 2r (Equation 3).
Non-bonding energies are calculated for the interaction of the backbone and other residues with rotamers. 12-1 of hydrogen bonds between all polar hydrogens and possible acceptors
A zero potential was used. Satisfactory, but in order not to lose the solution, which may show some VdW collisions, the Lennerd-Jones repulsive energy was truncated at a value of 30 Kcal / mole for a particular atom pair.

[0177]

[Equation 9] The function was used in mean field form. Bower et al. (J. Mol. Biol. 199).
7; 267: 1268-1282), and any rotamer is given a probability in the backbone-dependent rotamer library, as implemented in the SCWRL algorithm. The energy of interaction from Equation 3 is multiplied by the probability (p) assigned from the relevant rotamer (the sum of rotamer probabilities is 1 for each residue). All rotamer-rotamer interactions, rotamer-backbone interactions, and backbone-backbone interactions were considered. A subset of residues with at least one atom at a distance of 10 Å from the native loop was included as a “template”. If the atom in Equation 3 is from the backbone, then its probability is p = 1. The binding energy terms included elongation (Equation 4), bending (Equation 5), and twist energy (Equation 6). The elongation energy (Equation 3) is calculated between the carbonyl carbon and the nitrogen closing the loop (d ₁ between residue 2 and residue 3 in FIG. 15).

[0178]

[Equation 10] _{"K b"} parameter for controlling the stiffness of the coupling spring (bond spring), _{r 0} defines the equilibrium length. A value of k _b = 100 was assigned to soften the energy function. Bending energy is calculated as (
Equation 5).

[0179]

[Equation 11] The “k _θ ” parameter controls the stiffness of the angle spring, while θ ₀ defines its equilibrium angle. Parameters specific to angle bending are
It is assigned to each of the linked triplets of atoms based on the type of atom. Two types of triplets were used. The first triplet was Cα-NC (d ₂ in Figure 15), where Cα is part of the previous residue. The second triplet is Cα-C-N
Where Cα and C are some of the previous residues) (d ₃ in FIG. 15).

[0180]   The twist energy is modeled by a periodic function (Equation 6).

[0181]

[Equation 12] The "A" parameter controls the amplitude of the curve, the n parameter controls the periodicity of the curve, and φ shifts the entire curve along the axis of rotation (τ). Torsional rotation-specific parameters are assigned to each of the bonded quartets based on the atom type. The angular twist energy (ie, the energy “price” of the deviation from the flat (180 °) amide bond) between the Cα-N—C—Cα atoms (d ₄ in FIG. 15) is then calculated. .

All results were evaluated by comparison with loop coordinates from high resolution X-ray crystallography by applying the coordinate root mean square deviation algorithm (cRMS). Such comparisons can be made by other methods (van Vlijmen and Karpl.
us, J. Mol. Biol. 1997; 267: 975-1001 and D.
eane and Blundell (Proteins 2000; 40: 135.
-144)) was only done on the backbones N, Cα, and C to allow comparison.

Results and Discussion The above tests can be used to see if the new stochastic search method is also applicable to loop construction, and to reconstruct structurally known loops of various sizes. The purpose was to see if. The example used was a transmembrane protein. The only extensive experimental example is bacteriorhodopsin containing seven transmembrane helices, which was recently studied by high-resolution crystallography (Leucke et al., J. Mo.
l. Biol. 1999; 291: 899-911). A search was applied to this structure (X-ray results at 1.55 Å resolution, PDB file 1c3w). The six loops of bacteriorhodopsin are shown in Table X. Loop 3 (CD, intracellular) and Loop 4 (DE, extracellular) contain 2 and 1 residues, respectively, and are not interesting test cases. In loop 5 (EF, intracellular) the coordinates are not included in the entry and therefore the quality of the results cannot be clearly evaluated. The rest of the loop:
1 (AB, intracellular), 2 (BC, extracellular), and 6 (FG, extracellular) are attractive test cases with 4-16 residues. To avoid bias, 1c3w. Did not include the pdb entry. The RMS value was 0.28-2.46 (Table XI) and the average value was 1.35. It is encouraging that the algorithm can, on the one hand, produce very small RMS in small loops and quite large RMS values in the case of very large loops.

Subsequently, Deane and Blundell (Proteins 2000;
40: 135-144) and van Vlijmen and Karplus.
(J. Mol. Biol. 1997; 267: 975-1001), an attempt was made to compare the effectiveness of stochastic loop predictions for globular proteins. Deane and Blundell used the ab initio loop construction method. Their algorithm selects polypeptide fragments from computer generated databases. Each fragment is defined by a representative set of 8 pairs (φ; ψ). This set of fragments is scored and classified using the RMS fit with the anchor region and the knowledge-based energy function. van Vlijmen and Karplus are 2
A search against a database of 130 loops from one protein was used. The best loop among many candidates was CHARMM applied to the backbone and C (β) atoms (Gunsteren and Karplus, JC).
omp. Chem. , 1980; 1: 266-274) non-bonding energy function (
(Without electrostatics). The method of the invention was tested on the longest 7 of their 11 loops. Table XII compares our results with those reported by the above two methods. Vlijmen and Karpl
Average RMS value of 2.3 (0.3 to 5.2) for us, Deane and Blu
Average RMS value of 2.1 (1.3 to 3.2) in the case of
MS value was 1.86 and range was 1.06 to 2.99. These lower average RMS values clearly demonstrate the superior quality of the method of the invention.

Furthermore, the algorithm yields a large collection of low energy loop conformations, which can be used to evaluate loop properties (eg flexibility) as well as from crystallography in the case of reconstructing known loops from PDBs. Can be further used to compare with the loop temperature factor.

During this study, some basic questions were raised. The first issue involved the accuracy of the approximation using standard bond lengths and angles. For this purpose, the method of the invention was used for the first loop of bacteriorhodopsin (see below). The RMS value between the predicted backbone and the experimentally obtained backbone was 0.280. The dihedral angles obtained from the true experiment were added to the angle database and the remaining dihedral angles were deleted. Therefore, the only option of the method of the present invention was to build according to the dihedral angle obtained from the experiment. It is expected that an RMS value of 0 will be obtained if the remaining angles and bond lengths are about the same as those obtained from the experiment. However, an RMS value of 0.204 occurred. This shows that such an approximation has a slight but nonnegligible effect.
This must be taken into account, especially when creating large loops where the accumulation of errors can distort the results.

The second problem is that φ differs by less than 2 ° (for both angles) from another pair of the same residue.
It was related to the accuracy of the approximation excluding the ψ angle pair. The above test was repeated with minor changes (ie increasing all dihedral angles obtained from the experiment by 2 °). Surprisingly, an RMS value of 0.198 resulted. Repeating the same test with all dihedral angles reduced by 2 ° gave an RMS value of 0.220. Such a slight difference indicates that the approximation is adequate.

The global optimization of the loop structure is
This is a difficult task due to the strong dependencies between variables. That is, one φ or ψ
Changes in angle can cause significant conformational changes throughout the loop. Therefore, the question raised concerns whether the method of the present invention, which successfully locates the protons and side chains, can generate a population of loops according to the geometric criteria defined in equation 2. did. Again, the first loop of bacteriorhodopsin was used. This is not extremely large for exhaustive searches, while still presenting combinatorial challenges. The 10,000 "minimum cost function" conformations of Equation 2 are shown in FIG. Both searches began with 54,330,000 combinations. The same global minimum was achieved by two different search approaches. The first 66 conformations had the same cost value (RMS). The worst error (RMS, 2018th solution) was 3.36%, and the difference in cost value was 0.006721Å. This test demonstrates that the method of the invention can effectively search for small and large loops and obtain significant results for these loops. These results strongly support the certainty of the method of the invention in solving this type of biomolecular problem.

The quality of the results should be investigated with the goal of achieving a negligible RMS compared to that obtained from the experiment. Here, the basic assumption is that the smaller the energy, the better the RMS. Like any other tool, R
MS has its own limitations. The user has to consider which atom should be superposed. The RMS values were compared in Table XIII, where different atoms were selected for superposition. Calculation of the loop RMS between N, Cα and C gave an average RMS value of 1.86. Addition of carbonyl oxygen increased the average RMS value to 2.10. Addition of protein residues attached to the loop increased the average RMS value to 2.62. The inclusion of these two residues is
One might consider reducing the RMS value because these coordinates are "correct". However, the opposite phenomenon is observed at least in our test case. In other words, when superimposing the atoms of the loop, RMS ignores the rest of the protein, ignoring geometric factors such as bond length and dihedral angle between loop and protein. Small RMS values do not necessarily indicate that the inner loop shape is acceptable, that the predicted loop has the proper shape (the loop is not left open) and does not collide with "known" portions of the protein. You should make sure. This phenomenon can be explained by the RMS superposition mechanism. The RMS function translates and rotates the predicted loop to overlap the known loop, ignoring the rest of the protein. Given the protein structure made up of wires, the "wire loop" can bend (by prediction) without changing the inner loop coordinates. Therefore, RMS = 0.0 can be achieved for residues m to n. On the other hand, this curvature can cause a large deviation from the protein structure, so that attaching a “correctly predicted loop” to the protein results in deviations of other residues from that protein position. , RMS increases significantly.

Section 4: Examples of Other Biological Problems The previous section provides detailed test results that illustrate the effectiveness of the present invention in solving a number of difficult biological problems. This section describes many other such problems and how such problems can be solved by the present invention.

Homology modeling. Homology modeling construction of unknown protein structures based on proteins found from X-ray or NMR studies requires "insertions" and "deletions" and mutations of peptide fragments compared to known structures. (
The homologous part of the target (to be constructed) is superposed one residue at a time over the homologous part of the known protein. Other portions may differ in length and are regularly encountered in loops of known proteins, beta hairpins, and random coil portions. Difference in length (
Each such manipulation requires a re-evaluation of the backbone coordinates in the non-homologous part, due to "insertions" and "deletions") and at least the position of the side chains near the modified part of the structure. Building the model using a fixed backbone that does nothing initially helps any scheme of mutations in the known protein structure. Substantial advances in solving this serious problem have already been achieved by the method of the invention.

In the present invention, by applying the above approach to loop construction, with or without simultaneous side-chain position determination, the "insertion" and "deletion" The impact can be dealt with.

Because of this problem, some additional challenges, such as the use of various force fields to assess energy, and alternative rotamer libraries with or without statistical weighting Can also be incorporated into the present invention.

Ring closure of peptides, peptidomimetics, and other cyclic structures. Many studies have shown the potential therapeutic importance of cyclic ("conformation-limited") peptides, rather than biologically labile linear peptides (Hruby, L.
if Sci. 1982; 31, 189; Alstein et al., J. Am. Biol.
Chem. 1999; 274: 17573). Theoretical study of the conformation of such peptides (Keasar and Rosenfeld, Foldin
g and Design 1998; 3: 379; Tieleman et al., Bi.
ophys. J. 1999; 76: 1757) is largely dependent on the "correct" ring closure option due to the high barrier between peptide conformations.
Many such cyclic peptides have been studied in solution by NMR (Ba
ysal and Meirovitch, Biopolymers 1999; 5.
0: 329).

Cyclization of active peptides and other linear molecules enhances binding to biological receptors by the expected reduction in entropy loss, enhances stability to digestion, and specificity and This is one of the selection methods such as for strengthening the selectivity. Alternative preliminary modeling for ring closure can greatly aid the design of such an annular structure. It is a function of many variables such as ring size, bond length, bond angle, and other factors.

This problem is very similar to the problem of loop structure prediction for the present invention. Cyclic peptides are generally smaller than loops, and a very small amount of "freedom" can be introduced into the backbone and side chain conformational flexibility. Also, a relatively small increase in Phi and Psi (backbone) angles is required for an exhaustive search for ring closure options.

Flexible docking of drug candidates to the active site. Computational methods for predicting binding of a ligand to its target are generally based on the search for the most stable binding conformation of the complex (Stryndaka et al., Nature Struct.
Biol. 1996; 3: 233). In theory, this corresponds to the crystallographically observed binding conformation (Rosenfeld et al., Annu. Rev.
． Biophys. Biomol. Struct. 1995; 24: 677; C
lark and Westhead, Comput. Aided Mol. Des
． 1996; 10: 337; Abagyan and Totrov, J .; Mol.
Biol. 1994; 235: 983; Trosset and Scheraga.
, Proc. Nat. Acad. Sci. USA 1998; 95: 8011).
. However, low energy conformations other than the global energy minimum of the complex may contribute to binding affinity. The new docking algorithm incorporates this assumption (Head et al., J. Phys. Chem. 1997.
101: 1609). However, DOCK, AUTODOCK, FLEX
Most of the flexible docking software, such as X, GOLD, and others, does not consider the flexibility of the target protein or biomolecule, but the active site protonation (pKa) state or potential changes in water content. I haven't.

Flexible docking is essential for testing the ability of most molecules to bind to biomolecular targets at different positions and variable binding methods. In flexible drug-protein interactions, structural changes (primarily conformations) can occur at both drug-protein interaction sites (active sites of enzymes, binding sites of receptor proteins).

In the context of the present invention, this is an extension of the problem of side chain localization and protein loop structure determination, but it requires the movement of the drug 6 degrees of freedom (translation + rotation) with respect to the active site of the biomolecule. There is a difference. The present invention must handle both side chain localization and loop (backbone change) predictions described above, but there is a further need to optimize the relative position of the biomolecular target and the ligand. ,
Optimization is applied simultaneously to both biomolecular targets and ligands.

These additional degrees of freedom can be arbitrarily introduced as variables, but there are special requirements. Moreover, optionally and more preferably, the problem is analyzed according to the method of the invention by adding an additional variable that is the relative distance of entities (eg, active sites of drugs and biomolecules).

Variables therefore include distance and angle variables for a total of 6 additional variables for translation and rotation. The present invention must handle both side chain localization and loop (backbone change) predictions described above, but there is a further need to optimize the relative position of the biomolecular target and the ligand. , Optimization is applied simultaneously for both biomolecular targets and ligands. Thus, the variables preferably include a total of six such variables: distance and angle variables for three types of translation along the X, Y, Z coordinate axes and three types of rotation about the same angle. included.

Structural comparison of flexible molecules. Traditional RMS approach and other superposition methods (Lemmen et al., Pac. Symp. Biocomput. 1999; 4;
82) is not suitable for comparing a very large range of conformations of flexible molecules.

Such a comparison allows the assessment of the probability that different molecules can bind to the same biomolecular site / target. Two different molecules may show similar binding affinities for the enzyme active site or receptor. In order to discover candidates for the "bioactive conformation" of both such molecules, the methods of the invention allow structural discrimination between such molecules to be optimized.

This problem is another conformational search of the present invention, where the function or quantitative parameter to be minimized is the RMS difference between the spatial positions of selected atoms in the two molecules. .

Construction of molecules from fragments. This is a classical problem in structure-based drug design (Krygeer et al., Structure 1999; 7: 2.
97), with considerable attention (Mizutani et al., J. Mol. Biol. 19).
94; 243: 310; Tomioka and Itai, J .; Comput. A
ided Mol. Des. 1994; 8: 347; Leach and Lewi.
S.J. Comp. Chem. 1994; 15: 233). Although several excellent approaches have been developed to study the affinity of molecular fragments for biomolecules (eg GRID (Wade et al., J. Med. Chem. 199).
3; 36: 140; Wade and Goodford, J .; Med. Chem.
1993; 36: 148; Boobbyer et al. Med. Chem. 198
9; 32: 1083)), combining fragments into molecules that can be drug candidates requires a huge amount of computational work. Again, the population of ligands should be followed rather than just one.

Although this indicates that some positions of molecular fragments are known from other studies,
Combining these fragments into specific and selective ligands is a major task in drug design, a complex task. A considerable number of programs (eg GRID) (Wade and Goodford, J. Med. Chem. 19).
93; 36: 148-156) may indicate the best interaction site between a molecular fragment (eg hydroxyl, amine, carbonyl, etc.) and the active site of a known protein structure. However, the number of potential combinations of such fragments is very large and requires an optimization process to construct drug candidates or lead compounds. In this case, the synthesizable structure and molecular weight,
The process must also be guided by chemical knowledge in order to achieve structures that are limited by lipophilicity and the like.

This problem differs from all other problems because chemical knowledge of synthesis and molecular stability must be introduced into the structural evaluation. In this case, the variables are fragments whose position in the "active site" has been previously optimized, as well as the assembly of the fragments into a viable complete structure, the interaction energy with the "active site" and the internal and hydration energies. The spatial position and orientation of the molecular “linking” fragments (eg, aliphatic, cycloaliphatic, and aromatic) used to assess

Folding of small proteins. Recent Short Review on "Energy Landscape in Non-biological and Biological Molecules" (Fra
unfelder and Leeson, Nature Structure. Biol
． 1998; 5: 757), the authors conclude that "the folding protein is evolutionarily selected so that its energy landscape resembles a funnel." These funnels are often (Wales et al., N.
Nature 1998; 394: 758), with a steep shape with a small barrier inside, but other funnel shapes exist. Chaperonin (Horovitz
Curr. Opin. Struct. Biol. 1998; 8: 93), there may be many accessible conformations near these "global minimums" that may be important for kinetic and functional studies. Search for these funnels (Dill and Chan
, Nature Struct. Biol. 1997; 4:10) became one of the important issues in modern biology, the "protein folding" problem. mainly,
On-lattice simulation (Sh
aknovitch, Curr. Opin. Struct. Biol. 1997
7:29) after years of studying small models, more modern simulations have attempted to fold peptide fragments or entire proteins (Dobson and Karplus, Curr. Opin.
Struct. Biol. 1999; 9:92). A very long simulation (1 μs) was recently developed (Duan and Kollman, Science).
1998; 282: 740), allowing folding of a 36 residue protein. Monte Carlo method (Hansmann and Okamoto, C
urr. Opin. Struct. Biol. 1999; 9: 177) and other stochastic dynamics (Sanderowitz and Still, J. Compu.
t. Chem. 1998; 19: 1294) is still valid. These energy-based methods do not readily find the native folding of proteins with more than 35-40 residues.

For the past 20 years, protein folding has been a key issue in biophysics. The method of the invention can be applied to a relatively small set of proteins with 50-80 residues depending on the primary structure of the protein. This approach is
In addition to the "global" minimum, many other low energy conformations that are close to the energy of the global minimum and that contribute to all the characteristics of proteins can be generated.

In a small protein of about 50 residues, the variables were the Phi and Psi angles along the backbone (6 or 12 rotations, each angle differing by 60 ° or 30 ° respectively) and It is a side chain rotamer. For the backbone only, with 6 rotations for each of the Phi and Psi angles, the problem size is 6 ¹⁰⁰ or about 10 ⁶⁶ . If there are additional rotamers to be simultaneously configured, the magnitude of the problem increases to about 10 ¹⁰⁰ . Therefore, the resulting calculations, which can be complex, can be done using the method of the invention.

It is understood that the above description is only intended to serve as an example, and that many other embodiments are possible within the spirit and scope of the invention.

[Table 1]

[Table 2]

[Table 3]

[Table 4]

[Table 5]

[Table 6]

[Table 7]

[Table 8]

[Table 9]

[Table 10]

[Table 11]

[Table 12]

[Table 13]

[Brief description of drawings]

[Figure 1]   3 is a flowchart of an exemplary method according to the present invention.

[Fig. 2]   1 is a schematic block diagram of an exemplary system according to the present invention.

[Figure 3]   Flowchart for hydrogen allocation algorithm.

FIG. 4 shows a molecule containing two carbonyls, one sp ² amide and one hydroxyl, which together form one population. The two carbonyls (1,2) act as acceptors, the hydroxyl donates one non-trivial hydrogen (3) and two non-trivial lone pairs of electrons (4,5), and the amide one trivial hydrogen. Donate (6). The atom 3 and the lone electron pair 4, 5 are one segment because they are bound to the same oxygen.

5A illustrates an exemplary initial two-dimensional matrix for the system of FIG. Hydroxyl hydrogen (3) can form hydrogen bonds with either carbonyl (1,2),
The hydroxyl lone pair (4,5) can form hydrogen bonds with the trivial hydrogen (6).

FIG. 5B is a diagram of a two-dimensional matrix after refinement. Since the two lone electron pairs of hydroxyl degenerate, one of them can be omitted (5 → 6). The omitted lone pair is automatically added after hydrogen and the first lone pair are located.

FIG. 5C is a diagram of a three-dimensional matrix created using the two-dimensional matrix to maintain all possible combinations. Each combination is evaluated and the best combination results.

FIG. 6 is a diagram showing an example of a “large” system. An initial two-dimensional matrix in large biological systems (eg proteins). Since it is considered that an attempt to create a three-dimensional matrix exceeds the capacity of a computer, the two-dimensional matrix is refined by removing high energy components.

FIG. 7 shows a “test” protein with 1186 amino acids. Of the 1186 amino acids, 13 are serines (labeled as a CPK model) (13 segments) and 1173 are glycines (0 segments). The probabilistic search starts with a total of 5.02 × 10 ¹⁰ combinations and 2 after 204 iterations.
． It became 7 × 10 ³ . After that, these were comprehensively evaluated. A global minimum for the position of hydrogen has been found.

FIG. 8 is a graph showing the natural logarithm of (total number of possible combinations) against the number of iterations in a pure “stochastic approach”. Five proteins are shown.

FIG. 9 is a graph of energy distribution in the 1st to 4th iterations for 5PTI (A), 5RSA (B), 2MB5 (C), INTP (D).

FIG. 10: Ribbon diagram of trypsin (INTP) and its polar residues. Many polar hydrogens form hydrogen bonds with water molecules, but the coordinates of any water molecule are not included in the PDB file.

FIG. 11 In finding the lowest energy conformation of 10,000,
Crambin (46) as a test example comparing a complete exhaustive search with a stochastic search.
The figure which showed the model of (amino acid residue). The cranbin backbone is shown with a ribbon. Atoms other than hydrogen are shown by the ball-and-stick model.

FIG. 12 is a diagram showing a comparison between a stochastic search and an exhaustive search when finding the lowest energy conformation for 1 to 10,000 conformers. The deviation (%) between the two searches is the bottom curve.

FIG. 13A: Detectable E. The figure which showed the percentage of the angle in E. coli ribonuclease HI. Of the 115 dihedral angles, 7 angles are not found in the rotamer library.

FIG. 13B. E. coli detected by the stochastic algorithm. coli ribonuclease HI
The figure which shows the percentage of the angle in.

FIG. 14: Single residue α values for 2 to 29 possible rotamers that are likely to be removed. An α (Δ) value is assigned for each number of rotamers. The greater the number of rotamers, the smaller α. The percentage of certainty has been calculated for each given rotamer number and α (□).

FIG. 15 shows an example of a 6-residue (0-5) loop. Residues 0-5 are part of the transmembrane helix. A search is performed for the conformation of residues 1-4. The method of the present invention is used to search the loop conformation space to find all possible loop closure conformations defined in Equation 2.

FIG. 16 is a diagram showing the definition of dihedral. The ψ of residue n in the construction strategy is the ψ of the previous residue on the N-terminal side.

FIG. 17 shows 10,000 “lowest cost function” conformations in a 4-residue study. Probabilistic and exhaustive searches have reached the same global minimum.
The first 66 conformations are in agreement.

─────────────────────────────────────────────────── ─── Continued front page    (81) Designated countries EP (AT, BE, CH, CY, DE, DK, ES, FI, FR, GB, GR, IE, I T, LU, MC, NL, PT, SE, TR), OA (BF , BJ, CF, CG, CI, CM, GA, GN, GW, ML, MR, NE, SN, TD, TG), AP (GH, G M, KE, LS, MW, MZ, SD, SL, SZ, TZ , UG, ZW), EA (AM, AZ, BY, KG, KZ, MD, RU, TJ, TM), AE, AG, AL, AM, AT, AU, AZ, BA, BB, BG, BR, BY, B Z, CA, CH, CN, CR, CU, CZ, DE, DK , DM, DZ, EE, ES, FI, GB, GD, GE, GH, GM, HR, HU, ID, IL, IN, IS, J P, KE, KG, KP, KR, KZ, LC, LK, LR , LS, LT, LU, LV, MA, MD, MG, MK, MN, MW, MX, MZ, NO, NZ, PL, PT, R O, RU, SD, SE, SG, SI, SK, SL, TJ , TM, TR, TT, TZ, UA, UG, US, UZ, VN, YU, ZA, ZW (72) Inventor Gold Bram, Amiram             Israel 93501 Jerusalem Sim             Sean 20 (72) Inventor Glick, Meir             Israel 59503 Bat Yam She             Shet hayamin 2 F-term (reference) 5B056 BB64 BB65                 5B075 ND20 UU18

Claims (26)

[Claims] 1. A method for searching every corner of a combinatorial space, the space featuring a plurality of combinations, each combination consisting of a plurality of elements, wherein the steps of the method are performed by a data processor. (A) providing measurable quantitative parameters for each combination to determine whether the result of searching the combinatorial space in every corner is successful; and (b) the combinatorial space. Selecting a plurality of combinations in the above to form a selected combination; (c) calculating a value for the quantitative parameter for each of the plurality of selected combinations; ) Determining the effect of each factor on the value of the quantitative parameter, and (e) reducing the effect according to the effect. Both holds one combination, the method comprising the steps of providing a result of searching the combined space throughout. 2. The step performed before step (a), further comprising the step of determining a structure for a plurality of combinations of combination spaces such that there is interaction between the elements. The method of claim 1, wherein 3. The element according to claim 2, wherein each element is a variable having a value, and the quantitative parameter is calculated according to the value of the variable for each combination and the interaction between the variables. The method described. 4. The method of claim 3, wherein the quantitative parameter is a cost function. 5. The method of claim 4, wherein each variable has one discrete value for any particular combination. 6. Step (e) comprises: (i) rejecting a value if the value does not consistently improve the cost function; and (ii) each combination characterized by the value for the variable. The method according to claim 1 or 5, further comprising: rejecting. 7. Step (e) comprises: (iii) determining whether the number of remaining combinations is below a minimum number, and (iv) if the number of remaining combinations is above the minimum number, Step (
The method according to claim 6, further comprising repeating steps c) to (e) at least once. 8. The step (e) further comprises: (v) evaluating each remaining contribution according to a parameter if the number of the remaining combinations is less than the minimum number.
The method described in. 9. The method of claim 8, wherein step (v) is performed by an exhaustive search of the remaining combinations. 10. Step (i) comprises: (1) creating a plurality of combinations by randomly giving values to the variables; and (2) calculating the values for the cost function. (3) determining that the effect of the value is a negative effect if the value is found in a plurality of combinations having a value for the cost function that is below a predetermined minimum value, The method of claim 6, wherein the method is performed. 11. The cost function, wherein the value for the variable is found in a plurality of combinations having a value for the cost function below a predetermined minimum value, and the value for the variable is above the predetermined desired value. 11. The step (3) is performed if found in a combination having the value for and the value for the cost function is determined to be below the predetermined minimum value. the method of. 12. The method according to claim 6, wherein each value is relative to a position of a polar proton in a biomolecule, and the cost function is a calculated minimum energy value for the polar proton. . 13. The combination comprises: (i) parameterizing the atoms of the biomolecule; (ii) a hydrogen atom and a lone pair of electrons; a trivial hydrogen atom, a polar hydrogen atom, and a non-trivial lone pair of electrons; 13. The method of claim 12, determined according to the steps of: categorizing; and (iii) adding trivial hydrogen atoms to each combination. 14. The cost function is a non-bonded energy function paired as follows:
number [Equation 1] (In the formula, A_{i, j}Is the repulsion parameter for two atoms (i, j), B _{i, j} Is the polarizability attraction parameter, and q₁Is the partial charge and r_{i, j}Is an atom
Distance, and ε is the permittivity) Method according to claim 13, characterized in that it is calculated according to: 15. Each combination of the combination space includes a position with respect to a side chain of an amino acid in a protein, and the cost function is a calculated minimum energy value with respect to the side chain of the amino acid. The method according to 6. 16. Each combination is formed from rotamers for each side chain and the step of forming a combination comprises eliminating rotamers that collide with the backbone of the protein. 15. The method according to 15. 17. Each element is a rotamer, and the effect of the element on each combination is to sample a plurality of combinations, and for each population, measure the quantitative parameter across the sampled plurality of combinations. 16. Method according to claim 15, characterized in that it is determined by examining the distribution. 18. The combination according to claim 6, wherein each combination in the combination space includes a structure for a loop of a protein, and the cost function is a minimum energy calculation value for the structure of the loop. Method. 19. The method of claim 18, wherein the structure for the loop is determined according to multiple angle pairs between residues of the protein. 20. A value for each angle pair is randomly selected to create each combination, each combination having an associated value for the cost function, and an associated value for which the value exceeds a predetermined threshold. 20. The method of claim 19, wherein the value is removed if it contributes only to the combinations that it has. 21. The step (e) evaluates each combination according to the presence or absence of collision between side chains of amino acids in the loop, and rejects the combination if the collision exists. 21. The method of claim 20, further comprising: 22. Each combination of the combinatorial space includes structures for all loops in a target protein, the cost function being a minimized energy calculation for the structures of the loops, and The structure is determined according to multiple angle pairs between residues of the protein, step (a) comprising providing a predetermined structure for a known protein having homology to the target protein, and step (e) ) Between the side chains of amino acids within the loop within the target protein and between the remaining portion of the target protein and the loop after the loop has been evaluated by comparison with the structure of the known protein. Evaluate each combination according to the presence or absence of collision, and discard the combination if the collision exists. 7. The method of claim 6, further comprising the step of allowing said structure for said target protein to be determined according to said structure for said known protein. 23. A value for each angle pair is randomly selected to form each combination, each combination having an associated value for the cost function, and the value being a predetermined threshold value for the value. 23. A method according to claim 22, characterized in that it is eliminated if it only contributes to combinations having a higher associated value. 24. Each combination of combinatorial spaces includes positions for multiple moieties within a cyclized molecule, the structure of the cyclized molecule is determined from the structure of a linear molecule, and the cost function is the linear function. 7. The method according to claim 6, wherein a calculated minimum energy is obtained for the position of the portion by comparison with the structure of the molecule. 25. Each combination of combinatorial spaces comprises an assembly of molecular fragments to form one molecule, said assembly for connecting each molecular fragment in place to at least one other molecular fragment. 7. The method of claim 6, wherein the cost function is a calculated minimization energy for the predetermined position of the molecular fragment in the structure of the assembly. 26. Each combination of the combination space comprises at least one part of the structure of the first and second entities, each part comprising a variable for rotation about an angle between each said part and said first part. Is defined by a relative distance between an entity and the second entity, the relative distance being defined by a variable for transformation along a coordinate axis, the cost function being the first entity and the second entity. 7. The method of claim 6, wherein the minimized energy calculation is for interactions between, for the distance, and for at least a portion of the first entity and the second entity.