BACKGROUND OF THE INVENTION There is a need for predicting the evolution of existing and proposed chemical or biochemical mechanisms, in particular those mechanisms involved in living cell processes. There exists high profile cellular modeling efforts founded on classical bulk chemical analysis or switching network theory, such as: 1) the E-cell project led by Masaru Tomita of Keio University [1-3]; 2) The BioSpice project for cellular network simulation led by Adam Arkin at Uni- versity of California at Berkeley [4, 5] ; 3) the whole-cell metabolic simulation efforts of the Metabolic Engineering Work- ing Group led by Bernhard Palsson [6-8] and George Church [9, 10]; 4) The Systems Biology Project of Lee Hood [11]; and 5) the Dutch Silicon Cell Project [12]. Several of these efforts have attracted the attention of the popular press (The Scientist, May 10,1999 ; Time Magazine, Aug.

6, 2000; Scientific American, Aug. 2001).

From elementary physics, we know that many physical phenomena are most conveniently explained by using different sets of equations (or models) for differ- ent levels of resolution that explain the bulk properties of particulate matter.

Atomic-scale physics is b@9 @xplained by Schrodinger's equation, macro-scale physics is best explained by Newton's equations of motion and grand-scale phys- ics is best dealt with by Einstein's theory of General Relativity. In a completely parallel manner, existing methods undertake to simulate biochemical systems at only one of three different levels of"resolution"using three fundamentally differ- ent, but well tested computational approaches.

1) A continuum model (resolution 10-100 nm) that is most suited for modeling metabolic (i.e. small molecule) fluxes and flows ; 2) a particulate model (1-10 nm scale) which is most suited for modeling mac- romolecular interactions and cellular superstructure ; and 3) an atomic model (0. 1-1 nm scale) which is most appropriate for modeling single macromolecules or small numbers of macromolecules and catalysis over a limited time (10-12 s to 10-9 s) and volume (100 nm3).

It is desirable to provide a simulation architecture and set of methods that are fully and continuously salable both in spatial resolution and in time. There is a need to not only track the chemical state but also the position of each individual chemical particle. It is further desirable to provide a model which is implicitly 4- dimensional and jointly represents the spatial and temporal distributions of cell activity.

By and large, existing efforts have adopted two main approaches, namely, the continuum approach or the network approach. The former is based on the fun- damental assumption of a chemical continuum where the concentration of each chemical species varies continuously and"smoothly"across the reaction space.

This basic assumption makes it possible to represent biochemical systems as large sets of differential equations and draw on the rather extensive mathemati- cal theory which exists to solve such systems. However, this formalism is ill- suited to describe the low concentrations and discontinuous spatial distributions of molecules that can exist at the intra-cellular level. In such cases, simulations based on the continuum approach tend to become highly unstable and often fail to converge towards a solution. As a result, current approaches are mostly lim- ited to one-dimensional representations where the cell model is formulated as a set of ordinary differential equations capturing only the temporal evolution of the system and sacrificing all spatial resolution.

An alternative formalism has been proposed and is being used by Arkin et al.

[4, 5]. This network approach describes biochemical reaction systems by casting such systems in the framework of conventional switching circuits. Once again, this is a one-dimensional representation where several key simplifying assump- tions are made and where the evolution in the spatial distribution of molecules is not described or captured by the model. In this approach, biochemical systems are modeled by using biologic components in a manner that is that is analogous to the modeling of electrical circuits by electrical logic components. Developed in the field of electrical engineering, this switching circuit formalism draws on a rich foundation of analysis of modeling and analysis techniques. Although modifica- tions to this approach have made it possible to avoid some of the problems as- sociated with continuum modeling, the fact remains that the network approach largely lacks in chemical formalism and the scalability of this approach is at the very least unclear.

It is, therefore, desirable to provide a framework that captures the discrete and probabilistic nature of the chemical aspects of molecular behavior found in a cell while also formally representing the spatial distribution of atoms, molecules and molecular assemblies and other components of the cellular machinery.

SUMMARY OF THE INVENTION In recent years, the life sciences have witnessed an explosion of technologies enabling the high-throughput generation of data in the areas of genomics, pro- teomics and metabolomics. Concomitant with the development of these tech- nologies has been the emergence of data management and data interpretation technologies collectively referred to as bioinformatics.

The data that is generated from genomics, proteomics and metabolomics is largely a qualitative inventory of cellular constituents. Consequently, bioinformatics has primarily evolved as textual based method for statistically based comparisons that draw static relationships between elements within and <BR> <BR> <BR> between different inventories. Taken together, these approaches have shed considerable light on 1) the components list of the cellular inventory, 2) what the components look like, 3) how the components fit together and, 4) how the components function. What these approaches fail to address is a dynamic description of a living system. What is ultimately desired and required to move biology to the next logical plane are high throughput methods for quantifying cellular activity, whose data can then be used to for simulating living systems in silico.

Ideally, these"virtual"organisms will be manipulated and visualized at different levels of molecular resolution and will be able to respond, adapt and evolve to fit their virtual environments.

The present invention provides a system and method for simulating the behavior of processes associated with cellular life. Any specified real or proposed set of chemical or biochemical reaction mechanisms can be simulated as they evolve both in space and in time. Said simulations include sets of reaction mechanisms describing metabolism, signaling pathways, biomolecular interactions, cell growth and division, post transcriptional modifications, protein translation, gene tran- scription, DNA replication and other systems of coupled reaction mechanisms associated with cell life.

The present invention is quite flexible and can be endowed with the ability to op- timize existing cellular processes and evolve new ones through virtual natural se- lection. Unlike similar attempts at creating a simulation environment for cell bio- chemistry, the present approach has several distinct and significant advantages : 1. The simulation output is highly graphical in nature and amenable to various visualization technologies. Individual reaction elements, be it atoms, molecules or groups of molecules, are displayed as they traverse the physical reaction space and collide with each other, prompting reaction events according to a set of conditional probabilities ; 2. The graphical user interface allows for rapid prototyping of complex biological pathways by the user, without requiring extensiv@ programming skills ; 3. Both the architecture of the data structure and that of the computational meth- ods embodied in the present invention are fully and continuously salable. This scalability applies to both the size scale and the time scale.

4. By its very nature the models formulated using this system and its methods is implicitly 4-dimensional, providing a representation of the biological system which is distributed both in space and evolving in time.

5. The computational engine uses a novel combination of logical and probabilistic methods and operations to simulate the motion of individual chemical particles as well as the changes in chemical state of these particles as they participate in chemical reaction events.

6. By using this salable probabilistic approach to simulating movement and re- action events, the present simulation environment offers a seamless transition between the stochastic world of the very small and the apparent deterministic na- ture of the macroscopic world.

7. By using populations composed of arbitrary numbers of individual chemical reaction agents, this simulation approach is highly amenable to use with mas- sively parallel computational platforms and environments.

The present invention provides a system and method for predicting the evolution of existing and proposed chemical or biochemical mechanisms, in particular, those mechanisms involved in living cell processes.

What distinguishes the present invention from all of the above simulation efforts is its inherent 4-dimensional nature as well as its scalability in both time and space.

In this invention, a framework is established that captures the discrete and prob- abilistic nature of the chemical aspects of molecular behavior found in a cell while also formally representing the spatial distribution of atoms, molecules and mo- lecular assemblies and other components of the cellular machinery.

The discrete nature of this computational approach is inherently robust to discon- tinuities in the spatial and temporal distribution of chemical or biochemical reac- tions. As such, the simulation system embodied in this invention will continue to provide stable and accurate solutions at very low concentrations or when these atoms, molecules or molecular assemblies are dispersed according to irregular spatial patterns. By the same token, describing a chemical reaction as the result of interactions occurring between sets of discrete physical objects is a formula- tion that is highly amenable to massively parallel computation.

