Whole-cell modeling framework in which biochemical dynamics impact aspects of cellular geometry

ARTICLE IN PRESS

Journal of Theoretical Biology 244 (2007) 154–166 www.elsevier.com/locate/yjtbi

Whole-cell modeling framework in which biochemical dynamics impact aspects of cellular geometry Ivan V. Surovsteva,b, Jeffrey J. Morganc, Paul A. Lindahla,b, a

Department of Chemistry, Texas A&M University, Spence and Ross Streets, P.O. Box 300012, College Station, TX 77843-3255, USA b Department of Biochemistry and Biophysics, Texas A&M University, P.O. Box 300012, College Station, TX 77843-3255, USA c Department of Mathematics, University of Houston, Houston TX, 77204-3008, USA Received 25 May 2006; accepted 20 July 2006 Available online 7 September 2006

Abstract A mathematical framework for modeling biological cells from a physicochemical perspective is described. Cells modeled within this framework consist of at least two regions, including a cytosolic volume encapsulated by a membrane surface. The cytosol is viewed as a well-stirred chemical reactor capable of changing volume while the membrane is assumed to be an oriented 2-D surface capable of changing surface area. Two physical properties of the cell, namely volume and surface area, are determined by (and determine) the reaction dynamics generated from a set of chemical reactions designed to be occurring in the cell. This framework allows the modeling of complex cellular behaviors, including self-replication. This capability is illustrated by constructing two self-replicating prototypical whole-cell models. One protocell was designed to be of minimal complexity; the other to incorporate a previously reported well-known mechanism of the eukaryotic cell cycle. In both cases, self-replicative behavior was achieved by seeking stable physically possible oscillations in concentrations and surface-to-volume ratio, and by synchronizing the period of such oscillations to the doubling of cytosolic volume and membrane surface area. Rather than being enforced externally or artificially, growth and division occur naturally as a consequence of the assumed chemical mechanism operating within the framework. r 2006 Elsevier Ltd. All rights reserved. Keywords: Systems biology; Cell cycle; Biokinetics; Interfacial reactions; Growth and division synchronization

1. Introduction Living cells are extraordinarily complex machines that exhibit remarkable behaviors including, most predominately, self-replication. Although much progress has been made in understanding the mechanism(s) underlying this behavior, further efforts on both experimental and theoretical fronts will be required to deepen our understanding of it. Towards this end are efforts to simulate global cellular behavior starting from the molecular level of detail (Palsson, 2000; Tomita, 2001; Kitano, 2002; Mori, 2004; You, 2004). Ideally, such whole-cell models would Corresponding author. Department of Chemistry, Texas A&M University, Spence and Ross Streets, P.O. Box 300012, College Station, TX 77843-3255, USA. Tel.: +1 979 845 0956; fax: +1 979 845 4719. E-mail address: [email protected] (P.A. Lindahl).

0022-5193/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2006.07.020

integrate massive amounts of biological information into compact, unambiguous and testable hypotheses which could be used to explore the entire panoply of cellular relationships in both healthy and diseased states (Ideker et al., 2001; Jansen, 2003; Lindon et al., 2004; Rao et al., 2005). Much remains before useful comprehensive mathematical models of cellular behavior with predictive power can be realized. Obviously, the vast number of unknown kinetic parameters (rate coefficients and concentrations) and ambiguities in chemical mechanisms occurring within a cell hinder progress, but there is also no consensus as to what modeling framework will ultimately prove most effective. By framework, we mean a group of general assumptions and procedures that can transform a physicochemical model of a cell into a set of mathematical expressions. Such a framework would be applicable to all

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

such models, regardless of the particular set of chemical reactions and components assumed. Given the complexity of biological cells, it seems intuitively obvious that cell models should be enormous (in terms of numbers of reactions and components), and indeed numerous groups are working on constructing ‘‘genome-scale’’ models. Palsson and co-workers have analysed very large metabolic systems at steady state in constant volume using constraint-based modeling and fluxbalance analysis, without decomposition of reaction fluxes to concentration dependence (Fo¨rster et al., 2003). Understanding the behavior of large metabolic models in terms of robustness, sensitivity, and regulation will also be required of cell models (Morohashi et al., 2002). Metabolic Control Analysis (Hofmeyr and Westerhoff, 2000) and Biochemical Systems Theory (Ni and Savageau, 1996) have been effective in analysing perturbation of systems from steady state. Tyson, Novak and others have focused on the cell division cycle, using the classical modeling approach based on deterministic kinetics, differential equations and continuous variables to generate chemical oscillations (Novak and Tyson, 2003; Qu et al., 2003; Tecarro et al., 2003). Arkin, McAdams and others have shown the importance of stochastic kinetic in modeling complex genetic systems (Arkin et al., 1998). Loew, Mendes, Tomita, Ortoleva and others have developed software environments (the Virtual Cell, Gepasi, E-cell, Karyote, etc.) in which large-scale userdefined biochemical systems involving multiple cellular compartments can be modeled without the need for users to fully understand the underlying mathematics (Mendes and Kell, 2001; Sayyed-Ahmad et al., 2003; Slepchenko et al., 2003; Takahashi et al., 2004). An excellent review on cell modeling has recently appeared (Ishii et al., 2004). During the process of self-replication, a ‘‘newborn’’ cell grows by importing nutrients and converting them into all cellular components. When some critical size is reached, the cell divides to form two ‘‘daughter cells’’ indistinguishable from the newborn parent. With unlimited nutrients, this process repeats itself ad infinitum, giving rise to a periodic cell cycle in which component concentrations oscillate in time and are coordinated to changes in cell geometry. Thus, modeling the cell cycle will involve connecting changes in chemical concentrations to changes in cellular geometry. Different types of concentration oscillations in biochemical system have been identified and kinetically modeled (Muller et al., 1998; De Wit, 1999; Zhdanov, 2000; Kurin-Csorgei et al., 2005; Madsen et al., 2005). Oscillatory behavior associated with the cell cycle in particular is coordinated with periodic changes of cell shape/geometry. In a previous study (Morgan et al., 2004), we developed a framework which allows cell growth and division to be modeled. However, periodicity within that framework was enforced by specified invariant periodic changes in cell geometry. Particular geometrical changes were required to occur in concert with cell growth. These assumptions regarding the time-dependent shape of the cell were

155

