Method and Means for Synthesis-Based Simulation of Chemicals Having Biological Functions

1. Field of the Invention

The invention concerns the simulation, by a programmed, general purpose digital computer, of chemicals that have biological functions, and more specifically, concerns computer-implemented simulation of such chemicals based upon their synthesis and their assayed biological activity.

2. Description of the Related Art

Currently, there is no method or means for simulating the relationship between a chemical produced in a laboratory, a biological function of the chemical as experimentally measured in an assay, and the outcome of the experiment where the chemical is assayed. The ability to simulate such an experiment could vastly improve prediction of the outcome, enabling the use of the information developed in simulation to better identify and develop chemicals that will produce desired effects on selected biological receptors.

Currently, chemicals having biologically efficacious functions are synthesized from chemical libraries that contain huge numbers of chemicals. Based upon clinical structures, experience, intuition, or strictly random sampling, sets of chemicals are selected from a chemical library, and they are combined in chemical reactions subject to known reaction conditions to produce product chemicals. The biological function of the product chemicals is measured experimentally in an assay procedure, with the degree to which a product chemical achieves a desired biological function being indicated by a numerical value. After synthesis, thousands, or tens of thousands, of product chemicals may be assayed, with the assay results being determined and ranked by known means.

These large numbers of product chemicals are synthesized in order to discover one, or a few product chemicals that react with a protein to produce a desired effect, such as one that would inhibit or impair the progress of a disease. Assume that the disease is inhibited by attachment of an enzyme to a known location on the protein. The known location is called a "receptor site". Assume that attachment of a product chemical to the protein at receptor sites of interest inhibits enzyme attachment, thereby strongly indicating that the product chemical would also inhibit the progress of the disease. In this regard, attachment of the product chemical to a receptor site of interest would mimic the therapeutic effect of a (possibly scarce) naturally occurring compound. In this case, the product chemical would comprise a synthetic substitute for the scarce compound.

In any event, an assay may be conducted to determine whether, and to what extent, the product chemical binds to receptor sites of interest by coating a well or a container with the protein and then washing the protein with a solution containing the product chemical. Excess solution is emptied from the coated well and the compound for which a block, or substitute, is sought is "tagged" with another compound and, so tagged, is brought, in solution, against the protein. The extent to which the product chemical binds is evaluated by introducing another chemical that "develops" the tagging compound. In this regard, development activates the tagging compound in the sense that the tagging compound shows a certain color. Introduction of the protein to a spectrophotometer enables the detection of chemical activity by the measurement of color. Thus, if the tagged compound is not blocked, it will bind to the protein at the receptor sites and its presence will be indicated by the color of the tagging compound. If altogether blocked, the color of the tagging compound will not be detected at all or, if partially blocked, a smaller amount of the color will be detected.

Color-based assay techniques are well known, and apparatus are available to conduct photoassay of a great number of separate reactions. One such device is a spectrophotometer available from Molecular Devices that automatically conducts photoassay of a plurality of reactions in an array of reaction wells. Such a spectrophotometer can be coupled to a programmed, general purpose digital computer to process the results of an M.times.N array of reaction wells. For example, a MacIntosh brand computer available from the Apple Corporation, running the SOFTMX application program, also available from Molecular Devices, provides an output that presents photoassay results of a plurality of normalized assay result values, referenced to maximum and minimum values of a tagging color.

While the synthesis/assay process is well known and widely used, it is time-consuming and very manpower-intensive. To the extent that the process is iterative, because it is manually performed, desirable results emerge slowly, and patterns of desirable activity by product chemicals may be difficult to detect.