By portraying chemical and physical interactions between physical objects as stochastic events, this invention makes it possible to accurately describe the phenomenology of enzyme-based kinetics and thermodynamics, as well as other detailed phenomena without addressing these phenomena explicitly. Simply in- creasing the number of chemical reaction particles present makes it possible to <BR> <BR> <BR> <BR> seamlessly move from the stochastic behavior of very small-scale events to an accurate representation of the deterministic behavior of the macro-scale envi- ronment.

The present invention comprises a method and system for simulating a living cell process of an organism. The system creates a simulation space in the computer that includes a three-dimensional volume described in Cartesian coordinates.

The system also includes the dimension of time whereby the simulation of a living cell process is developed over a user-defined period of time. The simulation space is sized according to the living cell process to be simulated and to the particles of predetermined chemical species to be introduced into the simulation space during the course of the simulation.

Each living cell process is described as a chemical pathway which, in turn, is a composition of one or more reversible chemical reactions that have at least one reactant and at least one product, both of a known chemical species. Each reversible chemical reaction has a statistical probability of occurring in both the forward and reverse direction determined by stochastic methods. The parameters of every possible chemical reaction and of all chemical species involved in the simulation are stored in computer memory. In addition, particles of predetermined chemical species are introduced in the simulation space. The particles may represent an influx of atoms, molecules or molecular assemblies unrelated to the organism or may represent the organism itself including the cell membrane of the organism or particles situated within the cell.

The simulation is first initialized by defining the physical volume of the simulation space based on the process being simulated and the particles to be introduced.

The simulation then defines an initial position and velocity for each particle to occupy the simulation space. In operation, the simulation is organized into a series of iterations of predetermined periods of time as defined by the user. For each iteration, the movement of each particle in the simulation space is predicted based on stochastic methods and navigational cues based on and applied to the position and velocity of the particle at the start of the iteration period. In addition, each possible chemical reaction that can take place between neighbouring particles is determined. The reactions that are likely to occur are then simulated using a set of stochastic methods by removing those particles that would be consumed by the reaction and by inserting those particles that would be produced by the reaction. An inventory is then taken of each particle in the simulation space that includes information regarding the chemical species, population, position and velocity of each particle now existing in the simulation space. This information is then used in the subsequent iteration. This procedure is repeated for each iteration until the simulation is completed.

At predetermined times during the simulation, information or data regarding the particles in the simulation space is taken and placed in data files. The information or data is also used to generate a visual representation of the simulation occurring in the simulation space on a computer display monitor.

The simulation tool may be configured as an object-oriented computer program or as a serially written computer program. Preferably, the simulation tool is configured to operate on a multi-processor computing platform. The object- oriented approach permits particles, chemical pathways and chemical reactions to all be described as objects, each having attributes and properties specific to that object. This approach has the advantage of each object having all relevant information associated with it but has the disadvantage of increasing the computational time to perform a simulation. The serial coding approach, conversely, permits the attributes and properties of chemical species and chemical reactions to be organized into shared tables of information. This approach also streamlines the computations and reduces the time required to perform a simulation.

During the initialization routine, the number of particles situated in the simulation space and their position are determined prior to running the simulation. In on preferred embodiment, this information is used to then divide up the simulation space into a number subspaces, one or more subspaces being assigned to one or more computer processors. Furthermore, the boundaries of each subspace are determined in such a way that the computations conducted for each said subspace require essentially the same computational resources. In this fashion, the computational load of performing the simulation throughout the whole simulation space is distributed across the multiple processors as evenly as possible. After each iteration, this process is repeated to re-size the subspaces so that they each require substantially the same amount of computational effort.

This is done to optimize the computational load on each processor as particles will migrate throughout the simulation space during the simulation thereby changing the distribution and density of particles in each subspace. If this is not done, the computation speed of the simulation will slow down as some processors will become bogged down with a large number of computations for their respective subspace while other processors will be lightly loaded or idle should their respective subspace become lightly populated or depleted of particles. By redistributing the size of each subspace after each iteration, the computational load on each processor can be kept substantially the same thereby maximizing the computational resources of the computing platform running the simulation and minimizing the simulation run time.

Broadly stated, one aspect of the present invention is a method for simulating at least one living cell process of an organism in a simulation space on a computer having at least one processor, said process including at least one of the group consisting of metabolism, catabolism, enzyme function, biochemical kinetics, gene and protein expression, physiological regulation, molecular biology and cell morphology, the method comprising the steps of storing a plurality of environment configuration parameters in a computer-readable medium accessible by the computer, said plurality of said parameters including the physical parameters of said simulation space and the chemical and physical characteristics of at least one said living cell process of said organism; initializing the simulation space on the computer in accordance with said stored environment configuration parameters ; defining said at least one living cell process in the simulation space as at least one chemical pathway comprising of at least one predetermined reversible chemical reaction between a plurality of predetermined chemical species and storing the chemical characteristics of each of said at least one predetermined reversible chemical reaction and the attributes of each of said predetermined chemical species in a computer-readable medium accessible by the computer, said characteristics including the statistical probability of the occurrence of said at least one chemical reaction in both the forward and reverse directions; populating the simulation space with a plurality of predetermined particles selected from the group consisting of atoms, molecules and molecular assemblies, and storing the characteristics of each of said predetermined particles in a computer-readable medium accessible by the computer, said characteristics including the chemical species, population, position and velocity of each of said particles in the simulation space; simulating at least one chemical reaction between neighbouring particles in the simulation space over at least one iteration of a predetermined time period if the probability of said at least one chemical reaction occurring exceeds a predetermined threshold; determining the particles in the simulation space consumed by said at least one simulated chemical reaction following each said at least one iteration; determining the chemical species and population of each of said particles in the simulation space produced from said at least one simulated chemical reaction following each said at least one iteration; determining the position and velocity of each of said at least one particle in the simulation space following each said at least one iteration; and generating data representing at least one of the group comprising of the chemical species, population, position and velocity of each particle in the simulation space following each said at least one iteration.

Broadly stated, another aspect of the present invention is a method for simulating at least one living cell process of an organism in a simulation space on a computer having at least one processor, said process including at least one of the group consisting of metabolism, catabolism, enzyme function, biochemical kinetics, gene and protein expression, physiological regulation, molecular biology and cell morphology, the method comprising the steps of defining the physical parameters of the simulation space and storing said physical parameters in a computer-readable medium accessible by the computer; initializing the simulation space on the computer with said stored physical parameters; defining at least one living cell process as at least one chemical pathway object and storing the attributes of said at least one chemical pathway object in a computer-readable medium accessible by the computer ; associating at least one sub-reaction object with each of said at least one chemical pathway object, each of said at least one sub-reaction object representing a predetermined reversible chemical reaction and having a plurality of reaction properties, said properties including the statistical probability of the occurrence of said reversible chemical reaction in both the forward and reverse directions, and storing said reaction properties in a computer-readable medium accessible by the computer ; associating at least one substrate object with each of said at least sub-reaction object, each of said at least one substrate object comprising of a predetermined atom, a predetermined molecule or a predetermined molecular assembly, and storing the attributes of each of said at least one substrate object in a computer-readable medium accessible by the computer; defining a plurality of particle objects in the simulation space, each of said particle objects being a member of a predetermined parent class having predetermined attributes whereby each of said particle objects inherits all of the attributes of its parent class, each of said particle objects associated with a predetermined chemical species having attributes whereby each of said particle objects inherits at least one attribute of its associated chemical species, each of said particle objects representing an atom, a molecule or a molecular assembly, and storing said attributes of each of said particle objects in a computer-readable medium accessible by the computer; computing the motion of at least one particle object in the simulation space over at least one iteration of a predetermined period of time and storing the position and velocity of said at least one particle object in a computer-readable medium accessible by the computer; determining a list of possible chemical reactions that can occur between said particle objects in the simulation space following each said at least one iteration and storing said list in a computer-readable medium accessible by the computer ; determining which of said possible chemical reactions are iikeiy to occur ; simulating at least one of said possible chemical reactions likely to occur if the probability of said possible chemical reaction occurring exceeds a predetermined threshold; determining the particle objects consumed and the particle objects produced by said at least one simulated chemical reaction and the storing the chemical species, population, position and velocity of at least one particle object in the simulation space following each said at least one iteration in a computer-readable medium accessible by the computer ; and generating data representing the chemical species, population, position and velocity of each particle object in the simulation space following each said at least one iteration.