external (and thus insensitive) to any assumed chemical mechanism. In this paper, we describe a novel cell-modeling framework in which changes in particular aspects of cell geometry, namely cell volume and membrane surface area, arise naturally from the chemical dynamics inherent in an assumed chemical mechanism. When appropriate kinetic parameters and the remaining determinants of cell geometry are assumed, self-replicating behavior of cells can be modeled. We illustrate the novel features of this framework with two models: a minimal model of a hypothetical protocell, and a model of the eukaryotic cell cycle that is capable of autonomous synchronized oscillations of cellcycle effectors and of cell geometry. 2. Development of the modeling framework 2.1. Cellular compartments and components We consider an individual cell in a nutrient-rich aqueous environment surrounding the cell. The concentrations of environmental components remain constant, such that they collectively compose an isotonic solution (307 mM) with an osmotic pressure of penv ¼ 7.80 atm. The cell contains an aqueous cytosol with volume V1, and a non-aqueous plasma membrane with surface area S1 and it may also include additional compartments (e.g. organelles). In general, each lth compartment would have a volume Vl encapsulated by its respective surface membrane with area Sl (1plpL). A cell would consist of m chemical species C1,y,Cm. Species that exist in more than one region are viewed as different distinguishable species. For the ith reagent, ni denotes moles of i while Ci denotes concentration. Species in aqueous regions have volume concentrations C i ¼ ni =V l , whereas species in membrane regions have surface concentrations C i ¼ ni =S l . Each volume is a well-stirred chemical reactor, such that components are distributed uniformly and transportation effects can be neglected. Membranes are two-dimensional liquids with a third dimension that is infinitely small. Membrane components are oriented; they can be exposed to either face of the surface (e.g. membrane-bound receptors), or to both volumes on either side of the surface (e.g. permeases). Other details are given in Supplementary Material A. 2.2. Reaction rates of non-interfacial and interfacial reactions The m cellular components react with each other and with environmental nutrients according to s modeldependent reactions, which collectively constitute the chemical mechanism. Reaction j (1pjps) is described as m X i¼1

aij C i !

m X i¼1

bij C i ,

(1)

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

156

where aij and bij are stoichiometric coefficients of reactants and products, respectively. According to the law of massaction, the rate of the jth reaction is Ym a (2) Rj ¼ kj i¼1 C i ij , where kj is a rate coefficient. The rate of change (or flux) of the ith component in the lth compartment due to the jth chemical reactions within this compartment is 1 dni ¼ ðbij  aij ÞRj þ Fi V l dt

or

1 dni ¼ ðbij  aij ÞRj þ Fi , S l dt (3a)

where Fi denotes all other terms that contribute to the flux of component i. Eq. (3a) is appropriate when all reactants and products are located in the same region (Vl or Sl). For such non-interfacial reactions, the rate of the reaction and the chemical flux of Ci due to the reaction are determined with respect to the same region. The situation is different for interfacial reactions where one or more reactants or products are in different regions. Generally, an interfacial reaction will involve membrane components reacting with volume components on the surface. Interfacial reactions involving two volumes separated by a membrane are also possible. As shown in Supplementary Material B, the chemical flux of a surface component involved in an interfacial reaction is described by Eq. (3a) while analogous changes of a volume component is described by 1 dni ¼ wl 0 l ðbij  aij ÞRj þ Fi , V l dt

(3b)

where wl0 l is a geometric parameter defined as wl0 l ¼ Sl0 /Vl (i.e. the surface-to-volume ratio). The total chemical flux of Ci due to all s reactions, denoted as Fi, is determined by a weighted sum of the all reaction rates: Fi ¼

s X

mij Rj ,

(4a)

j¼1

where weights mij ¼ wl0 l (bijaij) if the ith component exists in a volume region and is involved in interfacial reactions— otherwise mij ¼ bijaij. As a result, a component ni will satisfy an equation of the form 1 dni ¼ Fi V l dt

or

1 dni ¼ Fi S l dt

dC i =dt ¼ F i . The situation is fundamentally different when reaction volume varies in time. In this case changes in volume must be taken into account (Aris, 1969), leading to the expressions dC i dC i ¼ F i  al C i or ¼ F i  bl C i . (5) dt dt Here al and bl are specific growth rates of volume Vl and surface Sl, respectively, defined by al ¼

1 dV l V l dt

2.3. Systems with changing volumes Within the framework, the chemical content of the cell is described through concentrations. In the case of constant volumes and surface areas (i.e. a non-growing cell), Eqs. (3) can be applied to calculate changes in concentrations as

bl ¼

1 dSl . S l dt

(6)

The Fi terms on the right-hand sides of Eq. (5) reflect changes in concentration due to chemical reactions, while the terms containing a and b represent changes due to volume or area changes (e.g. dilution due to cell growth). Changes in the number of molecules ni are only due to chemical reactions, not to dilution effects. Under steady-state conditions all reaction rates are constant, and all al and bl functions become constants in time. This allows Eq. (6) to be integrated, yielding V l ðtÞ ¼ V l ð0Þ expðal tÞ and Sl ðtÞ ¼ S l ð0Þ expðbl tÞ. Thus, for a cell that is simply growing at steady state, all volumes Vl grow exponentially, with al serving as the exponential growth rate-constant. Similarly, under these conditions, all surfaces Sl grow exponentially in accordance with exponential growth rate constant bl. Under non-steady-state conditions, al and bl are time-dependent and so Vl(t) and Sl(t) can change in complex ways; this is the situation required for modeling self-replication behavior, as illustrated below. 2.4. Mechanism of cell growth At this point, functions al and bl are undetermined. In other words, the volume of each region Vl and the area of each surface Sl are free to vary in time but there is no mechanism by which they can vary. For any volume, this deficiency is rectified by requiring that the osmotic pressure in the volume equals that of the environment at all times. As the simplest approximation, the ideal gas law can be used to calculate volume of lth compartment (R, T are gas constant and temperature, respectively), yielding the relationship

(4b)

for volume and surface components, respectively. In general, the total chemical flux Fi depends on all concentrations and reaction rates, and depending on the component, it may also depend on one or more geometrical parameters (w’s).

and

Vl ¼

RT p

ql X

ni .

(7)

i¼ql1 þ1

For surfaces, we assume that the entire area Sl must be occupied by components (e.g. phospholipids and proteins); i.e. there can be no solvent or empty space in membranes. As such, the area of a membrane surface will be determined by summing the areas of all components within it, Sl ¼

ql X

si ni

(8)

i¼ql1 þ1

with unit-areas si (surface area per one molecule). We will refer to this as the area-filling-tile assumption. Eqs. (7) and

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

(8) yields the conservation laws ql X

si C i ¼ 1

(9)

i¼ql1 þ1

for membrane components, and ql X

Ci ¼

i¼ql1 þ1

p ¼ C0 RT

(10)

for volume components. Eq. (10) indicates that the collective concentration of components in the volume must remain fixed and equal to that of the isotonic solution C0 at all times. The growth rates for volume and membrane of lth compartment could be calculated from Eqs. (7) and (8) by derivation: al ¼

ql RT X F i, p i¼q þ1

(11)

l1

bl ¼

ql X

si F i .

(12)

i¼ql1 þ1

Thus, the mathematical formulations for the volume of aqueous regions and the area of membranes (Eqs. (7) and (8)) allow us to connect particular geometrical properties of the cell with the chemical processes within it. Indeed, changes in volumes of all compartments in the cell are dictated by the rates of chemical reactions occurring in the cell. Geometric parameters wl 0 l can be calculated from the appropriate Sl0 and Vl values, or alternatively, from the appropriate al0 and bl, via the equation dwl 0 l ¼ wl 0 l ðbl 0  al Þ. (13) dt Eqs. (5)–(6) and (11)–(13) form a system of autonomous differential equations. These equations completely describe the dynamics of the cell within the framework: dC i =dt ¼ F i  a1 C i ; .. .