Manifestly, automated means and methods for simulating the direct synthesis/assay procedure have the potential of vastly increasing the speed of the process and quickly identifying patterns or trends of desirable results. However, to date only the evaluation phase of the assay step has been automated, and that only to a rudimentary degree. As yet, no effective mode of simulating the synthesis/assay process has emerged. Attempts have been made to create simulation models based upon molecular structure of product chemicals. However, such attempts have yet to find an effective model that links how a chemical is made and what its structure is, once synthesized. Accordingly, a manifest need exists for a method and means that simulate product chemicals in a way that directly reveals how they may be made and what their biological functions are.

This invention provides a method and a means that simulate product chemicals in terms of how they are made and the biological effects that they produce, thereby enhancing the efficiency with which drugs may be discovered and developed. Before the conception of this invention there was no way to mathematically represent a relationship between a chemical produced in a laboratory, a biological function of the chemical as experimentally measured in an assay, and the outcome of the experiment where the chemical is assayed. The invention is based on the inventors' mathematical model of the outcome of an assay experiment performed on a chemical synthesized according to a prescription of starting chemicals and reaction conditions.

An essential object of this invention is therefore to provide an alternative or an improvement to the prior art structure-based simulation model.

The invention is based upon the critical realization that synthesis and assay of a chemical may be statistically described as a physical system. Importantly, the synthesis of a chemical is modeled as a preparation procedure that is evaluated by a registration method. The preparation procedure is specifically modeled as a procedure by which a chemical is produced, while the registration procedure is modeled as an assay. Critically, there exists a function .mu. that underpins the description of a physical system modeled as a combination of the preparation and registration procedures.

The invention is embodied as a machine or as a machine-executed method for simulating chemical function using a database in which data representations of chemicals, data representations of synthesis procedures, and a set of functions representing assay procedures are stored. In this regard, a first set of first data structures representing chemicals of the chemical library is obtained from the database. A second set of second data structures representing industrially performable reactions by which chemicals represented by the first data structures are combined is also obtained from the database. The first and second sets are combined to produce a third set of data structures that represent experimentally performed, or hypothetically specified synthesis procedures in which chemicals in the chemical library and reactions in the reaction library are involved. The feasibility of a synthesis procedure represented by a third data structure is simulated by subjecting the third data structure to a function of the set of functions, such that the function produces a simulation value. A product chemical is synthesized according to the specified synthesis procedure and assayed to produce an assay result. The assay result is compared with the simulation value with reference to a tolerance limit, and another function may be selected from the set of functions based upon the comparison.

In this manner, a model simulating the synthesis and assay of large numbers of product chemicals can be evaluated and adaptively reconfigured by testing against actual synthesis/assay results produced from a much smaller set of product chemicals.

FIG. 1 illustrates an iterative process for composing a function .mu. that predicts the feasibility of a specified synthesis procedure in addition to the biological function of the produced chemical;

FIG. 2 is a block diagram of a database used in the invention;

FIG. 3 illustrates a tree that represents a reaction sequence by which a product chemical is synthesized;

FIGS. 4A and 4B illustrate how a synthesis tree representing synthesis of a product chemical may be iteratively described; and

FIG. 5 illustrates a combination including a programmed, general purpose digital computer, the database of FIG. 2, and a synthesis/assay procedure by which the invention is practiced.

With reference to FIG. 1, assume a collection of available chemicals, denoted by C, and a collection of reactions which are industrially performable, denoted by R. A chemical compound, w, which is synthesized from a set of chemicals C'C through a list of reactions R'R is denoted by w(C', R').

A traditional synthesis/assay process that leads to drug discovery is described by the pseudocode in Table 1. The synthesis and assay in this process are actually-performed procedures. With large C and R, this enumeration-styled process is very time- and fund-consuming.

With our invention, regularities between , C and R are found. We represent the relationship between them by the mathematical formula of equation (1).

In equation (1) ' represents the result of a simulated assay of a simulated product chemical w. Based on quantum theory, the mathematical function .mu. (hereinafter, the "prediction function") is likely to be decomposed into a set of sub-functions, .mu..sub.i, i,=1, . . . , n, each of which describes one aspect of this problem (a quantum dimension): ##EQU1##