Broadly stated, another aspect of the present invention is a method for simulating at least one living cell process of an organism in a simulation space on a computer having at least one processor, said process including at least one of the group consisting of metabolism, catabolism, enzyme function, biochemical kinetics, gene and protein expression, physiological regulation, molecular biology and cell morphology, the method comprising the steps of defining at least one living cell process in the simulation space as at least one chemical pathway comprising of at least one predetermined reversible chemical reaction having at least one reactant and at least one product, and storing the chemical characteristics of said at least one predetermined chemical reaction in a computer-readable medium accessible by the computer, said characteristics including the chemical species of the at least one reactant consumed in said at least one reversible chemical reaction, the chemical species of the at least one product produced in said at least one reversible chemical reaction and the probability of occurrence of said at least one reversible chemical reaction in both the forward and reverse directions; storing the attributes of a plurality of predetermined chemical species in a computer-readable medium accessible by the computer ; defining a plurality of particles capable of existing in the simulation space, said particles selected from the group consisting of atoms, molecules and molecular assemblies, each of said particles capable of being a reactant or a product of said at least one reversible chemical reaction ; populating the simulation space with a plurality of predetermined particles, and storing the characteristics of each of said predetermined particles in a computer-readable medium accessible by the computer, said characteristics including the chemical species, population and initial position and velocity of each of said predetermined particles in the simulation space ; building the simulation space in a computer- readable medium accessible by the computer, and initializing the simulation space in accordance with the chemical and physical characteristics of said at least one chemical pathway and the characteristics of said predetermined particles ; simulating at least one reversible chemical reaction of said at least one said chemical pathway of said at least one living cell process in the simulation space over at least one iteration of a predetermined time period if the probability of said at least one chemical reaction occurring exceeds a predetermined threshold, said simulated at least one chemical reaction comprising of at least one particle ; determining the probable position and velocity of at least one particle in the simulation space following each said at least one iteration, and storing the position and velocity of said at least one particle in a computer- readable medium accessible by the computer; determining the chemical species and population of at least one particle in the simulation space following each said at least one iteration, and storing the chemical species and population of said at least one particle in a computer-readable medium accessible by the computer; and generating data representing the chemical species, population, position and velocity of said at least one particle in the simulation space following each said at least one iteration.

Broadly stated, another aspect of the present invention is a system for simulating a living cell process of an organism in a simulation space on a computer, the system comprising a computer having at least one processor, said computer controlled by an operating system; a graphical user interface operatively connected to said computer via said operating system ; a computer-readable medium operatively connected to said computer via said operating system ; and a computer program product recorded on said computer-readable medium having computer program logic for simulating at least one living cell process of said organism, said logic capable of being run on said at least one processor via said operating system, said process including at least one of the group consisting of metabolism, catabolism, enzyme function, biochemical kinetics, gene and protein expression, physiological regulation, molecular and cell morphology.

BRIEF DESCRIPTION OF THE DRAWINGS FIGURE 1 is a block diagram of a first embodiment of a system and method for simulating one or several biochemical pathways in accordance with the present invention.

FIGURE 2 is a schematic representation of a graphical user interface used to specify the biochemical processes to be simulated, in accordance with another illustrative aspect of the invention.

FIGURE 3 is a block diagram of the object types and the hierarchy of these ob- ject types used in the first embodiment of the present invention.

FIGURE 4 (comprising component FIGURES 4A and 4B) is a block diagram of the simulation method used in the first embodiment of the present invention.

FIGURE 5 is a diagrammatic example of data entry, data processing and visuali- zation, illustrating possible particle motion and particle reaction outcomes.

FIGURES 6A and 6B show computer monitor screen displays of the pathway window, which allows the user to specify the biochemical reactions to be used in the simulation, grouped together as pathways.

FIGURE 7 shows a computer monitor screen display of the simulation window, which allows the user to observe the displacement of the individual atoms, mole- cules and molecular assemblies in the modeled cellular process.

FIGURES 8A and 8B show computer monitor screen displays of the graph win- dow, wherein user-requested graphs of population statistics are displayed graphically and updated dynamically with simulated results generated by simula- tion methods in accordance with one embodiment of the present invention.

FIGURES 9A and 9B show diagrammatic examples of the central data tables list- ing the relevant attributes of particle agents, the chemical reactions, and the chemical species in the serial application of the simulation methods used in the preferred embodiment of the present invention.

FIGURE 9C shows diagrammatic examples some of the different types of output files, including detailed snapshot files, chemical assay summary files, and simu- lation computer runtime performance files generated by the serial application of the simulation methods used in the preferred embodiment of the present inven- tion.

FIGURE 10 is a block diagram of the serial implementation of the simulation method described in the preferred embodiment of the present invention.

FIGURE 11 is a block diagram of the initialization phase of the serial implementa- tion of the simulation method described in the preferred embodiment of the pre- sent invention.

FIGURE 12 is a block diagram of the checkpoint strategy and re-initialization from existing data in the serial implementation of the simulation method de- scribed in the preferred embodiment of the present invention.

FIGURE 13 is a block diagram describing the management of the chemical state of each particle agent in the serial implementation of the simulation method de- scribed in the preferred embodiment of the present invention.

FIGURE 14 is a block diagram describing the execution of chemical reactions in the serial implementation of the simulation method described in the preferred embodiment of the present invention.

FIGURE 15 is a block diagram of the navigation strategy used in the serial im- plementation of the simulation method described in the preferred embodiment of the present invention.

FIGURE 16 is a block diagram describing the computation of the position and ve- locity of each particle agent in the serial implementation of the simulation method described in the preferred embodiment of the present invention.

FIGURE 17 is a block diagram describing the computation of the conditional probability distribution of the possible new locations for each particle agent in the serial implementation of the simulation method described in the preferred em- bodiment of the present invention.

FIGURE 18 is a graph of a number of objects, as a function of iteration or time step for a set of reversible enzyme-catalyzed reactions involving enzyme E, in- hibitor 1, substrate S, intermediate product ES and final product P, for two differ- ent combinations of initial substrate and initial inhibitor concentrations, [So] and [lo], as an example of simulation methodology according to the preferred em- bodiment of the invention.

FIGURE 19 is a graph of Initial Rate Vo, as a function of initial concentra- tion/number of objects [So], for the same reaction system as employed in Figure 18, showing that the data generated by the simulation adhere to behavior de- scribed by the Michaelis-Menton equations.

FIGURE 20A shows the competitive inhibition reaction system involving enzyme E, inhibitor 1, substrate S, intermediate product ES, enzyme inhibitor complex El, and final product P, as employed in the simulations for which data are shown in FIGURES 18 and 19.