1pipq1 ;

dC i =dt ¼ F i  aL C i ; dC i =dt ¼ F i  b1 C i ; .. . dC i =dt ¼ F i  bL C i ;

qL1 þ 1pipqL ; qL þ 1pipqLþ1 ; .. . q2L1 þ 1pipm:

dSl =dt ¼ bl Sl ;

1plpL;

dV l =dt ¼ al V l ;

1plpL:

157

volume and surface. If the unit surface areas si and all reaction rates kj are known, then the system (14) will have a unique solution for each given set of initial data. However, we are interested in a much more complex problem than simply solving Eqs. (14). A cell model should obviously include a chemical mechanism in which all components of the cell are either directly or indirectly derived from nutrients. Such models should also exhibit the property of self-replication. Consider an idealized behavior of a self-replicating spherical cell (Fig. 1). A newborn cell would grow as an expanding sphere and then divide by constriction about its middle, yielding two newborn daughter cells having the same geometry as their newborn parent. The time for this cell cycle will be called t. The daughter cells would repeat this cycle such that all concentrations and the surface-tovolume ratio would oscillate with period t. For a cell requiring a single w function (derived from one compartment consisting of an outer membrane and the cytosol), cell growth would correspond to a decrease in w (since the surface-to-volume ratio declines with expansion) while constriction/division would correspond to an increase of w, ultimately returning to its initial value. In other words, a chemical model of a self-replicating cell should posses a stable periodic solution for Eqs. (13) and (14a). Ideally, the volume of the cytosol and area of the membrane (Eqs. (14b)) would double during period t, i.e. growth would be synchronized with division. Thus, our fundamental problem (P) can be stated as follows: Given a value t40, determine unit areas si and rate coefficients kj, so that there exists a solution Ci(t), Vl(t), Sl(t) to Eqs. (14) so that C i ðtÞ ¼ C i ð0Þ, V l ðtÞ2V l ð0Þ, S l ðtÞ ¼ 2Sl ð0Þ, and C i ðtÞ, Vl(t), Sl(t) contain physically realistic quantities (1pipm and 1plpL). The meaning of the phrase physically realistic will depend upon specific information associated with the chemical system as well as inherent physical constraints. Our strategy for solving problem (P) involves the use of an averaged system associated with the chemical concen-

(14a)

(14b)

2.5. Requirements of a self-replicating cell model For a given chemical mechanism of a cell, i.e. a list of reactions occurring within the cell, the set of Eqs. (14) determine the dynamics of the cell, including changes of

Fig. 1. Geometry of the assumed cell model. Dimensions of cell images were obtained from the simulations described in the text. Time (t) in min, cell radius (R, in mm), and intersphere separation (x, in mm) for images (A–D) are (0, 0.58, 0}, {38, 0.64, 0}, {67, 0.55, 0.41} and {100, 0.58, 0.58}, respectively. Expanded membrane region illustrates the assumed invariant packing density.

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

158

trations in Eqs. (14). First note that if we have a solution to (P), then the average values of al and bl are equal ln 2/t. The averaged geometric parameters wl 0 l can be denoted by wl 0 l (v¯ denotes the set of all independent wl 0 l ). Using these averaged values, we make a first approximation of Eqs. (14) by decoupling them into the systems dC i =dt ¼ F i ðvÞ  V l ðtÞ ¼ 2t=t V l ð0Þ

ln 2 Ci; t and

for 1pipm,

(15a)

S l ðtÞ ¼ 2t=t Sl ð0Þ for 1plpL. (15b)

We refer to Eqs. (15) as the averaged system associated with Eqs. (14) and a steady state solution of this system as an idealized steady state. Then we proceed as follows: (i) Determine realistic ranges for the chemical concentrations Ci, the unit areas si and the rate coefficients kj. (ii) For randomly selected, physically realistic chemical concentrations Ci and given unit-surface areas, solve the linear system (15) for the reaction constants kj. That is, find kj so that C ¼ Ci is an idealized steady state. This system will typically be underdetermined, and as a result, this process will serve to reduce the number of independent entries in ki. (iii) Determine whether the values obtained in part ii are physically realistic. If not, then repeat part (ii). (iv) Use the remaining parameters to steer the system towards a Hopf bifurcation to a family of nontrivial asymptotically time periodic solutions with period approximately t for the concentration portion of the averaged system. If the search fails, then return to part (ii). Use the method of continuation to identify a family of time periodic solutions to Eqs. (13) and (14a) associated with the values from the averaged system with period t and doubling time t. (v) Use the independent parameters to steer through this family to find a solution to problem (P). If the search fails, then return to part (ii). If the phrase ‘‘physically realistic’’ is quantified , then it is possible to repose steps (ii)–(iii) above as a minimization problem. We illustrate this process on a simple system in the section below.

each identical in size and composition to the original newborn parent cell. The framework operates with the volume and surface of the cell as variables; these aspects are necessary but not sufficient determinants of cell geometry. Thus, if all other parameters required to completely define cell geometry are specified—call this the geometry shell—the entire cell shape can be derived. The only remaining parameters needed to completely define the geometric changes depicted in Fig. 1 would be the radius of the parent and two daughter cells, and the inter-sphere distance (R and 2x, respectively, in Fig. 1C). As shown below, these last two parameters can be obtained as output from a chemical mechanism modeled using our approach. We assumed a newborn cell consisting of cytosol (with volume Vcyt) encapsulated by a membrane (with surface Smem) and located in an aqueous-based environment containing a single nutrient N (Fig. 1A). In this system, there is a single w ¼ Smem/Vcyt. The aqueous-based cytosol contains three symbolic components M, L, and P, which loosely represent the {metabolome+genome}, phospholipid-ome, and proteome of the cell, respectively. The membrane consists exclusively of two components, Pmem and Lmem, with unit-surface areas sP and sL . The five reactions constituting the assumed chemical mechanism are shown in Fig. 2. Pmem catalyzes the import of N and conversion into M. P catalyses, to second order, the conversion of M into L. Stoichiometric coefficients m1 and m2 reflect the numbers of M molecules required to generate one P and L, respectively; however, these coefficients are not taken as exponents in rate expressions. The set of ODEs (Eqs. (13) and (14)) describing the dynamics for this particular cell model is presented in Eqs. (16) of Table 1. These equations fully determine the dynamics of the cell once all parameters of the model are specified (Supplementary Material C). 3.2. Search for appropriate periodic solutions There is no guarantee that any set of selected parameter values would result in periodic behavior, or that such behavior, if achieved, would be physically possible. We investigated by a random search of the parameter space whether the model described by Eqs. (16) was capable of periodic behavior under any set of concentrations and rate

3. Construction of a minimal whole-cell model R2

3.1. Defining features We first illustrate this framework by seeking a model with a minimal number of reactions and components capable of satisfying requirements formulated as (P) above. Such a mechanism might be considered to be a minimal whole-cell (MWC) model of the hypothetical self-replicating protocell illustrated in Fig. 1. Here a newborn spherical cell first expands (during a growth phase of the cell division cycle) and then divides into two spherical daughter cells,