This mathematical function .mu. is composed based on actual experimental data E obtained by synthesis and assay, where, E={e.sub.i =(w.sub.i, .sub.i).vertline..sub.i =Assay(w.sub.i), i=I, 2, . . . }. The composition of the function .mu. is depicted by FIG. 1 and described by the pseudo code in Table 2. If the computation results obtained by .mu. match the experimental results within a range of tolerance that lies between 0 and a tolerance limit, it is determined that a prediction function has been identified that may reliably predict drug discovery. In this regard, .mu. is a "prediction" function in that it predicts the feasibility of a specified synthesis procedure in addition to the biological function of the chemical(s) produced according to the procedure.

Now, assume composition of a reliable prediction function .mu.. Before real synthesis and assay experiments are performed respectively in a chemical laboratory and assay laboratory, the invention composes a drug which is most likely to be synthesizable and have a particular biological effect of interest. This process, picking out the most promising drug composition solution from the vast fabricated chemical compound library, is an optimization process, and is illustrated by the pseudocode given in Table 3.

Actually, during the utilization of the prediction function, the prediction function will also be optimized due to access to more experimental data becoming available (E expends during this process).

2. Composition of the database

In order to manage the data and information used in this invention in a way sufficient for prediction function development and for drug discovery optimization, we have designed a database 10 that is illustrated in FIG. 2.

The database 10 comprises four libraries: a chemical library 12, a reaction library 14, a synthesis library 16, and an analysis function library 18.

All of these libraries are updatable and extendable either by known maintenance means, by a user, or both.

3. The data structures of the entities in the database

The data structures of entities of the database directly determine the time and space (memory) efficiency of the database.

(1) Chemical library (C) 12:

The chemical library 12 is a first set of first data structures, each representing an entity in the library, wherein each entity in the chemical library is a chemical, represented by a list of properties, such as molecular geometry, bond length and bond angles, molecule weight, and so on.

The chemicals can be classified into many categories in terms of the type of tags in molecular structure.

.A-inverted. c.sub.i .epsilon. C, can be represented as a first data structure, which may be treated as an instance of an object, as follows:

(2) Reaction library (R) 14:

The reaction library 14 comprises a set of second data structures, each representing an entity of the library wherein each entity of R may represent necessary conditions for an industrially performable reaction.

.A-inverted. r.sub.i .epsilon. R, is represented as a second data structure, which may be treated as an instance of an object described as follows:

(3) Synthesis library (w) 16:

The synthesis library 16 comprises a set of third data structures, each representing an entity of the library wherein each entity w .epsilon. W represents a reaction sequence which can be represented as a tree-shaped graph as shown in FIG. 3.

This synthesis tree can be expressed in an iterative way by means of an object defined as follows:

if tree.sub.-- state=0, this tree is actually a chemical reagent (c) instead of a synthesis tree. Otherwise, it is a sub-tree. For example, FIGS. 4A and 4B illustrate a synthesis tree's iterative description.

(4) Function library 18:

As to a hypothetically specified synthesis w, the invention predicts its feasibility and other properties through mathematical analysis based on the quantified parameters of the synthesis and involved chemicals. In other words, a mathematical function is composed, as follows:

where ' is a random variable representing some kind of probability of something.

The prediction function .mu. can be decomposed into a collection of relatively simple sub-functions (.mu..sub.1, .mu..sub.2, . . . , .mu..sub.p):

Each sub-function, say .mu..sub.i, takes a synthesis tree (w) as its argument. All sub-functions .mu..sub.i are pre-defined and assigned a unique name that enables use of the sub-function in the invention. The data structure for this library is defined as follows:

A dedicated language may be provided for composition of a formula expression using the functions in the function library 18. For example, assume composition of a formula as follows:

.mu..sub.i .epsilon. function library

.alpha..sub.k is a coefficient