FIGURE 20B shows graphs of the initial reaction velocity and the reciprocal of the initial velocity, as plotted against the initial substrate concentration (top graph) and its reciprocal (bottom graph), respectively, for the competitive inhibi- tion reaction of FIGURE 20A.

FIGURE 21A shows the reaction system involving enzyme E, inhibitor 1, sub- strate S, intermediate product ES, final product P, and inactive enzyme- substrate-inhibitor complex ESI, representing the basic set of reactions involved in producing uncompetitive inhibition, as used in simulation, in another aspect of the invention.

FIGURE 21 B shows graphs of the initial reaction velocity and the reciprocal of the initial velocity, as plotted against the initial substrate concentration (top graph) and its reciprocal (bottom graph), respectively, for the uncompetitive inhi- bition reaction of FIGURE 21A.

FIGURE 22A shows the reaction system involving enzyme E, inhibitor 1, sub- strate S, intermediate product ES, final product P, inactive enzyme-inhibitor com- plex EI and inactive enzyme-substrate-inhibitor complex ESI, representing the basic set of reactions involved in producing noncompetitive inhibition, as used in simulation, in another aspect of the invention.

FIGURE 22B shows graphs of the initial reaction velocity and the reciprocal of the initial velocity, as plotted against the initial substrate concentration (top graph) and its reciprocal (bottom graph), respectively, for the noncompetitive in- hibition reaction of FIGURE 22A.

FIGURE 23 schematically depicts the cellular fluxes and reaction system for a hypothetical example illustrating the methodology of the present invention, in one embodiment thereof.

DETAILED DESCRIPTION OF THE INVENTION The present invention is a system and method that demonstrates that <BR> <BR> <BR> populations of biochemical reaction elements acting as cooperative computing agents in a virtual environment can accurately represent the phenomenology of cellular processes at multiple levels of resolution. The object-orient and serial encoded embodiments are described below.

A) Object Oriented Implementation Simulation Structure-Object Oriented Implementation The basic constructs used to define the structure of the variable and process space for the living cell simulator are listed in Figure 1. As discussed earlier, the simulation space is populated by a number of discrete objects each representing atoms, molecules or molecular assemblies (101). These particle objects belong to classes representing their chemical species or substrate class (102) and are listed in a central object repository (100). Interactions between these classes of objects are managed through associative relationships encoded in reaction or sub-reaction objects (103). Once again building on a hierarchy of classes, reac- tion objects are grouped together into major pathways (104) based on their con- tribution to specific metabolic functions. Class objects are also listed in a central repository (105).

The state of these objects is altered according to a number of processes. Ob- jects which are to be used actively in the simulation must be configured by way of a number of object definition procedures (106). Further more, the runtime pa- rameters for the simulation, for example duration, output requirements, etc... must be set be applying methods incorporated in a number of initialization proce- dures (107). The actual main simulation run is managed by the simulator's main engine (108) which includes 2 major components : the object motion engine (109) and the object reaction engine (110). These engines process active parti- cle objects which are listed in the main object list (111). As the simulation is conducted output data is being written to an output data array (112), the contents of which are to be post-processed into various output formats by the graphical output (113) and text output (114) procedures. Execution of the simulation is managed as another process by the main operating system (115) which governs activity of the central processing unit or CPU (116). Interaction between the user and the computer system as well as with the simulation itself is provided by the user interface (117).

The avenues available to the user for interaction with the simulation are repre- sented graphically in Figure 2. The hierarchy of graphical interfaces and their relationship with objects in the simulation space are also presented in Figure 2.

The main graphical window (200) provides visual access to the three basic com- ponents of the user interface, namely the pathway window (201), the simulation window (202) and the graph window (203). The pathway window (203), also shown in Figure 6, allows the user to specify the biochemical reactions to be used in the simulation, grouped together as pathways (104). This interface ele- ment also allows users to adjust a number of run-time parameters required to manage the execution of the simulation. These are transmitted to the simulation space (204) where the data is manipulated by calling upon the reaction engine (110), the object motion engine (109) and the random number generator (206).

The methods embodied in these software modules process substrate (102) and reaction (207) data, continuously adjusting the state of the simulation space (204). Similarly, results which are generated by post-processing populate the summary information space (205) which serves as a data source to produce the user-requested graphs in the graph window (202), also shown in Figure 8. Fig- ure 8 (Figures 8A and 8B) shows population statistics that are displayed graphi- cally and updated dynamically with simulated results generated by simulation methods of the present invention. Detailed results present in the simulation space (204) are a) so presented graphically in the simulation window (201), where the user can observe the displacement of the individual atoms, molecules, and molecular assemblies. An example of this window is presented in Figure 7, which depicts a simulation output screen showing the instantaneous position and chemical state of each particle in the simulation space, updated dynamically with simulated results generated by the simulation methods of the present invention.

Simulation Procedure-Object Oriented Implementation The data manipulated and processed by the methods serves to define the ob- jects and object hierarchy described in Figure 3. At the highest level of aggrega- tion, entire biochemical or chemical pathways are grouped under a global Path- ways Object List (300). Each element of this list is a separate group of path- ways, the contents of which are found in the corresponding pathway object list (301). Each individual pathway object (302) is described by a number of attrib- utes including a number of distinct chemical or biochemical reactions or subreac- tions (303). The reactions describe interactions between a number of participat- ing chemical species or substrate objects (304). Chemical substrate objects are described in terms of a number of chemical and physical attributes. In addition, chemical substrate objects contain references to all atoms, molecules or molecu- lar assemblies which bear this chemical identity. These particle objects (305) are basic elements in the simulation and are described in terms their position in the simulation space as well as by a number of attributes inherited from their parent class the substrate object (304).

The execution of a simulation run is described in overview as a process flowchart in Figures 4A and 4B. A given simulation run must start with a general initializa- tion of the simulation space (400). This step includes assigning values to run- time parameters as defined by the user in the pathway window (203). If the cur- rent iteration number does not match or exceed the user-defined maximum num- ber of iterations, and the stop flag is not set to a logical true, the interrupt condi- tion (401) is not met and the simulation continues. Execution proceeds by com- puting the motion of each object in the simulation space using the object detec- tion engine (402), or ODE, which calls and executes the method ode. MoveObject. As dictated by conservation of momentum, this method will at- tempt to move a given object one spatial increment in the direction in which it is currently traveling. Should such a move result in a collision with a neighboring particle, this collision is resolved and a new direction computed by the method ode. Find-Resolution.

Once a new position has been obtained the reaction engine is called (Find New Reactions ()) to manage reactions which might occur between adja- cent particles (403). This engine first selects a type of reaction or reaction case (404) which might apply to this chemical type by drawing on the information stored in the subreaction table SubReactioninfo (SRI). Four specific reaction cases have been defined. Cases 0) and 1) involve the conversion from single reactant to one or more products and the reverse reaction where multiple prod- ucts are converted back to a single reactant respectively. Cases 2) and 3) in- volve the forward conversion of multiple reactants to multiple products and the reverse reaction respectively. Once a reaction case has been selected, a list of reactions of this type is extracted from the SRI table using the method reac- tant. getPossibleRxn () and a specific subreaction selected (404). For a given case and specific subreaction, the method then conducts a search of the local neighbourhood of the host particle object and identifies adjacent particles (405) using ode. GetNeighboursOf. The method CanReactionOccur is then called to evaluate the neighborhood and determine if the adjacent particles identified will support the selected subreaction (406). In conjunction with this, since chemical reactions are treated as stochastic events in this invention, a probability function is applied to compute the probability that this reaction will indeed occur (4-07) should the neighborhood prove suitable. Where the neighboring particles do support the selected subreaction (406) and the probability exceeds a given criti- cal value (409), the reaction is given the go ahead. The approved reaction is not processed immediately ; instead it is added to a queue by the method Queue. adareadant). Shoula coneitions (406) and (4107) not be met, the algo- rithm returns to the list of potential subreactions of the selected case type and the next candidate reaction is selected and evaluated as described above. All reac- tions accepted into the queue are processed in a subsequent step by the method Resolute Reactions (Queue).

