Converting Macromolecular Regulatory Models from Deterministic to Stochastic Formulation Pengyuan Wang, Ranjit Randhawa, Clifford A. Shaffer, Yang Cao, and William T. Baumann Virginia Tech, Blacksburg VA

The Fundamental Goal of Molecular Cell Biology

The Cell Cycle

Cell Cycle Control Mechanism

Modeling Techniques  One method: Use ODEs that describe the rate at which each protein concentration changes  Protein A degrades protein B: 

d

[ B ]  

c

[ A ]

dt

… with initial condition [

A

](0) = A 0 .

Parameter

c

determines the rate of degradation.

Sometimes modelers use “creative” rate laws to approximate subsystems

Simulation: Budding Yeast Cell Cycle

Expermental Data 1 2 3 4 wild type (daughter)

clb1 clb2 clb1 clb2

1X

GAL-CLB2 clb5 clb6

5 6

clb5 clb6 GAL-CLB5 sic1

7

sic1 GAL-SIC1

8

hct1

9 10

sic1 hct1 sic1 GAL-CLB5

first cycle second cycle mass at birth 0.71 0.71 mass at SBF 50% 1.07 (71’) 1.07 0.65 0.73 0.61 0.66 0.80 0.73 0.71 0.71 0.52 1.10 1.07 (65’) 0.93 1.00 (73’) 1.07 1.08 No SBF 0.74 Table 6. Properties of

clb, sic1,

and

hct1

mutants mass at DNA repl. 1.15 (84’) mass at bud ini. 1.15 (84’) mass at division 1.64 (146’) T G1 (min) changed parameter 84 Comments CT 146 min (time of occurrence of event) 1.16 1.16 No mit

Surana 1991 Table 1, G2 arrest.

1.19 1.30 (99’) 0.92 0.82 (37’) 1.38 1.17 0.72 0.73 No repl 1.19 1.17 (80’) 0.96 1.06 (83’) 1.17 1.18 No bud 0.76 1.50 1.70 (146’) 1.41 1.52 (146’) 1.86 1.69 No mit 1.20 105 99 73 38 94 82 k' s,b2 = 0 k" s,b2 = 0 k' s,b2 = 0.1 k" s,b2 = 0 k' s,b5 = 0 k" s,b5 = 0 k' s,b5 = 0.1 k" s,b5 = 0 k' s,c1 = 0 k" s,c1 = 0 k' s,c1 = 0.1 k" s,c1 = 0 k" d,b2 = 0.01 k' s,c1 = 0 k" d,b2 = 0.01 k' s,b5 = 0.1 k" s,b5 = 0 k' s,c1 = 0

Surana 1993 Fig 4, 1X GAL-CLB2 is OK, 4X GAL-CLB2 (or 1X GAL-CLB2db) causes telophase arrest. Schwob 1993 Fig 4, DNA repl begins 30 min after SBF activation. Schwob 1993 Fig 6, DNA repl concurrent with SBF activation in both GAL-CLB5 and GAL-CLB5db. Schneider 1996 Fig 4, sic1



uncouples S phase from budding. Verma 1997 Fig3B, Nugroho & Mendenhall 1994 Fig 2, most cells are viable. Schwab 1997 Fig 2, viable, size like WT, Clb2 level high throughout the cycle. Visintin 1997, telophase arrest. Schwob 1994 Fig 7C, inviable.

First cycle OK, DNA repl advanced; but pre-repl complexes cannot form and cell dies after the first cycle.

Putting it Together

Chen/Tyson Budding Yeast Model  Contains over 30 ODEs, some nonlinear.

 Events can cause concentrations to be reset.

 About 140 rate constant parameters  Most are unavailable from experiment and must set by the modeler

Fundamental Activities of the Modeler     Collect information  Search literature (databases), Lab notebooks Define/modify models  A user interface problem Run simulations  Equation solvers (ODEs, PDEs, deterministic, stochastic) Compare simulation results to experimental data  Analysis

Modeling Process

Stochastic Simulation Motivation  ODE-based (deterministic) models cannot explain behaviors introduced by random nature of the system.

  Variations in mass of division Variations in time of events  Behavior of small numbers (RNA, DNA)  Differences in gross outcomes