Automated database processing may automatically apply this formula to a candidate syntheses tree in the synthesis library 16, and compare the computational results with experimental results. The formula can be dynamically modified to approach the experimental data.

FIG. 5 is an illustration of an industrial application of the invention in which a programmed, general purpose digital computer 50 having standard hardware and software components, including a memory 52, is conventionally coupled to a storage facility that may include, for example, a direct access storage system including components 54a, 54b, 54c and 54d. The components 54a-54d of the storage facility are used to store the storage library illustrated in FIG. 2. The storage facility is shown as comprising individual components 54a-54d only for purposes of illustrating the library architecture of the database shown in FIG. 2. Relatedly, the storage facility component 54a includes the chemical library 12, the storage facility component 54b has stored in it the reaction library 14, the storage facility component 54c provides storage capacity for the synthesis library 16, while the storage facility component 54d stores the functions that make up the function library 18.

The programmed, general purpose digital computer ("computer") 50 is programmed in a conventional way to execute one or more computer programs that embody the pseudocode of Tables 2 and 3. In this regard, Tables 2 and 3 are embodied as a composition procedure 58 and a simulation engine 60 that may comprise, for example, one or more computer software programs that include program code in compiled, assembled and linked form that cause the computer 50 to perform the functions embodied in the pseudocode of Tables 2 and 3.

Embodiment of the invention in the industrial application of FIG. 5 includes programming the computer 50 with an assay application 62 that receives and processes experimentally-derived assay results from an actually-performed synthesis/assay procedure 64, 74.

The synthesis/assay procedure 64, 74 includes the steps of selecting chemicals from a chemical library 66, and subjecting the chemicals to reaction conditions 68 in a synthesis laboratory to produce a plurality of product chemicals 70. In FIG. 5, the plurality of product chemicals are represented as an M.times.N array. In an assay laboratory, the assay step 74 includes subjecting the product chemicals 70 to assay conditions 78, with assay results being measured by a spectrophotometer 80 from which assay results are provided at 82 as an input to the assay application 62.

In operation, the function .mu. is composed according to procedures set forth in Table 2. Table 2 takes experimental results derived from synthesis/assay of one or more product chemicals. Each of these results is represented by a respective assay result .sub.i. Each assay result is a measure of biological activity of a particular product chemical represented by a synthesis tree w.sub.i. In this regard, during the iterative composition of the function .mu., each w.sub.i comprises a third data structure 94 that is composed in the memory 52 by combining first data structures 90 obtained from the chemical library 12 with second data structures 92 obtained from the reaction library 14, where the first data structures 92 represent chemicals of the chemical library 66, while the second data structures represent reaction sequences 68 by which chemicals represented by the first data structures are combined to produce product chemicals whose assays are represented by the input 82 to the assay application 62. Thus, during composition of the function .mu., there is, (for each specified synthesis procedure w.sub.i) a particular assay result .sub.i comprising actual experimental results produced by the synthesis/assay procedure 64, 74. Thus, in line (4) of Table 2, a calculated assay result '.sub.i is obtained and stored at 96 in the memory 52 for the third data structure instance of w.sub.i that corresponds to a synthesis procedure that produces a product chemical w for which the assay result .sub.i has been obtained. In line (5) of Table 2, the absolute value of the difference between the simulated assay value '.sub.i and the product of e.sub.i and .sub.i is compared to a tolerance limit. If the absolute value of the difference is not greater than the tolerance limit, the function .mu. is deemed reliable. It may equivalently be stated that if the absolute value of the difference is less than or equal to the tolerance limit, the function .mu. is deemed reliable. Otherwise, the function is modified by selection of another function from the function library 18 and another iteration of the composition procedure in Table 2 is executed.