P

R4

Pmem

m1

N

R1

M m2 R3

L

R5

Lmem

Fig. 2. Kinetic mechanism of the MWC model. Rates for reactions (R1–R5) are given in Table 1. Dotted lines indicate catalytic influence.

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

159

Table 1 Chemical reactions, rates and ODEs for the MWC and TN-MWC models MWC Model d M=d t ¼ wR1  m1 R2  m2 R3  aM d P=d t ¼ R2  wR4  aP d L=d t ¼ R3  wR5  aL d Pmem =d t ¼ R4  b Pmem d Lmem =d t ¼ R5  b Lmem d w=d t ¼ wðb  aÞ d V cyt =d t ¼ aV cyt

(16)

TN-MWC model d M=d t ¼ wR1  R2  R3  R6  R8  R12  aM d P=d t ¼ R2  wR4  aP d L=d t ¼ R3  wR5  a L (18) d Pmem =d t ¼ R4  bPmem d Lmem =d t ¼ R5  bLmem d Clnb=d t ¼ wR6  wR7  aClnb d Cdh1=d t ¼ R12 þ R14  R13  wR15  aCdh1 d Cdh1a=d t ¼ R13  R14  aCdh1a d Cdc20a=d t ¼ R9  wR11  aCdc20a d Cdc20=d t ¼ R8  R9  wR10  aCdc20 d w=d t ¼ wðb  aÞ d V cyt =d t ¼ aV cyt

R1 ¼ k1 Pmem N R2 ¼ k 2 M R3 ¼ k 3 P 2 M R4 ¼ k 4 P R5 ¼ k 5 L a ¼ ðwR1  ðm1  1ÞR2  ðm2  1ÞR3  wR4  wR5 Þ=C 0 b ¼ sP R4 þ sL R5 R1 R2 R3 R4 R5 R6 R7

¼ k1 Pmem N ¼ k2 M ¼ k3 M ¼ k4 Cdc20 P ¼ k5 Cdc20 L ¼ k6 P2mem M ¼ k7 Cdh1a Clnb ðClnb=J 8 Þn R8 ¼ k 8 M 1 þ ðClnb=J 8 Þn Cdc20 R9 ¼ k9 Clnb J 9 þ Cdc20 R10 ¼ k10 Cdc20 R11 ¼ k11 Cdc20a R12 ¼ k12 M Cdh1 R13 ¼ k13 Cdc20a J 13 þ Cdh1 Cdh1a R14 ¼ k14 Clnb J 14 þ Cdh1a R15 ¼ k15 Cdh1 a ¼ wðR1  R4  R5  R7  R10  R11  R15 Þ=C 0 b ¼ sP R4 þ sL R5

coefficients. For each set of parameter values randomly chosen within specified limits, the Jacobian matrix for Eqs. (16) was calculated using the resulting k-values and selected concentrations. Eigenvalues of this matrix were assessed to determine whether a Hopf bifurcation would give rise to a stable periodic solution (Jordan and Smith, 1999). Locating Hopf bifurcations is only one method for identifying periodic solutions, and our use of it undoubtedly biased our search. Nevertheless, this search method executed rapidly, and was used to assess the probability of the model to possess a stable periodic solution. For the model given above, the results of 105 trials revealed 1917 successful attempts, each indicated by a point in the 3D plot of Fig. 3. Numerical integration of Eqs. (16) with parameters from this random search showed that solutions approached a stable periodic orbit, as expected. There are other considerations besides periodicity in selecting appropriate solutions. First, the ability to mimic growth and division of real cells requires that the growth rate (given by averaged a and b values) matches the division rate (given by the frequency of oscillations in w). In other words, the period of w oscillations must be synchronized with the time of volume doubling (note that periodic behavior of the w does not guarantee such synchronization). Secondly, the amplitude of w oscillations (i.e.

wmax–wmin) should remain within physically reasonable limits. Such limits depend on the assumed cell geometry. For a spherical newborn cell that enlarges and then divides, this suggests that wmax should not exceed the sum of w for two spheres, each with Vcyt ¼ 1/2. Similarly, wmin should be greater than that of one sphere with V cyt ¼ 2. The parameters found using the random search routine had to be adjusted to satisfy these additional requirements. Both subcritical and supercritical types of Hopf-bifurcations were revealed by exploring the dependence of the system’s behavior upon increment change of parameters. The latter type is characterized by ‘‘soft’’ excitation of oscillations which allowed the amplitudes to be adjusted. Set of parameters exhibiting this behavior were adjusted to satisfy additional requirements (regarding the period and amplitude of w(t)) described above (see Supplementary Material D). The resulting system matched growth and division rate, and had an appropriate range of changes in surface-to-volume ratio. As such, it represents a physically possible periodic solution. Since equations for concentrations and w form a complete subset in Eqs. (16) and do not depend explicitly upon Vcyt or Smem, there is a freedom in choice of an initial point of the cell cycle (t ¼ 0), at which Vcyt and Smem are assigning to their initial values. For a variety of cell shapes,

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

160

2.0

45000

Vcyt , Smem

(m1,m2.χav,ρ,P0,L0)

Vcyt

1.5

Smem 1.0 5.3

P0

χ, μm-1

30000

15000

5.0 4.7

12000

0

0

(a)

25

50

75

100

v

4000 6.7

m2

0

Fig. 3. Search routine used to find states capable of periodic behavior. A set of parameters [m1, m2, w0, r, P0, L0] was chosen at random within defined limited ranges. Rate coefficients were obtained by solving the system assuming idealized steady-state conditions. The corresponding Jacobian matrix was constructed and evaluated. Each point (m2, w0, P0) in the 3D plot represents a solution with a positive eigenvalue, projected from the 6D parameter space [m1, m2, w0, r, P0, L0]. Parameter r ¼ (m1  Pmem)/(m2  Lmem) reflects the density of Pmem in the membrane. Further details are provided in Supplementary Materials.

w will decline as the cell grows and increase as it divides. This constrained our selection of w0 to the region where dw/dto0. The particular point selected within this region will influence the geometry of the cell during its celldivision cycle (since Vcyt and Smem indeed depend upon the variables of the subset). To finalize the model and generate cell shape dynamics, this initial point should be specified. 3.3. Dynamics and analysis of the MWC model The dynamics of the final MWC model which resulted from these adjustments are shown in Fig. 4. Vcyt, Smem and w are plotted between t ¼ 0 and t (panel 4a). As required, volume and surface area double with each cell division cycle, but neither grows purely exponentially nor purely linearly. We find it interesting that Woldringh et al. (1993) reported nonexponential growth in budding yeasts. To the best of our knowledge, there is no mathematical model to account for this. However, a population of such cells in a culture, viewed over many cell cycle periods, would appear to grow exponentially. The geometry shell assumed with our model allows for various scenarios in which two fully overlapping spheres of radius R0 and inter-sphere separation 2x ¼ 0 could grow and divide to form two separate cells with radius R0 and x ¼ R0 . The volume and area are related to geometry according to the relationships   4 x2 V cyt ¼ pR3 þ 2px R2  3 3   x and S mem ¼ 4pR2 1 þ . ð17Þ R