For all queued reactions, the method will examine the spatial implications of cre- ating product objects. In cases where the space required to position product par- ticles is not locally available, the reaction is ignored and the method proceeds to evaluate the next instance in the queue. In cases where space is available to position new products, these products are generated and reactants removed, in accordance with the reaction equations originally entered by the user (410, 411), by the method AddNewReactants. Example applications of this method are dis- played in a diagram format in Figure 5. This diagram depicts various ways in which particle objects will fuse to form one product or split to generate multiple product objects.

Following the completion of a reaction, particle object flags are reset (412) and updated positions are provided to the graphical output engine (413). The simula- tion window and summary graph windows are also both updated (414,415). A new iteration number is assigned and a test conducted (401) to determine if the simulation is to continue. In cases where the stop flag assumes a value corre- sponding to a logical true, the simulation is halted and a post-processing step is called (416) to produce a file containing the numerical output. When this post- processing is complete, execution of the simulation is terminated (417).

B) Serial Code and Distributed Computing Implementation <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> Simulatio Structure-Serial Code and Distributed Computing Brmplement-<BR> <BR> <BR> <BR> <BR> <BR> <BR> glen The basic constructs used to define the structure of the variable and process space for the living cell simulator as implemented in the serial embodiment are listed in Figures 9A and 9B. As discussed earlier, the simulation space is popu- lated by a number of discrete particle or agents each representing atoms, mole- cules or molecular assemblies. In the serial code embodiment, particle agents are now represented in the code by entries in a central list instead of code ob- jects. Each row entry in the master table part data tbl (1001) contains a concise description of a single particle agent using a minimal number of attribute fields.

These attributes include but are not limited to the active status of the particle agent, the position in 3D space, the velocity components in 3D space, and the chemical identity or species of the particle agent. Data describing the attributes of each different type of agent is contained in a single central table speciesattr (1002). This table includes but is not limited to data on the relative populations of each chemical species, the dimensions of the particles of said species type, the fraction of agents of a given chemical type initially located in various regions of the simulation space, as well as the molecular weight and other attributes of par- ticle agents of this chemical species type. Each row describes in this central species attribute table lists the numerical attributes associated with a specific chemical species. The names of each of the species are listed in a separate character array or table identified as species tags (1003) in Figure 9A.

As with the object-oriented implementation, existing particle agents can be con- sumed and new particle agents generated through the execution of chemical re- actions. Interactions between the various chemical species or types of particle agents are specified by the stochiometric coefficients stored in the table reacttbi (1004) shown in Figure 9B. The assignment of chemical species state to new particle agents generated by a chemical reaction as well as the number of said particle agents generated are dictated by the statement of a specific chemical reaction. Each possible chemical reaction is a row entry in table reacttb) (1004).

The probability that a specific reaction will occur given that the appropriate local conditions are met is determined by a set of two threshold probability values listed the table p_react (1005). If a local region of space around a particle agent contains particle agents of the necessary chemical species types, in the propor- tions specified in reacttbi, then the reaction will proceed from the left to the right with the forward probability Pf specified in the first column of the corresponding row of the table p_react. Similarly, the same reaction will proceed from the right to the left hand side of the reaction equation with the reverse probability Pr speci- fied in the second column of the corresponding row of the table preact.

The general outline of the steps involved in the definition and execution of a simulation problem by the serial embodiment of the method is presented dia- grammatically in Figure 10. The user defines a simulation problem by assigning values to parameters in a number of input files. These can be divided into files describing the characteristics of the simulation problem and environment (1101), the external stimuli applied to the simulated system (1102), and files describing the computer hardware and software environment (1103). Characteristics of the simulation problem include but are not limited to the initial population of particle agents, their location and chemical species type as well as the chemical reac- tions dictating the interaction between these chemical species types and the probabilities of said interactions. Examples of external stimuli include but are not limited to the description pulse trains specifying the rate of addition of new parti- cle agents, their location and chemical species types. This feature allows the user to simulate among other things the application of various drug therapies to the model system. Finally the hardware and software environment is described by parameters which specify a serial or a parallel mode of execution, the number of computer processors available, new problem initialization or re-start from a previous execution, as well as the location of required input files and the destina- tion for report files containing simulation results. Several types of results are generated by each simulation run. Typical output data files are shown diagram- matically in Figure 9C. At user defined intervals, the simulation will write a table of results to a snapshot file (1006). These files list, among other things, the loca- tion and type of each and every particle agent present in the simulated system at the specified iteration. These files serve as a basis for subsequent post- processing and graphical rendering ; producing a detailed 3D map of the simu- lated system. In addition, summary information about the overall concentrations of each chemical species or particle agent type is written as a row entry into a single summary report file (1007) or table at user specified iteration intervals.

Finally, information describing the computational performance of the simulation run is written upon completion of the simulation to a single runtime performance file (1008). Performance statistics recorded in this file include for example the total execution time or wall time, the ratio of time used in computation versus the idle and process to process communication time as well as other metrics of com- puter code performance.

Simulation Procedure-Serial Code and Distributed Computing Implemen- tation A general overview of the sequence of operations carried out during the course of a simulation is described in Figure 10. The simulation is started by launching the master code file (2001) which first processes data describing the computational environment. Information describing the in silico biological system to be simulated is then imported by execution of the subroutine inp_read (2002).

Another separate subroutine, get-pulse (2003), will import data describing The scheduling, type of stimuli and other details of external stimuli to be applied to the in silico system during the course of the simulation. Using information on the size of the problem and the computational environment available, the subroutine gridsizecaic (2004) is used to compute the size of the coordinate system for the simulated 3D system and reconcile the population sizes requested for each chemical species type with the user-specified limits on the overall volume fill. If the mode of execution specifies a parallel or distributed execution, then the user- requested number of parallel slave processes is launched. The overall simulation space is divided into sub-regions, and the locations of the borders of these sub- regions are computed and, together with the details of the problem, are communicated to each parallel slave process (2006). These problem characteristics include but are not limited to the population sizes and spatial distribution of particle agents of each chemical species type. Simulation methods and operations are carried out in parallel by the multiple slave processes. These slave processes are interrogated periodically to extract summary chemical assay results which are the collated and written to a report file (1007) by the subroutine reportchemcomp (2007).

Each of the slave processes carry out a number of operations and apply a num- ber of methods in parallel with one another. A slave process will start problem execution by initializing its region of the overall simulation space as defined by the spatial limits it has been communicated. The subroutine grid_init (3100) will compute the appropriate number of particle agents of each chemical species type and locate each agent in a specific part of the sub-region. Additional details of the initialization procedure are presented in Figure 11. Particle agents are first located on the cell membrane surface, if applicable, by subroutine rigid_membrane (3110). As a product of this subroutine and the set of methods used therein, a closed surface is computed in the shape and size specified by the user. Individual particle agents are then located on this geometrical closed sur- face until the surface is fully populated. This cell membrane is chemically active and the particle agents can participate in chemical reactions. This results in a chemically active outer shell for a given virtual cell enabling selective transport of agents of specific chemical species into and out of the cell. Particle agents des- tined for the intra-cellular space of the sub-region are positioned by subroutine init inner (3120) and those destined for the extra-cellular space are created by initouter (3130). A central particle agent data table (1001) is created in each sub-space and the location and other attributes of each agent are recorded as new entries. If the user has requested the simulation continue from a previous state then this exhaustive initialization of the main particle agent data table is by- passed. Additional details of this procedure are shown in Figure 12. Instead of initializing from an empty subspace, data stored in the snapshot fils (1006), re- corded for each slave process at the specified time step, are imported by the subroutine snap_inp_read (3200). This data is then used to produce the initial version of the main particle agent data table as well as inter-process transfer ar- rays and reset the required iteration counters. These transfer arrays are used to capture particle agent diffusion as will be described later in this specification.