When a reliable function .mu. is found, this is deemed the "simulation function". The simulation function is used by the simulation engine 60 to execute the procedure of Table 3. In Table 3, a subset .GAMMA. of the solution space C.times.R is obtained, with members of the solution space combined as described above to produce the subset as a set of third data structures, each in the form, preferably, of the synthesis tree described above, which are subjected to the simulation function. In this regard, each synthesis tree represents a product chemical, while the simulation function represents an assay performed on the set of product chemicals represented by the third data structures. For each synthesis tree, Table 3 obtains a solution to equation (1) and compares it to tolerance value .zeta.. For every synthesis tree that produces an equation (1) result that exceeds the tolerance value, the simulated chemical is synthesized, assayed, and if the assay result exceeds the tolerance value, the product chemical is deemed successful.

During the iterative composition of the function .mu. (Table 2) and during the discovery process of Table 3, the synthesis library 16 is built, with each third data structure stored in the library representing a product chemical. Each third data structure may be stored with a value representing the actual experimentally-determined assay result for the represented product chemical. Eventually, the synthesis library 16 will be deep and broad enough to be mined for function trends leading to the rapid discovery of product chemicals that exhibit desirable functions.

1. Quantum Theory and Structures of Molecules

Quantum theory originated from the study of atoms or small molecules, and other small physical systems. The formulation of the theory took several decades of active research to get complete.

The structures of molecules can be computed using the equations of Quantum Theory. The mathematical equation representing the theory (most notably the Schrodinger equation), although difficult to study, has been approached by various methods (e.g., electron orbital theory, Hartree Fock theory), giving meaningful approximations to the molecular structures. Actually, exact solutions are found only in very simple cases, in which only a handful of electrons are involved.

The study of structures is still a very active field of research. It has proved to be very useful in predicting the chemical reactions among the micro-molecules.

2. Structure-Based Prediction

There are many works in the direction of a structural approach in predicting the ligand binding of G-protein coupled receptors. Here the ligand is usually a micro-molecule, but the protein-receptor consists of a formation of macromolecules.

There have been attempts to set up a Hamiltonian for the macrosystem. The Hamiltonian for the macrosystem could be useful in predicting qualitative properties or the values of statistical physics quantities.

Assume that the Hamiltonian formula is correctly worked out. Trying to exactly solve a system containing thousands of electrons under one Hamiltonian not only is beyond the present (or near future) mathematical computational ability, but also could be meaningless. One reason is that is could take too long for the system to go from an arbitrary state to the ground state (or some more or less stable state), and the process could be disturbed by many small random factors. So even if the ground state is the most probable state, the probability density in its neighborhood could still be too small to count.

Another way to attack the problem is to study chemical reactions between a ligand and local region of the macromolecule. That is, set up a local Hamiltonian which involves a ligand, and the local part of the macromolecule. One may develop a parallel theory in the local region. However, we should point out that although this is similar to the Hamiltonian of microsystems, the problem could be substantially more difficult as there is no clear cut limit for the local region.

There have been none to rare reports of positive results in these directions.

3. Alternates to the Structural Approach

It is an open question if the quantum theory of atoms and molecules, reasonably correct for small systems, also works for macrosystems.

Since the "potential ligands" are small molecules, their structures can be worked out reasonably well by theoretical computations (using the method of X-ray crystallography). Many invariants could be defined: from the geometric structures of the configuration of atoms, the wave (and spin) function of electrons, to the measures of electromagnetic vibrations of atomic nuclei, and functions derived from these.

Assuming that the structural approach works, some of these structural parameters may be correlated with the experimental outcome of the prediction function .mu.. If consistent correlation can be detected, then regression models could be used in predicting the outcome of .mu. (with probability) for other potential product chemicals.

Many parameters may be used for the prediction function. Following are examples of a few, out of hundreds of structural or non-structural parameters.

(i) Let I denote the measure of inertness of the potential product chemical. Consider, as an example, the following simple linear regression model:

where c.sub.1 and c.sub.2 are constants, .epsilon. is a random variable (the error term). The constants c.sub.1 and c.sub.2 can be determined by some routine calculation, based on the method of best square mean approximation.