0.50

R

0.25

growth

division

x 0.00 (b)

0

25

50

75

100

M Concentrations

3.3 χa

R, x, μm

8000

0.0

P L Pmem Lmem 0

(c)

25

50 Time, min

75

100

Fig. 4. Dynamics of the MWC model (Eqs. (16), Table 1): (a) Cytosol volume, membrane area, and w during a cell cycle period. Both Vcyt and Smem are plotted relative to initial values, (b) sphere radius R and intersphere separation distance x were calculated from numerical w data assuming that cell geometry can be described as two spheres that overlap completely at t ¼ 0 min, then grow and separate to yield two completely separated spheres at t ¼ 100 min (Eqs. (17)). The file MWC.mov is a Quicktime animation of the cell cycle geometry using the R and x data, (c) concentrations of each component in the model. Minimum:maximum concentrations (in mM) for M, P, and L were 194,000:214,000, 4,800:5,600; 90,000:107,000, respectively. Similar ranges (in mm mM) for Pmem and Lmem were 3.745:37.32 and 8115.7:8116.4, respectively.

For each [Smem, Vcyt] pair, these equations were used to obtain corresponding values for R and x (Fig. 4b). Ideally, growth and division phases are expected to be distinct, such that R would initially increase and x would remain at 0, followed by a phase in which x would increase to R0 and R would return to R0. Indeed, plotted profiles have two distinct phases. The growth phase, evident between 0 and 40 min, shows increasing values of R and insignificant changes in x. The division phase, located between 40 and 100 min, exhibits declining values of R and increasing x. Within the growth phase, R and x are complex numbers, with imaginary part small relative to real part. Only the real part of these numbers are plotted. This situation arises because Smem within this phase is slightly less than what

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

would be physically possible for a sphere of volume Vcyt. Back-calculating values for Smem and Vcyt from the real parts of these complex R and x values yielded discrepancies with actual Smem and Vcyt values of less than 6%. We regard this as physically acceptable, in that real membranes stretch by comparable amounts (Evans and Skalak, 1979). An animation of these computationally generated geometry changes is provided as Supplementary Material. Time profiles for all chemicals concentrations within the minimal protocell are shown in Fig. 4c. Concentrations oscillate in a cell-cycle-dependent manner, with P and L in phase, and M essentially counter-phase (keeping total cytosol concentration constant, as required). Pmem and Lmem are similarly counter-phased, as required by Eq. (9). A requirement of oscillatory chemical systems is that they be antagonistic, with dominance by one component at one time causing the dominance of another component at a subsequent time, which eventually leads to re-establishing the dominance of the original component ad infinitum. The MWC model can be simplistically viewed in those terms. The major flow is through the pathway N-M-L-Lmem (i.e. the system expends most of its efforts to build the phospholipid membrane), while the M-P-Pmem pathway is a minor regulating offshoot. The catalytic action of P and Pmem on the two reactions of the major pathway can be considered as valves that control flow, with the one controlling the downstream reaction (L-Lmem) more sensitive than that controlling the feed into the system (because of the higher order-dependence). At high P, the more sensitive valve opens, causing Lmem to increase. This causes Pmem to decrease by dilution. This attenuates the major flow, causing Lmem to decline and Pmem to rise. A more accurate analysis of the behavior of this system would require evaluating the influence of each term of the set of ODEs at each point along the cell cycle. Within the constraint of this five-component skeleton model, with Pmem always catalyzing the import of N into the cytosol, we evaluated various related mechanisms differing only by the component assumed to be a catalyst for other steps and by the order-dependence of that catalytic influence (Fig. 5 and Table 2). This suggests that periodic solutions are obtained only when P or Pmem catalyzes the M-L reaction, and that they do so with an order-dependence greater than 1. Simply including a second-order dependence at any arbitrary step in the network (e.g. Table 2, model B) is not sufficient. We

Fig. 5. Different types of chemical networks evaluated via the random search routine described in the text. See Table 2 for results.

161

Table 2 Chemical mechanisms (Fig. 5), evaluated for Hopf bifurcation Order dependence a A B C D E F G H I J K L

1 2 1 2 1 1 1 2 1 1 1 1

b 0 0 0 0 0 0 1 1 2 3 0 0

c 0 0 1 1 2 3 0 0 0 0 0 0

Attempts d 0 0 0 0 0 0 0 0 0 0 2 0

e

Total

0 0 0 0 0 0 0 0 0 0 0 2

7

10 107 107 107 107 106 107 107 105 105 107 107

Positive

%

0 0 0 0 2345 720 0 0 1917 2986 0 0

0 0 0 0 0.024 0.072 0 0 1.92 2.99 0 0

suspect that connecting the ‘‘P-branch’’ of the network to the ‘‘L-branch’’ with a catalytic influence is important, as this might create the antagonism needed for periodic behavior. We also fixed m1 and m2 stoichiometric coefficients to 1, a potential simplification, but stable physically possible periodic behaviour (SPPPB) was not forthcoming. There are undoubtedly more complex mechanisms that will give rise to SPPPB, but these were not explored. Although we cannot prove it mathematically, the models shown in Fig. 2 (with P or Pmem catalyzing the M-L reaction) appear to be the simplest chemical mechanisms exhibiting the periodic dynamics (in both concentration and geometry) required for self-replicating protocells having distinct and synchronized growth and division phases. 4. Cell model incorporating the Tyson–Novak cell-cycle mechanism In real eukaryotic cells, cell cycle events are temporally divided into G1, S, G2, and M phases. Passage from one phase to the next is controlled by the periodic ordered activation and inactivation of cyclin-dependent kinases (Cdks) (Murray, 2004). In proliferating cells these events are synchronized with cell growth. Synchronization is achieved through cell size checkpoints, mechanistically undefined controls which allow cells to delay cell cycle progress until they attain a certain critical size (Polymenis and Schmidt, 1999; Jorgensen and Tyers, 2004). A number of mathematical models and their non-linear analyses have been developed to understand the complexity of the cell cycle regulating network (Tyson and Novak, 2001; Qu et al., 2003) and these have afforded insight into this most fundamental biological process (for review see Ingolia and Murray, 2004; Jorgensen and Tyers, 2004). In kinetic mechanistic models of the cell cycle, cell-size synchronization is achieved by an assumed influence of the exponentially increasing cell size (or mass) on the rates

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

162