The occurrence of chemical reactions in the simulation space is now addresses by executing the methods and operations in the subroutine sim slave chem (3300). A more detailed description of this step is presented in Figure 13. Each particle agent in the data table is now examined to determine if the local condi- tions are capable of supporting one or several chemical reactions. First a map of neighboring particles is constructed by the search method neighbourTbl (3310).

Then each particle agent is processed by the methods in subroutine state-build (3320) to determine the quantity and type of neighbouring agents. If the local neighborhood for an agent is populated in such a way that a chemical reaction might be possible, then the probability that said reaction actually occurs is evalu- ated and the reaction carried out if appropriate. This is accomplished by execut- ing the methods in the subroutine local react (3330).

The main processing steps carried out in the subroutine local react are listed in Figure 14. For each agent being processed, chemical reactions which involve the said agent's chemical species type are identified and labeled as candidate reactions. The number of particle agents of chemical species type present in the local neighborhood is then compared to the proportions required by each of the candidate chemical reactions. Candidate reactions which pass this test are en- abled. Should more than one reaction be enabled, then a first reaction is se- lected at random. The reaction proceeds if the probability of occurrence gener- ated exceeds the threshold probability specified in the p_react table (1005). New product particle agents generated by the chemical reaction are added as new en- tries to the main particle agent data table (1001) and particle agents consumed are labeled with a delete status and later removed. New product particle agents are positioned in the local neighborhood space in either a random fashion for regular reactions or in a pseudo-random fashion for directional reactions involv- ing receptor or transport molecules.

Following the execution of all candidate chemical reactions deemed probable, the main data table (1001) is updated. This updated particle agent population is then submitted to the subroutine sim slave_nav (3400) so that the displacement of each agent in space may be computed. Details of this set of methods are pre- sented in Figures 15,16 and 17. First a map of neighboring particles is re- constructed by the search method neighbourTbl (3310), this time using the navi- gation horizon distance to determine neighborhood limits. Once again, each par- ticle agent is processed by the methods in subroutine state build (3320) to de- termine the quantity and type of neighboring agents. The new location and velocity of each agent is computed by the methods in the subroutine posca ! c (3410) summarized diagrammatically in Figure 16. This subroutine uses a condi- tional probability distribution which is computed by the subroutine move_prob (3410). This conditional probability distribution is computed based on local con- ditions and a number of navigational cues which are summarized in Figure 17.

Using the map of neighboring particles, a subset of allowable moves which do not violate spatial constraints is identified. Based on the mean temperature field specified by the user, and based on transport properties of the particle agent type, the probability of remaining stationary is computed. The probability that any movement whatsoever will occur is computed as the complement. The desirabil- ity of making a specific allowable move is then computed based on specific navi- <BR> <BR> <BR> <BR> gational cues and scaled into conditional probability values. These probabilities are combined and rescaled to produce a single cumulative conditional probability distribution cumproba.

Once again data in the main particle agent data table is updated in the subrou- tine pos-Cale (3410). Elements of this updated table are then recorded to a de- tailed snapshot file (1006) by the subroutine recordpointpos (3600). A chemical assay is then computed for the slave sub-space by subroutine caicsubpop (3500) which conducts a census of the population size for each particle agent type. The results for each subspace are returned by each slave process and the results compiled to produce an overall particle agent population profile. This overall population profile is periodically written to a chemical summary report file (1007).

Since particle agents are mobile and travel through the simulation space, meth- ods have been implemented to transfer control of particle agent from one slave process to another. Subsequent to the initial iteration, the location of each parti- cle agent is tested to determine if the agent still resides in the spatial region as- signed to the current slave process. Should it no longer reside in the current subspace, the particle agent is moved from the main data table to a transfer ar- ray. This transfer array is then communicated to the appropriate neighboring slave process using a numberlof MPI functions (4001). In addition to this transfer function, the subroutine balance slave (4002) is used to counter discrepancies in computational load between slave processes resulting from chemical reaction and diffusion phenomena. The exchange of particle agents with neighboring slave processes as a result of diffusion as well as the generation and consump- tion of particle agents resulting from chemical reaction produce fluctuations in the population size in each slave process subspace. The subroutine balance slave (4002) applies a control method to shift the borders of the subspace regions in such a way as to produce similar particle agent population sizes for each slave process.

Example Results: Competitive Inhibition in Enzymatic Reactions.

To illustrate the results which can be obtained from a simulation run and to sub- stantiate the applicability of this approach to predicting the behavior of biochemi- cal pathways, a simple example is presented. A basic set of reversible enzyme- catalyzed reactions is described in Figure 20A. This reaction system involves an enzyme E, an inhibitor l, an intermediate product ES and a final product P. The behavior of this reaction system was simulated using two different combinations of initial substrate and initial inhibitor concentrations, [So] and [Io]. Since the vol- ume of the simulation space is held constant in this example, the concentration is directly represented by the number of particle objects of each species. In both simulated cases, no intermediate product ES or final product P were present at the outset. The evolution in the number of molecules of product P is presented under both sets of initial conditions in Figure 18 and clearly adheres to the char- acteristic shape of a progress curve expected from first order chemical reaction systems. As expected, altering the initial balance between the number of sub- strate particles and the number of inhibitor particles directly affects the rate of generation of product P. These progress curves are not entirely smooth however and clearly show the underlying stochastic nature of the simulation. As the num- ber of particles participating in the simulation is increased an averaging effect will lead to the generation of smoother curves familiar to the macro-scale world of solution chemistry.

Using the same basic reaction system, a number of additional simulations were performed. Simulations were conducted at 4 different levels of initial substrate concentration [So] or rather initial particle count for substrate S. In addition, each of these was repeated at two distinct levels of initial particle count for inhibitor 1.

Simulation results are presented in Figure 19 and clearly demonstrate that the data generated by the simulation can adhere to behavior described by the Michaelis-Menten equations. It is important to note however that the simulation is not in any way limited to reactions and conditions which correspond to this type of global behavior. In fact, experiments have shown how this behavior is not conserved at low concentrations or under conditions which deviate from the basic assumptions used by Michaelis and Menten in their derivation.