(ii) Consider the invariants of the 3D geometric structure, such as diameter, volume, density, shape of its boundary, and functions of these. Examples: ##EQU2## These parameters may be used in the following model:

or

(iii) From the electron orbital theory, each electron is associated with four quantum numbers: n, l, m.sub.s, and s.

For each set (n, l, m.sub.s, s), define the following functions. Let N(n, l, m.sub.s, s) be the number of electrons whose quantum numbers are given by the quadruple, and for each of these N electrons, let E denote the total energy of the electron (including the potential energy with respect to the atoms), and let P be the distance from the boundary of the effective region of the electron to the boundary of the molecule.

Note that E and P could be regarded as a collection of variables indexed by the (n, l, m.sub.s, s) and k, where k=1, 2, . . . , (n, l, m.sub.s, s). Let J be the collection of the indices (n, l, m.sub.s, s, k), then, E.sub.i, P.sub.i will be used to denote all of these variables.

One may define an order of these indexes in J, according to, say, the energy E. In fact, we can assume J is fixed: just let E and P.sub.i be some adjustable small, negative numbers, if the electron does not exits. The assumption here is that electrons of relatively low energy level are less likely to be free.

This suggests another regression model: ##EQU3## where .alpha..sub.i, .beta..sub.i, c.sub.i and d.sub.i are all constants, to be determined by the least square method, and .epsilon. is a random variable.

4. Use of Synthesis Tree According to the Invention

The parameters introduced in section 3 may be defined for each component in the synthetic tree where each component is indexed by k .epsilon. K={1, 2, . . . , n}. For each k, we may use the method of 3 for a prediction model:

where .alpha..sub.k is a constant vector, and x.sub.k is a vector variable. We then combine these, and do a search: ##EQU4## Or, more generally, ##EQU5## where ›x.sub.k Xx.sub.i ! denotes the vector each of whose coordinates is the product of a coordinate of x.sub.k and a coordinate of x.sub.i.

The most general formulation is: ##EQU6##

We note that this last formula reduces to (4.1) once all but one variables are fixed. Thus, in a way, for each k .epsilon. K, the constants b.sub.k, .alpha..sub.k of (4.1) are related to the constants c.sub.k, c.sub.ik, . . . in (4.2)

5.0 Alternate Embodiment of Function Composition

A prediction function for .mu. may be constructed, using regression analysis, and based upon the experimental data, e, , e .epsilon. E. In this regard, a regression model, denoted by M, will have the following form: ##EQU7## where x.sub.i =x.sub.i (w)'s are the parameters used in the model, c.sub.i 's are the constants whose optimal values are determined by regression analysis (e.g. the method of least square means), and .epsilon. is the random error term. We assume a collection of such models.

The procedure in finding the prediction function .mu. then consists of two steps.

Step 1: For each regression model M, as in (5.0), minimize the following functional: ##EQU8## over

Where is the set of real numbers. Now, let (c.sub.1, c.sub.2, . . . , c.sub.n) denote the minimizer, and let: ##EQU9## Step 2: Compare Error(M) for different regression models M. Pick the one, say M, so the Error (M) is the smallest.

Then announce the choice of prediction function ##EQU10## where x.sub.1, x.sub.2, . . ., x.sub.N are the parameters used in M, and (c.sub.1, c.sub.2, . . . , c.sub.N) is the optimizer obtained in Step 1 for the model M.

In other words, the alternate emboidment of the procedure for function composition is:

Step 1: Optimize the coefficients (c.sub.1, c.sub.2, . . . , c.sub.n)

Step 2: Choose the best regression model.

Other embodiments and modifications of this invention may occur to those of ordinary skill in the art in view of these teachings. Therefore, this invention is to be limited only by the following claims, which include all such embodiments and modifications when viewed in conjunction with the above specification and accompanying drawings.