of certain regulating processes. This controls the progress of the cell from a stable state associated with G1 to either another stable state (Tyson and Novak, 2001; Chen et al., 2004) or to a limit cycle (Qu et al., 2003) associated with S–G2–M phases. At the end of a cell cycle simulation, defined by the attainment of some threshold condition, the governing cell size parameter is reset to half of its threshold value which returns the system to G1. This logical rule drives concentrations oscillations coupled with periodic cell size change, but it is external to the considered mechanism (i.e. it does not follow from the assumed chemical reactions). From a mathematical perspective this rule is tantamount to a forcing periodic term which actually drives oscillations in the model. In our previous cell modeling framework, we specified all aspects of geometry except cell size, and this was sufficiently severe as to enforce oscillatory behavior regardless of the assumed chemical mechanism (Morgan et al., 2004). What is required is a means of connecting the size of a cell to the chemical dynamics (changes in concentrations of components of a reaction network) occurring within it. Inspired by the behavior of our MWC model, we wondered whether this reset rule could be eliminated if the Tyson–Novak (TN) cell cycle mechanism (Tyson and Novak, 2001) was incorporated into the MWC framework. This was our motivation for designing the Tyson–Novak minimal whole cell (TN-MWC) model described below. 4.1. Design of the TN-MWC model For a complete description of the TN model, we refer readers to the original paper (Tyson and Novak, 2001). We employed a simplified version of this model which retains its essential behavior. At the core of the TN mechanism is a

R2

P

R6

N

R1

R4

hysteresis loop due to the antagonism between the Clnb:Cdk complex and the anaphase promoting complex (APC); APC destroys cyclins including Clnb, an activator of Cdk, while Clnb:Cdk inhibits APC activity by inactivating Cdh1, an auxiliary regulatory protein of APC necessary for cyclin degradation. Activation of Cdh1 is promoted indirectly by Cdc20, another regulatory protein of APC. Both forms of APC (with Cdc20 or Cdh1 associated) target cyclins for degradation. Cdc20 is synthesized in the S–G2–M phase, so it is assumed to be turned on by the Clnb:Cdk complex in a size-dependent manner and through a Hill-type function with nX2 (Tyson and Novak, 2001). The simplified version employed here is shown in Fig. 6 (red and blue portions). Forms of APC are indicated in red and are designed by the associated regulatory proteins (Cdc20 and Cdh1). This version exhibited stable periodic behavior (with t70 min) using previously reported parameters and the reset rule. Stable oscillatory behavior was abolished when the reset rule was removed. Next, we modified the MWC mechanism such that it could grow but not divide. To achieve this, we removed the feature of the MWC model most directly responsible for its oscillatory behavior; namely, the catalytic influence of P on the M-L reaction yielding the green portion in Fig. 6 (as in model A of Table 2). Using the idealized steady state (Eqs. (15), with t ¼ 70 min) we located kinetic parameters which produced a doubling time close to the period obtained of the TN-model described above. The two modified models were then combined as shown in Fig. 6 (red, blue and green portions). All cell cycle components were assumed to be synthesized from M and three components (Clnb:Cdk, Cdc20 and Cdc20a) were assumed to be specifically degraded. We omitted the

Pmem

R7

Clnb:Cdk

R10 R9

R8

R12

R3

L R5

Lmem

R11

Cdc20a

Cdc20

M

R13 Cdh1

Cdh1a R14

R15

Fig. 6. Kinetic mechanism of the TM-MWC model. Rates for reactions (R1–R15) are given in Table 1. Dotted lines indicate catalytic influence; doubleheaded arrows indicate second-order influence. Blue reactions represent a simplified version of the TN mechanism. Red components represent APC bound to various regulatory subunits in unactivated and activated states. Green reactions are from the MWC model. Brown catalytic influences were added to afford oscillatory behavior.

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

4.2. Dynamics and analysis of the TN-MWC model

Vcyt, Smem

2.0 Vcyt

1.5

Smem 1.0

χ, μm-1

3.3 3.0 2.7 0

(a)

25

R, x, μm

1.0

50

75

100

R

0.5

division

growth x

0.0 0

Concentration, μM

(b)

2

50

75

100

Cdh1a

Cdc 20a Clnb:Cdk 0

(c)

25

1

0

163

20

40

60 Time, min

80

100

Fig. 7. Dynamics of the TN-MWC Model (Eqs. (18), Table 1): Panels (a) and (b) are equivalent to those in Fig. 4(a) and (b). The file TN_MWC.mov is a Quicktime animation of the cell cycle geometry for this cell, (c) concentrations of each component in the model. Kinetic parameters for the simulations were k1 ¼ 0.0009, k2 ¼ 0.002, k3 ¼ 0.5  103, k4 ¼ 2.5, k5 ¼ 3000, k6 ¼ 0.32  1012, k7 ¼ 100, k8 ¼ 0.8042  106, J8 ¼ 0.15, n ¼ 2k9 ¼ 1, J9 ¼ 0.001, k10 ¼ k11 ¼ 3, k12 ¼ 0:16  106 , k13 ¼ 10, k14 ¼ 35, J 13 ¼ J 14 ¼ 0:04, k15 ¼ 4, sLmem ¼ 0:7200  106 , sPmem ¼ 0:72  104 , N ¼ 20000, C 0 ¼ 307000.

previous dependence of rate on size (in reactions 8 and 14 in Fig. 6) along with the corresponding reset rule. Instead, we assumed three additional catalytic dependencies (brown arrows in Fig. 6). Specifically, Pmem was assumed to be a second order catalyst for the synthesis of Clnb:Cdk from M, rationalized as connecting the oscillator engine with cell size through w and an interfacial reaction (similar to the P-branch influence on L production in the original MWC model, as discussed above). Cdc20 was allowed to catalyze the two reactions that increase membrane surface area (L-Lmem and P-Pmem) and hence to promote the division phase of the cell cycle, providing a feedback from the cell cycle engine to the rest of the cell. The corresponding set of ODEs which define the TN-MWC model are presented in Table 2. After minor adjustment, the system exhibited t ¼ 94 min and the stable periodic behavior shown in Fig. 7.

Time-dependent plots of volume, surface area and surface-to-volume ratio for the TN-MWC model are shown in Fig. 7a. In contrast to the MWC model, Vcyt increases smoothly in a near exponential manner while Smem increases more slowly for the first 50 min of the cell cycle, and then abruptly begins to grow rapidly thereafter. This causes the w function to have a more realistic sawtooth shape, rather than the simple sinusoidal shape generated using the MWC model. Corresponding values of R and x vs. time are given in Fig. 7b; they indicate that for the TN-MWC model, the process of cell division commences abruptly at 50 min. This seems to be more realistic, in that the proportion of the cell cycle period dedicated to growth in real cells is generally longer than the division phase and the division phase begins abruptly. An animation showing the growth and division of the TNMWC model is available as Supplementary Material. Notice the abrupt onset of the division phase relative, in comparison to the analogous behavior of the MWC model. Concentrations of cell cycle effectors vs. time are shown in Fig. 7c. Overall concentration dynamics are similar to those reported by Tyson and Novak (2001). The initial Clnb:Cdk concentration is kept low by elevated levels of Cdh1a, as is characteristic of the G1 phase. There is a slight difference in that the Cdh1a level increases slightly in G1 whereas in the original TN-model this concentration is almost invariant. The cell grows (Vcyt increasing as w declines) in this period. This is followed by an abrupt burst of Clnb:Cdk activity and a corresponding drop in Cdh1a, as is characteristic of S-phase. This triggers the increase in Cdc20a, as is characteristic of the M-phase. Finally, Cdc20a levels decline and Cdh1a abruptly increases, returning the system to G1. The burst in Cdc20 promotes synthesis of the cell membrane, as can be seen from surface area and surface-to-volume ratio time profiles. Phase trajectory of the TN-MWC model is presented as a projection onto the [Cdc20a, Clnb:Cdk] plane in Fig. 8a and onto the [w, Clnb:Cdk] plane in Fig. 8b. During this entire process, Vcyt and Smem double in synchrony with the period of w as achieved by adjustment of kinetic parameters. All aspects of this behavior, including concentration oscillations and abrupt triggering phenomena synchronized with growth and division, are fully determined by the assumed chemical mechanism of Fig. 7 without the need for a reset rule. 5. Discussion 5.1. Modeling framework Biological cells are chemical machines in which behavior is determined by complex molecular level interactions. Developing a comprehensive physicochemical model of the cell, which could mechanistically explain (and thus predict) behavior at the cellular level, will be an enormous

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