To further support the use of this approach in expioring and predicting the behav- ior of biochemical reaction systems, a number of simulations were conducted to demonstrate classical reversible enzyme inhibition. The basic types of reversible inhibition are competitive inhibition, uncompetitive inhibition, and noncompetitive inhibition. Competitive inhibitors are the most commonly found inhibitors in bio- chemistry and constitute a solid basis for the validation of this simulation ap- proach. A simple reaction system which displays competitive inhibition is pre- sented in Figure 20A. In such a system, the inhibitor I will bind only to the en- zyme E to form the complex El. Inhibitor will not bind to enzyme molecules which have already formed a complex with the substrate S. The behavior of such systems is best illustrated by plotting the reciprocal of the initial rate or 1/vo against the reciprocal of the initial substrate concentration or 1/ [So]. In a system which displays competitive inhibition, the value of the intercept on the vertical axis is conserved, that is 1Nmax is constant, for changing values of initial inhibitor concentration [lo]. The maximum reaction velocity Vmax is a well known metric used to characterize the kinetics of biochemical reactions. This metric is defined as the value of the initial velocity vo for a solution of enzyme E saturated in sub- strate S. Simulations were conducted using the invention described herein un- der a variety of initial substrate and inhibitor concentrations. As these simula- tions are stochastic rather than deterministic, repeat simulation runs were per- formed to mimic the averaging effect present in large-scale laboratory data. The results of individual runs are nonetheless presented in the figures. The results of these simulations were analyzed to determine the initial velocities vo obtained in each case. Both the initial velocity and the reciprocal of the initial velocity are plotted against the initial substrate concentration and it's reciprocal respectively in Figure 20B. Curves obtained for each of the 3 initial inhibitor concentrations are seen to cross on the y axis near the origin as one would expect. These simulation results indeed correspond to behavior expected of ciassicai competi- tive inhibition. Values obtained for the maximum reaction velocity Vmax converge to the same value as the initial substrate concentration tends to infinity.

A similar set of simulation experiments were carried out, this time using the reac- tion network presented in Figure 21A. This reaction network represents the ba- sic set of reactions involved in producing uncompetitive inhibition. Contrary to competitive inhibitors, uncompetitive inhibitors bind only to the enzyme-substrate complex ES, not to the free enzyme. In uncompetitive inhibition, Vmax is de- creased by the conversion of some of the enzyme molecules to the inactive form ESI. Another metric, Km, is also affected by uncompetitive inhibitors. The Michaelis constant Km is defined as the concentration of substrate present when the initial velocity vo is equal to one half the value of Vmax. Uncompetitive inhibi- tors also serve to decrease Km because the equilibria for the formation of ES is shifted by the binding of 1. Although both Vmax and Km are altered the ratio of Vmax to Km remains unchanged in uncompetitive inhibition. The result is a set of parallel lines, one for each initial level of inhibitor concentration [lo], on the double reciprocal plot. This is indeed the result obtained by the simulator as shown in Figure 21 B. Once again, the approach embodied in the present invention pro- duces results consistent with the known behavior of enzymatic reactions.

Finally, a reaction network which exhibits noncompetitive inhibition is presented in Figure 22A. Noncompetitive inhibitors can bind to free enzyme E or substrate S to form the inactive complexes EI or ESI respectively. This form of inhibition is characterized by an apparent decrease in Vmax (or increase in 1Nmax) while Km remains unchanged. This inhibition cannot be overcome by the addition of S and is relatively rare. As Km is unchanged, the resulting double reciprocal plot will consist of curves intersecting at the x axis to the left of the origin. Once again, simulations were conducted at 3 different levels of inhibitor 1. For each level of inhibitor, multiple levels of initial substrate concentration were used. In addition, for each combination of initial conditions, repeat simulations were per- formed to account for the stochastic nature of the simulation space. The initial velocities for each case were computed and used to prepare the Michaelis- Menten graphs shown in Figure 22B. As dictated by typical noncompetitive inhi- bition, the curves for each level of inhibitor concentration intersect on the x axis, recreating the expected behavior whereby Km remains unchanged.

Once again, it should be noted that these examples serve only to demonstrate that the simulation approach embodied in the present invention is quite capable of generating results which are fully consistent with certain well known aspects of the behavior of biochemical reaction networks. The capabilities, as well as the fields of application, of the cell simulator described herein are not in any way lim- ited to the simple examples presented here. Instead these examples serve only as demonstrations.

A Potential Application Problem Definition: The present example (Figure 23) shows two enzymatically catalyzed and re- versible metabolic reactions occurring within a cell (A+B>C & X+Y>Z). We as- sume that the presence of minimal and defined concentration levels for the prod- ucts C and Z are required for cell viability. The goal of this exercise is to maxi- mize the production of the metabolite C. In the interest of brevity, we have con- fined this example application to a pharmaco-kinetic problem in which only the concentration of small molecule metabolites that are exterior to the cell (A, B, X, Y &Z) may be manipulated experimentally or theoretically. It should be noted however that with the inclusion of genetic and protein engineering methodology, the existing application could easily incorporate the manipulation of protein spe- cific activity and concentration levels.

The reaction is regulated by the reactant X (a partial activator of E1), and the product Z (a partial inhibitor of E1). The equilibrating permeases EP1 and EP2 facilitate the diffusion of the reactants and products participating in these two re- actions back and forth across the cell membrane (EP1 : A, B, C, D, Y & Z EP2: X). To illustrate the utility of the present invention toward describing both the temporal and spatial aspects of cell physiology, we have localized the enzyme E1 to a specific region of the cell whereas E2 can diffuse throughout the cell.

This problem is defined by a number of physical and chemical parameters (ki- netic, thermodynamic and diffusion constants) as listed below. We assume herein that each of these parameters can be determined or approximated empiri- cally. We also assume that the shape of the cell is spherical and that its volume and membrane surface area are known, although this is not in any way a limita- tion imposed by the invention.

1. Concentrations of E1, E2, EP1 & EP2.

2. The thermodynamic and kinetic rate constants that define the be- havior of E1, E2, EP1 & EP2.

3. Starting concentrations of A, B, C, X, Y & Z both inside and outside the cell.

4. Diffusion rates of A, B, C, X, Y, Z, E1 & E2within the cell.

5. Diffusion rates of A, C, X, Y & Z outside the cell.

6. Diffusion rates of EP1 & EP2 within the cell membrane.

7. The position of E1 within the cell.

With the algorithm described herein, it is possible to recreate each of these pa- rameters by adjusting: 1) The size and shape of the reaction space, 2) the speed, number of particles and particle size of each chemical species, 3) their distribu- tion within the reaction space, and 4) the forward and reverse probabilities of a reaction occurring upon molecular collision.

Experimental Protocol: 1. All reactions are formatted into the'reaction manager'notation and entered as sub-reaction objects (303) into the pathway object list (301). This results in a collection of 48 sub-reactions (and their as- sociated forward & reverse probabilities) describing every step of each full reaction in mechanistic terms (Table 1).

2. The speed, concentration and distribution of each molecular spe- cies is adjusted and set at defined values to conform to empirically determined quantities.

3. The forward and reverse reaction probability values for each sub- reaction are adjusted and set to defined values in order to conform to empirically determined rates and thermodynamic constants de- scribing each overall reaction. This step can be carried out through manual input by the user or by taking advantage of an automated learning scheme.

Experimental Design: From Figure 23 it is clear that the rate at which C is synthesized can be in- creased by any combination of the following treatments: 1) increasing the con- centrations of reactants A and B therefore promoting the reaction in the forward direction, 2) increasing the concentration of the E1 activator X and finally 3) de- creasing the concentration of the E1 inhibitor Z. It remains unclear however what the optimal balance of concentrations exterior to the cell might be required in or- der to produce the maximal rate of production of species C without the loss of cell viability. Some of the complexity that underlies this problem may be de- scribed as follows : 1. The effectiveness of raising the external concentrations of A, B and X ( will be limited by the value of Vmax associated with their respec- tive permeases and enzymes.

2. Molecular species that use EP2 to transit back and forward across the membrane are defacto competitive inhibitors of each other.

3. The effectiveness raising the external concentration of the activator ( is counterbalanced by an increase in the production of the inhibi- tor Z.

4. Lowering the concentration of Y lowers the concentration of Z, however the latter cannot fall below the level required for cell viabil- ity.

5. In the interest of simplicity, this example has been presented as a closed system. In reality, the pools of some or all of the molecular species will be influenced by other reactions that either contribute or detract from their size or the activity of their constituents.

Clearly, these questions cannot be answered with conventional simulation tools and require the spatial and temporal discretization embodied in the present in- vention.

The simulation system and methods presented herein form a virtual environment which is envisioned as a predictive tool that will be capable of addressing many pragmatic scientific issues, including 1) the cellular consequences of gene muta- tion, 2) toxicity, 3) the direct visualization of cellular activity and 4) optimization of compound synthesis, 5) drug discovery and 6) genetic systems engineering to name just a few. Our early efforts will no doubt focus on single-celled organisms using these as proof-of-concept models that will generate the theoretical con- cepts and technologies for extending the virtual cell simulation to more compli- cated cell types and eventually, to multi-cellular organisms.

