Stochastic simulation, also known as kinetic Monte Carlo, is a numerical procedure for determining the dynamics of a continuous time Markov process.
A continuous time Markov process is a memoryless stochastic process that is used to describe all sorts of physical and chemical systems, including
biological systems. The solution of a Markov process includes two equivalent parts: a time-dependent probability distribution of states in state
space or trajectories of the state moving in state space over time. By computing an ensemble of trajectories, one can generate the distribution. Or,
by computing the distribution, one can sample it to obtain trajectories. If the state space of the Markov process is discrete, it is called a jump
Markov process and its probability distribution is governed by the Master equation. If the state space is continuous, it is called a continuous
Markov process and its distribution is governed by the Fokker-Planck equation. For all but the most trivial systems, the Master equation is
analytically unsolvable. For large dimensional systems (ie. many chemical species), the Fokker-Planck equation is impractical to solve. Stochastic
simulation is a way to generate trajectories of a Markov process and then to compute the distribution of all possible trajectories, effectively
sidestepping the problems with solving either the Master or Fokker-Planck equations. Of course, for some systems, stochastic simulation can be as
impractical, motivating the usage of hybrid stochastic simulation methods.

Biological systems can be described as a system of biochemical reactions, including processes such as metabolism, signal transduction, and
gene expression. The reactions may involve bond-breaking (as in the classic definition of a chemical reaction) or may only involve strong
non-covalent interactions between molecules, such as hydrogen bonding, salt bridges, or hydrophobic 'forces'. Either way, one may measure the
thermodynamics and kinetics of these reactions and quantitatively describe their rates. Traditionally, the dynamics of a system of reactions is
described using ordinary differential equations (ODEs). With ODEs, the trajectory of the state moves deterministically through time to some
attractor (ie. starting from the same initial condition, the trajectory will always be the same and no two trajectories may ever cross). However,
through both analysis of theory and experimental observation, it has been noted that ODEs do not adequately describe the dynamics of biological
systems. Instead, we describe the system as the above-mentioned Markov process. Now, given the same initial condition, different
trajectories of the state will occur and trajectories may cross paths in time. This seemingly simple change results in a large amount of new
phenomena, including populations of cells with differing phenotypes, cells which spontaneously change their phenotype, and oscillations which only
occur because of the 'noise' or stochasticity.

It is no surprise that there is a connection between jump Markov, continuous Markov, and deterministic processes. In fact, using theory, one can
state the necessary conditions in order to approximate a jump Markov process as a continuous one and, likewise, a continuous Markov process as a
deterministic one. By measuring the validity of these conditions, while the simulation is running , one can approximate only the
part of the system that meets these conditions. Our recent research shows that, after partitioning the system into multiple different mathematical
representations, one can self-consistently merge these disparate descriptions back together again. Transitions in the jump Markov process can be
converted to differential Jump equations (a type of stochastic differential equation). The effects of reactions described as a continuous Markov
process may also be described using stochastic differential equations. Likewise, the deterministic subset may be described using ODEs. All of these
equations are coupled and may be simultaneously numerically integrated forward in time using a variety of well-characterized stochastic numerical
integrators. It is this last part that separates our hybrid methods from others. We can characterize the error of the solution using the theory of
Wiener-driven SDEs and their numerical solution. Other hybrid methods sometimes rely on quite heuristic approximations which we can adeptly avoid.

Perhaps most importantly, these hybrid stochastic simulation methods can beorders of magnitude faster than the original
stochastic simulation algorithm. The exact speed up depends on the number of reactions classified as either jump Markov, continuous Markov, or
deterministic processes.

Perhaps most importantly, these hybrid stochastic simulation methods can be

Even when using hybrid stochastic
algorithms, simulations of realistic (ie. very large) biochemical
networks require serious computing power. In order to compute an
accurate probability distribution of the system dynamics, numerous
independent trajectories must be computed. A single trajectory of the
system dynamics may be computed in a short period of time (milliseconds
to minutes), but an accurate distribution requires at least 5,000
trajectories. Furthermore, it's highly desirable to simulate multiple
networks, each time changing one or more parameters, to determine
either the effect of that parameter on the network (sensitivity analysis)
or find a desired behavior (global optimization). Supercomputers are
quickly becoming common tools to the research community. We believe
their user-friendliness and widespread usage will only increase as in silico design yields more and
more fruitful results. Hy3S is designed to be used on cluster-based
supercomputers (not necessarily vector processors) with the
accompanying network and storage infrastructure, but the code will
function well enough on a fast personal computer.

The NetCDF
library and format is specifically designed to store huge amounts of
data in a platform-independent, self-describing optimized binary file.
Simulations of large networks produce an incredible amount of solution
data, requiring a data format that not only efficiently stores the
data, but also quickly retrieves it for analysis. NetCDF files may be
read and write in direct access mode, enabling access to hyperslabs of
non-contiguous blocks of data and reducing the memory requirements for
analysis. For example, even though a solution data file may contain
multiple trajectories of a system of multiple chemical species over
many time points, it is easy to read in only the data sampled from a
specific timepoint or about a certain chemical species. Instead of
reading in gigabytes of data, you can pick and choose which data to
read requiring much less time and using much less memory. The data may
be read from any operating system or architecture that has a C or
Fortran77/90/95 compiler. The API is available in additional
programming languages, including C++, Perl, Python, Java, MATLAB, and
Ruby.

SBML, or Systems Biology Markup Language, is a standard XML format for describing biochemical networks. Hy3S will eventually support the import and export of SBML files. We are waiting for the SBO classification terms in SBML v3.

SBML, or Systems Biology Markup Language, is a standard XML format for describing biochemical networks. Hy3S will eventually support the import and export of SBML files. We are waiting for the SBO classification terms in SBML v3.

To the best of our knowledge, this is the only stochastic simulation software that offers

- a fast hybrid algorithm,

- an easy to use GUI (usable on Windows, Linux, or Mac),

- parallelized code (MPI / Fortran95),

- useful features, such as cell replication events, gamma-distributed events, system perturbations, parameter variation, and initial condition variation
- the GPL license (LGPL is also available if you ask)
- developers who actually design biological systems themselves and
who will continue to add features useful to their own research

We're interested in creating simulation and analysis tools that enable computational biologists and engineers to
both study natural biological systems and design new synthetic ones. One goal is the ability to quickly simulate the stochastic dynamics of
*any* system of chemical or biochemical reactions, including non-linear ones with 'fast/continuous', 'fast/discrete', and
'slow/discrete' reactions. Another goal is to use those simulation results in a variety of design tools, such as global optimization and
sensitivity analysis. Of course, advances in the first goal will lead to faster results in the latter one. The mathematics of random dynamical
systems is relatively immature. Another goal is an advancement of the theory behind the numerical methods. The field is moving fast so stay tuned!