Linda Petzold and Yang Cao University of California Santa Barbara Daniel Gillespie

1. Multiscale Stochastic Simulation Algorithm with Stochastic Partial Equilibrium Assumption for Chemically Reacting Systems Linda Petzold and Yang Cao University of California Santa Barbara Daniel Gillespie Caltech / UCSB (consultant)

2. The Multiscale Problem in Systems Biology In the heat-shock response in E. Coli, an estimated 20 - 30 sigma-32 molecules per cell play a key role in sensing the folding state of the cell and in regulating the production of heat shock proteins. The system cannot be simulated at the fully stochastic level due to: • Multiple time scales (stiffness) • The presence of exceedingly large numbers of molecules that must be accounted for in SSA Khammash et al.

3. Computational Models of Chemical Reaction Systems • Discrete and stochastic - Finest scale of representation for well stirred molecules. Exact description via Stochastic Simulation Algorithm (SSA) due to Gillespie. We Are Here • Continuous and stochastic - The Langevin regime. Valid under certain conditions. Described by Stochastic Differential Equations (SDE). • Continuous and deterministic - The chemical rate equations. Described by ordinary differential equations (ODE). Valid under further assumptions.

4. Stochastic Simulation Algorithm • Well-stirred mixture • N molecular species • Constant temperature, fixed volume • M reaction channels • Dynamical state where is the number of molecules in the system

5. Stochastic Simulation Algorithm • Propensity function the probability, given , that one reaction will occur somewhere inside in the next infinitesimal time interval • When that reaction occurs, it changes the state. The amount by which changes is given by the change in the number of molecules produced by one reaction • is a jump Markov process

6. Stochastic Simulation Algorithm • Draw two independent samples and from and take the smallest integer satisfying • Update X

7. Hybrid Methods • Slow reactions involving species present in small numbers are simulated by SSA • Reactions where all constituents present with large populations are simulated by reaction-rate equations Haseltine and Rawlings, 2002 Matteyses and Simmons, 2002 What to do about other scenarios???

8. A Closer Look at the Multiscale Challenges for Heat Shock Response   The total “concentration” of 32 is 30-100 per cell But the “concentration” of free 32 is .01-.05 per cell DNAK RNAP  RNAP DnaK FtsH

9. Stochastic Partial Equilibrium Approximation • In deterministic simulation of chemical systems, the partial equilibrium approximation assumes that the fast reactions are always in equilibrium. These fast reactions are thus treated as algebraic constraints. • In stochastic simulation, the states keep changing. The stochastic partial equilibrium approximation is based on the assumption that the distributions of the fast species remain unchanged by the fast reactions. • Partition and re-index reactions so that are fast reactions and are slow reactions. Partition the species so that are fast species (any species involved in at least one fast reaction), and are slow species • The virtual system with only fast reactions and fast species is assumed to be in stochastic partial equilibrium, on the time scale of the slow reactions

10. Foundations of SPEA • Let the probability that, given , no slow reaction fires during . • Define . • We will not simulate the fast reactions. Instead we take attime as a random variable. The probability that one slowreaction will occur in is given by where mean over only.

11. Foundations of SPEA (continued) • Let be so small that in the next at most one slow reaction could fire. • Then (1) • Solving (1) we obtain (2) • Define the next slow reaction density function as the probability that, given , the next slow reaction will occur in and will be an reaction. Then (2) leads to (3) • Now apply the stochastic partial equilibrium approximation. Assume the virtual system is at partial equilibrium and is small compared to . Then can be taken as . • Thus

12. The Multiscale Stochastic Simulation Algorithm Given initial time , initial state and final simulation time , • Compute the partial equilibrium state for the fast reaction channels. This involves solving a nonlinear system arising from the equilibrium approximation and the local conservation laws for 2. For , calculate and . (This is the most difficult task) 3. Generate two random numbers and from . The time for the next slow reaction to fire is given by , where . The index of the next slow reaction is the smallest integer satisfying . 4. If , stop. Otherwise, update .

13. Down-Shifting Method The MSSA method captures the full distribution information for the slow species, but computes only the mean accurately for the fast species. However, the distribution information of the fast species can be easily recovered by taking a few small time steps whenever this information is needed. This ‘down-shifting’ works because the system is a Markov process!

14. Heat Shock Response Model • Stochastic Model involves 28 species and 61 chemical reactions. This is a moderate-sized model. • CPU time on 1.46 Ghz Pentium IV Linux workstation for one SSA simulation is 90 seconds. • CPU time for 10,000 simulations is more than 10 days. • 12 fast reactions were chosen for the SPEA. The fast reactions were identified from a single SSA simulation to be the ones that fired most frequently. These 12 reactions fire 99% of the total number of times for all reaction channels. • CPU time for the multiscale SSA • Without down-shifting: 3 hours for 10,000 runs • With down-shifting: 4 hours for 10,000 runs

15. MSSA Heat Shock Response Model Results Histogram (10,000 samples) of a slow species mRNA(DNAK) (left) and a fast species DNAK with (right) and without downshifting method (middle) solved by the original SSA method (purple solid line with ‘o’) and the MSSA method (red dashed line with ‘+’), for the HSR model. fast species with down shifting one slow species one fast species