Although a few preferred embodiments have been shown and described, it will be appreciated by those skilled in the art that various change and modifications might be made without departing from the scope of the invention. The terms and expressions in the preceding specification have been used therein as terms of description and not of limitation, and there is no intention in the use of such terms and expressions of excluding equivalents of the features shown and described or portions thereof, it being recognized as the scope of the invention as defined and limited only by the claims that follow.

Table 1. Hypothetical Application: Reactions & Sub-reactions formatted in Reac- tion Manager notation 1 Transport : P (forward) P (reverse) 1) X*Epi > XEpi (ezsterior) Nl N2 XEp, >X+Epl (interior) N3 N4 2) A Ep2>AEp2 (exterior) Ns N6 AEP2 > A+EP2 (interior) N7 N8 3) B*Ep2 >BEp2 (exterior) N9 Nlo BEp2 >B+Ep2 (interior) Nu N12 Conversions : 4) A*Ei >AEi N13 N14 AEi*B>ABEi Nis Nig ABEt >CEi N17 Nig CE, > C+Ez Nl9 N20 7) X*El > EIX N21 N22 A*EIX >AEIX N23 N24 AEiX*B>ABEiX N25 N25 ABEIX > CE1X N27 N28 CEIX > C+EIX N29 N30 8) Z*EI > EIZ N31 N32 A*EIZ > AEzZ N33 N34 AEIZ*B > ABEzZ N35 N36 ABEIZ > CEIZ N37 N38 CEIZ > C+EIZ N39 N40 7) X*E2 > XE2 N41 N42 XE2*Y >XYE2 N43 N44 XYE2 >ZE2 N45 N46 ZE2 >Z+E2 N47 N48 Shown are the seven reactions and 48 sub-reactions that define the system. The reaction no- menclature is defined as follows: 1. Adjacent letters: an assembly of molecules 2. Plus sign : indicates that molecules or assemblies to either side are autonomous.

3. Asterisk [*]: read as contact or collides (ie. A contacts B).

4. Greater than [>]: read as"goes to".

Implicit in this notation is that the reverse of reaction is read from right to left. The likelihood that a reaction occurs in either direction is determined by the probabilities associated with the forward (PF) and reverse reaction (PR).

Cited references 1. E-cell web site (www. e-cell. org) 2. Tomita, M., Ventner, J. C., Hutchison, C. et al. "E-CELL : software envi- ronment for whole-cell simulation"Bioinformatics. 15 : 72-84 (1999).

3. Tomita, M."Whole-Cell Simulation : A Grand Challenge of the 21st Cen- tury"Trends. Biotech. 19 : 205-210 (2001).

4. McAdams, H. H. and Arkin, A. Simulation of prokaryotic genetic circuits Annu. Rev. Biophys. Biomol. Struct. 27 : 199-224 (1998).

5. Rao, C. V. and Arkin, A. P."Control Motifs for Intracellular Regulatory Net- works"Annual Rev. Biomed. Eng. 3: 391-419 (2001).

6. MEWG web site (http : //www. epa. gov/opptintr/metabolic/index.htm#about) 7. Edwards, J. S. and Palsson, B. O."the E. coli MG1655 in silico metabolic genotype its definition, characteristics and capabilities". Proc. Natl. Acad. Sci.

97: 5528-5533 (2000).

8. Palsson, B. O. "What Lies Beyond Bioinformatics"Nature Biotechnology 15: 304 (1997) 9. Tavazoie, S. , Hughes, J. D., Campbell, M. J. , Cho, R. J. , Church, G. M."Sys- tematic determination of genetic network architecture"Nat Genet. 22 (3): 281-285 (1999).

10. Johnson, J. M. and Church, G. M."Predicting ligand-binding function in families of bacterial receptors"Proc Natl Acad Sci U S A. 97: 3965-3970 (2000).

11. Hood, L."Computing life : the challenge ahead". IEEE Eng Med Biol Mag.

20: 20 (2000).

12. Silicon Cell project web site (http ://www. bio. vu. nl/ hwconf/Silicon/index. html) 13. Endy, D. and Brent, R."Modelling Cellular Behaviour"Nature 409: 391-395 (2001).

14. Varner, J. D."Large Scale prediction of Phenotyp@"Biotechnol. Bio@ng.

69: 664-678 (2000).

15. Fell, D. A. "Understanding the Control of Metabolism" (Portland, London) (1996).

16. Schuster, S., Dandekar, T. and Fell, D. A."Detection of elementary flux modes in biochemical networks : a promising tool for pathway analysis and meta- bolic engineering"Trends Biotechnol. 17 : 53-60 (1999).

17. Schuster, S. , Fell, D. A., Dandekar, T."A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks"Nat Biotechnol. 18 : 326-332 (2000).

18. Cornish-Bowden, A and Cardenas, M. L."Complex networks of interac- tions connect genes to phenotypes". Trends Biochem Sci. 26 : 463-465 (2001).

19. Cornish-Bowden, A, Cardenas, M. L. "From genome to cellular phenotype- -a role for metabolic flux analysis ? Nat Biotechnol. 18 (3): 267-268 (2000).

20. Mendes P. "Biochemistry by numbers: simulation of biochemical pathways with Gepasi3". Trends Biochem Sci. 22 (9): 361-363 (1997).

21. Mendes, P. and Kell, D. B."MEG (Model Extender for Gepasi): a program for the modelling of complex, heterogeneous, cellular systems. "Bioinformatics.

17 (3): 288-289 (2001).

22. Sauro, H. M. "SCAMP : a general-purpose simulator and metabolic control analysis program. "Comput. Appl. Biosci. 9 (4): 441-450 (1993).

23. Jarnac web site (http://members. tripod. co. uk/sauro/ Jarnac. htm) 24. Goryanin, I., Hodgman, T. C. and Selkov, D. "Mathematical simulation and analysis of cellular metabolism and regulation"Bioinformatics. 15: 749-758 (1999).

25. Wang, W. , Donini, O., Reyes, C. M., Kollman, P. A."Biomolecular simula- tions: recent developments in force fields, simulations of enzyme catalysis, pro- tein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions" Annu Rev Biophys Biomol Struct. 30: 211-243 (2001).

26. Briggs, J. M."Prediction of rates of encounter of a small molecule with a protein using a brownian dynamics method in a continuum solvent environment" Cell Mol Biol Lett. 6 (2B) : 533-534 (2001). 27. Garcia-Olivares, A. , Villarroel, M. , Marijuan, P. C. "Enzymes as molecular automata: a stochastic model of self-oscillatory glycolytic cycles in cellular me- tabolism". Biosystems. 56 (2-3) : 121-129 (2000).

28. Wurthner, J. U., Mukhopadhyay, A. K., Peimann, C. J."A cellular automaton model of cellular signal transduction"Comput Biol Med. 30 (1) : 1-21 (2000).