Gillespie’s Stochastic Simulation Algorithm (SSA)   There is a population for each chemical species There is a “propensity” for each reaction, in part determined by population  Each reaction changes population for associated species  Loop:  Pick next reaction (random, propensity)  Update populations, propensities  Slow, there are approximations to speed it up

Question  Given an existing deterministic model, how do we convert it to a formulation capable of stochastic simulation?

  Can this be automated?

Is there a fundamental difference in representation?

 SSA is known to be CPU-intensive. How much computation resource is really needed to simulate the converted model stochastically?

Relation between the Two Formulations  In common: both models describe the same reaction network.

 Difference: the reaction rate equation is replaced by a propensity function describing how likely that the reaction will fire in next unit time.

 Connection: although they have different physical meanings, propensity function shares the same expression as corresponding reaction rate equation (written in number of molecules).

 Caveat: except for the “creative” rate laws

Missing Information

 Usually ODE models are written in terms of normalized concentrations.

 Thus they need to be converted to models in terms of number of molecules (population).

 Some information is missing  Characteristic concentration  Explicit definition of units  Volume of the container.

Conversion

 The relation between normalized concentration, real concentration and population of a species:

How Units are Used in the Model  Every parameter and species is assigned the correct unit, scaling factors.

 The conversion algorithm follows units to convert the model.

The Challenge  Assigning correct units to species and parameters is difficult because all the species, parameters, and reactions are connected by the whole reaction network.

 Once the modeler is forced to provide the “complete” specification, the conversion can be automated  Caveats:  “Creative” rate laws  Events

Events Need Extra Care

/*deterministic events*/ If (A>threshold) Then {event is triggered}.

/*stochastic events*/ If (Athreshold) Then {event is triggered; minimum=A}.

(Here “>” means rising above a threshold) (we ask for A truly rising from a low value, not happening to rise by oscillation.)  Except for events, all other parts of the model are automatically converted by JigCell.

Conversion Tool

  Part of the JigCell modeling suite Automatically checks unit consistency inside the model   Every two quantities (a parameter, a species, or the result of a sub-expression) connected by + or - in the rate law equation must have same units.

All species whose values are changed by the same reaction must have the same units.

 The unit of the result from the rate law equation must be equal to the unit of the reaction rate.

The Tool: Entering the Data

The Tool: Error Checking

The Tool: Error Correction

The Tool: Results – Reactions

The Tool: Results – Unit Types

Simulation Experiments: Setup

   Model:  A simplified cell cycle model  A full-sized budding yeast cell cycle model* Data:   38 of 45 species in full-sized model use realistic characteristic concentration found in the literature.

Cell volume is set to 50

fL.

Simulator:  StochKit, a C++ stochastic simulator integrated into JigCell, running SSA.

Distribution of Species on Converted Simplified Model  Ensemble result of 10,000 simulations at 200 minutes simulation time.



Simulations on the Converted Full sized Model

The same model (except events) can be simulated either deterministically or stochastically  The interesting cases are where they do not agree

Mass at Birth, Full-sized Model  Mean = 1.20, CV = 2.96%. (Compared with 1.21 from deterministic simulation)

Variance of Mass at Birth vs. Simulation Time vs. Population

Simulation Times Model Simplified Full-sized Wall Stochastic Time Total Avg./run 145 3862 12305 382267 1.23

38.2

Determini stic Time 0.029

0.311

 Even a single run of the stochastic simulation takes much more time than the deterministic simulation.

 Parallel computing is needed and feasible.

Effect of Random Number Generators SPRNG random()

Conclusions

 Improved support for the conversion process  The JigCell conversion tool  Deterministic and stochastic formulations are not fundamentally different  Deterministic modelers like to take short cuts  Real experience with stochastic simulations on meaningful models  Events  Runtimes  Approximation results

Future Work  Initial conditions distribution  Truly growing volume:  Our previous model had growing mass but fixed volume, which is not realistic  Change to growing volume will change the reaction rate (propensity function)  Simulations on mutants of particular interest