164

Clnb:Cdk, μM

0.4

M

S

0.2

G1

0.0

Cell Cycle Start

0.0

0.5

1.0 Cdc 20a, μM

(a)

1.5

Clnb:Cdk, μM

0.4 M

0.2 S G1

0.0 Cell Cycle Start

2.7 (b)

3.0

3.3

χ, μm-1

Fig. 8. Phase portrait of the TN-MWC model (Eqs. (18), Table 1). Arrows indicates Start of the cell cycle and flow directions.

challenge. Our objectives were to develop a modeling framework that might contribute incrementally towards this goal and to illustrate the properties of this framework by designing two whole-cell models. The framework consists of a set of simple (almost physically intuitive) assumptions that allow modeling certain limited aspects of real cells. These assumptions can be applied to any chemical reaction mechanism that satisfies the requirements of a cell model (namely that all components of the cell be derived from nutrients imported from the environment). The novelty of the framework is its ability to connect chemical dynamics, arising from the assumed chemical mechanism, to two important physical attributes of cell geometry, namely cell volume and surface area. The connection is made through accurate consideration of interfacial reactions. To satisfy the conservation of matter, this consideration entails a dependence of reaction fluxes upon a geometric parameter of the cell, namely the w-function, the surface-to-volume ratio. Surface area and volume can be calculated because certain assumptions are made regarding the physical properties of membranes and aqueous solutions (namely, the area-filling-tile assumption for membranes and the

osmotic pressure invariance for solutions). As the framework develops in the future, these simple assumptions could be replaced by more accurate expressions, perhaps including a non-ideal osmotic pressure equation or an expression that takes into account the ability of real membranes to stretch. Nor do our current expressions include terms that describe forces from cell walls, which in bacteria allow the cytosol to exist at high osmotic pressures. Indeed the type of protocell models that can be treated by the current framework is exceedingly simple relative to actual self-replicating cells. A framework capable of modeling real cells comprehensively would include spatial dimensions that would allow properties such as chromosome segregation, cytoskeletal movement, and concentration heterogeneities. However, such frameworks will take years of development and are beyond our interests here. The simplicity of our current framework allows for rapid numerical integrations and bifurcation analyses, as well as our focus on the relationship between the structure of a reaction network and its resulting dynamical behavior. Another important advantage of our framework is that the chemical and physical dynamics of the cell models constructed within it are governed exclusively by standard physical chemical laws and an assumed chemical network and kinetic parameters. This is achieved by treating volume and surface area as variables which affect (and are affected by) the other variables of the model. Once a geometry shell is assumed, models built within this framework are selfcontained and fully determined by initial values; no nonphysical reset rules or procedures external to the system need be applied. The difficulty of achieving desired self-replicating behavior using our framework is actually another advantage because it limits possible solutions. Perhaps the major problem in computational systems biology is that chemical systems are generally underdetermined such that a near infinite number of solutions exist. Three aspects of these systems contribute to this problem, including unknown rate-coefficients, component concentrations, and reaction mechanisms (which can be transformed into a set of rate law expressions and a stoichiometric matrix). Experimentally determining these parameters and relationships is arduous work even for simple systems involving only few reactions and components. Relative to frameworks that assume steady-state behavior and constant volumes, our framework severely reduces the number of solutions to those yielding stable, physically possible periodic behavior. 5.2. The minimal whole cell model The MWC model developed here appears to be the simplest hypothetical chemical mechanism which satisfies certain minimal formal requirements of self-replicating cells (see Problem (P) above), and can thus be viewed as a hypothetical protocell. The mechanism used to drive selfreplication in real cells is undoubtedly far more complex

ARTICLE IN PRESS I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

but is also poorly understood. By examining sizes of 6000 gene-deletion strains of Saccharomyces cerevisiae, Jorgensen et al. (2002) have systematically identified genes of 500 abnormally sized strains, indicating that a complex network is involved in the synchronization of the growth and division. At the heart of the problem is how cell growth is coordinated with cell division. There is general consensus that coordination is established internally, and that it is related to attaining a critical cell size or mass (Mitchison, 2003). However, this begs the question of how cells ‘‘know’’ their own size (Jorgensen and Tyers, 2004). Several different hypotheses have emerged to account for the sensitivity of intracellular reactions to cell size (Koch, 2002) and one of them, size-dependence through surface to volume ratio changes, is realized in our MWC model. Thus, despite the simplicity of this model, it may help understand essential aspects of the self-replicative mechanism of real self-replicating cells. One such aspect may be the inclusion of interfacial reactions. Such reactions may be required to synchronize growth with division in cell models. 5.3. The TN-MWC model The TN-MWC model represents a significant advance in modeling the cell cycle, because it allows cell size to be determined naturally from the physicochemical properties of the system, rather than being determined by an artificial/ phenomenological parameter that is external to the physico-chemical system and is used to force the system to oscillate. The TN-MWC model exhibits concentrations dynamics similar to those exhibited by the original TN model but in this case, dynamics are fully determined by the assumed chemical mechanism. Thus, the observed behavior is an internal essential feature of the mechanism. Because cell size is included within the models developed within our framework, mechanisms of cell-size checkpoints and their relationship to other cell cycle processes can now be explored, by developing different competing models within this framework and comparing their behavior to experimental results. The ability to merge two independent models within our framework is also important. Building a comprehensive molecular-level whole-cell model will inevitably involve merging independent models, but how (or whether) this can be done remains uncertain. Our ability to merge the TN and MWC models within the context of the framework developed here suggests that this framework might facilitate such mergers. 5.4. Conclusions The cell-modeling framework developed here is simple but it nevertheless has unique capabilities that make it useful for modeling certain aspects of growth and division processes in cells. These capabilities were illustrated by constructing a minimal symbolic self-replicating cell model

165

and a cell model that incorporates a previously developed cell-cycle mechanism. Synchronization of growth and division in our models was achieved because interfacial reactions were included. We believe that such reactions play critical roles in the mechanisms by which real cells sense their size and synchronize growth and division. Acknowledgements This study was supported by the Robert A. Welch Foundation (A-1170). Appendix A. Supplementary Materials Supplementary data associated with this article can be found in the online version at doi:10.1016/j.jtbi.2006. 07.020.

References Aris, R., 1969. Elementary Chemical Reactor Analysis. Prentice-Hall, Englewood Cliffs, NJ. Arkin, A., Ross, J., McAdams, H.H., 1998. Stochastic kinetic analysis of developmental pathway bifurcation in phage lambda infected Escherichia coli cells. Genetics 149, 1633–1648. Chen, K.C., Calzone, L., Csikasz-Nagy, A., Cross, F.R., Novak, B., Tyson, J.J., 2004. Integrative analysis of cell cycle control in budding yeast. Mol. Biol. Cell. 15, 3841–3862. De Wit, A., 1999. Spatial patterns and spatiotemporal dynamics in chemical systems. Adv. Chem. Phys. 109, 435–513. Evans, E.A., Skalak, R., 1979. Mechanics and Thermodynamics of Biomembranes. CRC Press, Boca Raton, FL. Fo¨rster, J., Famili, I., Fu, P., Palsson, B.O., Nielson, J., 2003. Genomescale reconstruction of the Saccharomyces cerevisiae metabolic network. Genome Res. 13, 244–253. Hofmeyr, J.H.S., Westerhoff, H.V., 2000. Building the cellular puzzle: control in multi-level reaction networks. J. Theor. Biol. 208, 261–285. Ideker, T., Galitski, T., Hood, L., 2001. A new approach to decoding life: systems biology. Annu. Rev. Genom. Hum. Genet. 2, 343–372. Ingolia, N.T., Murray, A.W., 2004. The ups and downs of modeling the cell cycle. Curr. Biol. 14, R771–R777. Ishii, N., Robert, M., Nakayama, Y., Kanai, A., Tomita, M., 2004. Toward large-scale modeling of the microbial cell for computer simulation. J. Biotechnol. 113, 281–294. Jansen, R.C., 2003. Studying complex biological systems using multifactorial perturbation. Nat. Rev. Genet. 4, 145–151. Jordan, D.W., Smith, P., 1999. Nonlinear Ordinary Differential Equations, third ed. Oxford University Press, Oxford. Jorgensen, P., Tyers, M., 2004. How cells coordinate growth and division. Curr. Biol. 14, R1014–R1027. Jorgensen, P., Nishikawa, J.L., Breitkreutz, B.J., Tyers, M., 2002. Systematic identification of pathways that couple cell growth and division in yeast. Science 297, 395–400. Kitano, H., 2002. Systems biology: a brief overview. Science 295, 1662–1664. Koch, A.L., 2002. Control of the bacterial cell cycle by cytoplasmic growth. Crit. Rev. Microbiol. 28, 61–77. Kurin-Csorgei, K., Epstein, I.R., Orban, M., 2005. Systematic design of chemical oscillators using complexation and precipitation equilibria. Nature 433, 139–142. Lindon, J.C., Holmes, E., Nicholson, J.K., 2004. Metabolomics and its role in drug development and disease diagnosis. Exp. Rev. Mol. Diagn. 4, 189–199.

ARTICLE IN PRESS 166

I.V. Surovstev et al. / Journal of Theoretical Biology 244 (2007) 154–166

Madsen, M.F., Dano, S., Sorensen, P.G., 2005. On the mechanisms of glycolytic oscillations in yeast. FEBS J. 272, 2648–2660. Mendes, P., Kell, D.B., 2001. MEG (model extender for Gepasi): a program for the modelling of complex, heterogeneous, cellular systems. Bioinformatics 17, 288–289. Mitchison, J.M., 2003. Growth during the cell cycle. Int. Rev. Cytol. 226, 165–258. Mori, H., 2004. From the sequence to cell modeling: comprehensive functional genomics in E. coli. J. Biochem. Mol. Biol. 37, 83–92. Morgan, J.J., Surovtsev, I.V., Lindahl, P.A., 2004. A framework for whole-cell mathematical modeling. J. Theor. Biol. 231, 581–596. Morohashi, M., Winn, A.E., Borisuk, M.T., Bolouri, H., Doyle, J.C., Kitano, H., 2002. Robustness as a measure of plausibility in models of biochemical networks. J. Theor. Biol. 216, 19–30. Muller, S.C., Mair, T., Steinbock, O., 1998. Traveling waves in yeast extract and in cultures of Dictyostelium discoideum. Biophys. Chem. 72, 37–47. Murray, A.W., 2004. Recyling the cell cycle: cyclins revisited. Cell 116, 221–234. Ni, T.C., Savageau, M.A., 1996. Model assessment and refinement using strategies from biochemical systems theory: application to metabolism in human red blood cells. J. Theor. Biol. 179, 329–368. Novak, B., Tyson, J.J., 2003. Modelling the controls of the eukaryotic cell cycle. Biochem. Soc. Trans. 31, 1526–1529. Polymenis, M., Schmidt, E.V., 1999. Coordination of cell growth with cell division. Curr. Opin. Genet. Dev. 9, 76–80. Palsson, B.O., 2000. The challenges of in silico biology. Nat. Biotechnol. 18, 1147–1150.

Qu, Z.L., MacLellan, W.R., Weiss, J.N., 2003. Dynamics of the cell cycle: checkpoints, sizers, and timers. Biophys. J. 85, 3600–3611. Rao, B.M., Lauffenburger, D.A., Wittrup, K.D., 2005. Integrating celllevel kinetic modeling into the design of engineered protein therapeutics. Nat. Biotech. 23, 191–194. Sayyed-Ahmad, A., Tuncay, K., Ortoleva, P.J., 2003. Toward automated cell model development through information theory. J. Phys. Chem. 107, 10554–10565. Slepchenko, B.M., Schaff, J.C., Macara, I., Loew, L.M., 2003. Quantitative cell biology with the virtual cell. Trends Cell Biol. 13, 570–576. Takahashi, K., Kaizu, K., Hu, B., Tomita, M., 2004. A multi-algorithm, multi-timescale method for cell simulation. Bioinformatics 20, 538–546. Tecarro, E.S., Obeyesekere, M.N., Auchmuty, G., 2003. Mathematical analysis of a 3-variable cell cycle model. Nonlinear Anal. Real World Appl. 4, 87–107. Tomita, M., 2001. Whole-cell simulation: a grand challenge of the 21st century. Trends Biotech. 19, 205–210. Tyson, J.J., Novak, B., 2001. Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions. J. Theor. Biol. 210, 249–263. Woldringh, C.L., Huls, P.G., Vischer, N.O., 1993. Volume growth of daughter and parent cells during the cell cycle of Saccharomyces cerevisiae a/alpha as determined by image cytometry. J. Bacteriol. 175, 3174–3181. You, L.C., 2004. Toward computational systems biology. Cell Biochem. Biophys. 40, 167–184. Zhdanov, V.P., 2000. Model of oscillatory pattern in cells: autocatalysis and transport via the cell membrane. Phys. Chem. Chem. Phys. 2, 5268–5270.