Computer models of bacterial cells: from generalized coarsegrained to genome-specific modular models

We discuss a modular modelling framework to rapidly develop mathematical models of bacterial cells that would explicitly link genomic details to cell physiology and population response. An initial step in this approach is the development of a coarse-grained model, describing pseudo-chemical interactions between lumped species. A hybrid model of interest can then be constructed by embedding genome-specific detail for a particular cellular subsystem (e.g. central metabolism), called here a module, into the coarse-grained model. Specifically, a new strategy for sensitivity analysis of the cell division limit cycle is introduced to identify which pseudo-molecular processes should be delumped to implement a particular biological function in a growing cell (e.g. ethanol overproduction or pathogen viability). To illustrate the modeling principles and highlight computational challenges, the Cornell coarsegrained model of Escherichia coli B/r-A is used to benchmark the proposed framework.


Introduction
A release of the 1000 th microbial genome is expected within the next few years [1]. This data will become a central bioinformatic resource in modern biology by providing access to thousands of accurate metabolic reconstructions of completely annotated genomes [1], genome-scale reaction networks, and detailed stoichiometric matrices and models [2]. Despite their dominance and fundamental importance, intrinsically static metabolic reconstructions and stoichiometric models are, by themselves, insufficient to explicitly relate dynamic nonlinear physiologic responses to static linear genomes. The predictive capability of stoichiometric models is limited to the calculation of instant phenotype snapshots under fixed defined media conditions. Therefore, such models cannot capture dynamic changes in metabolite concentrations, protein machinery, cell geometry, etc.
We develop a complementary modeling approach to relate genomic data to dynamic intracellular processes: generalized hybrid models [3]. Hybrid models start with a functionally complete coarsegrained model which explicitly links DNA replication, metabolism, cell division, and geometry to the external environment [4]. Such models can describe changes in energy and redox equivalents, RNA transcripts, transport, etc. A coarse-grained model must correctly predict robust bacterial reproduction [5], an obligatory test for all whole-cell models. A hybrid model of interest can then be constructed by embedding genome-specific detail for a particular cellular subsystem (e.g. glycolysis, nucleotide biosynthesis, etc), called here a module, into the original coarse-grained model. Recently, 'hybridmodularity' has been successfully demonstrated by developing a detailed nucleotide module [6]. The availability of detailed metabolic reconstructions [1] and genome-scale reaction networks [2] can significantly accelerate the development of genome-specific modules [3] which can then be reused in many large-scale computer hybrid models for a variety of completely annotated genomes.
Another advantage of generalized models is that they combine a detailed summary of the functionality required to sustain the cell's life with modest-size model's complexity. Such models are thus an ideal platform for the development of computationally tractable systems biology concepts and biomathematics approaches. In this paper, we discuss mathematical and computational approaches to evaluate the model's sensitivity and robustness. Specifically, a new strategy for the extension of Metabolic Control Analysis (MCA) [7] to limit cycles is introduced to identify which pseudomolecular processes should be delumped to implement a particular biological function in hybrid cell models. To illustrate the modeling principles and highlight computational issues, the updated Cornell Escherichia coli B/r-A model is used to benchmark the framework. Here v i is the rate of pseudo-reaction i.

The model and computational frameworks
The Cornell E.coli model, depicted in figure 1, represents a single cell of E.coli growing in a glucoseammonium medium. The model describes metabolism, DNA replication, and cell geometry [4]. The modeling principles include: (i) the aggregation of cellular compounds into a manageable number of lumped species, (ii) the use of pseudochemical reactions and accurate stoichiometry, and (iii) the evaluation of as many kinetic parameters as possible from independent measurements. A recently updated model includes 36 ODEs for metabolism and DNA replication, one algebraic approximation of the ribosomal protein biosynthesis, one algebraic equation to monitor the septum growth, and 31 discrete events describing instant changes in the model's parameters and state variables (e.g. changes in gene dosage, cell division, etc.). Dynamic systems describing a smooth evolution coupled with discrete transitions are called hybrid. The E.coli model is thus a hybrid differential-algebraic equation (HDAE), implemented in Matlab ® , C++, and Systems Biology Markup Language (SBML).
The E.coli model event network is depicted in figure 2, where nodes correspond to events and arrows indicate how events can imply one another. We find that DNA replication initiation is the most connected node (i.e. E 3 ) signifying its central role in the cell cycle. Although Matlab ® event detection is used, additional means are needed to identify "secondary" events induced by changes in the HDAE definition at each event. To catch all events, a general event detection algorithm has been developed. A stationary cell division cycle corresponds to a periodic solution of the HDAE. One way to study periodic solutions is to compute a first return or Poincaré map P(s, p), relating any two cell states over least period T, s t+T = P(s t , p). Here vector s includes masses, concentrations and numbers of molecules, and p includes model parameters (i.e. kinetic rates). Each event l is defined by the zero level of a scalar function F l (s, p), F l (s, p) = 0. We assume that s t never corresponds to any event, i.e.   An alternative approach to calculate periodic solutions is to solve a periodic multi-point boundary value problem (BVP) for a discrete closed orbit z = [(t 0 , s 0 , p 0 ), …, (t L+1 , s L+1 , p L+1 )], s 0 = s L+1 [8,9]. Here each (s l , p l ) is chosen just prior to event l. The unknown period T, event times t 1 ,…, t L , and states s l can be found using Newton solvers. Additional mesh points between events and a phase condition are needed to increase the accuracy and uniquely determine the periodic solution, respectively [9].
Typical time courses are shown in figures 3 and 4. We find that while changes in some species (e.g. ammonium ions, proteins) look 'smooth' between divisions, other species (e.g. different forms of DnaA molecules or RNA transcripts) experience complex behavior throughout the entire cell cycle.

Sensitivity and stability of the cell division cycle
Sensitivity analysis is an important tool for model evaluation as well as for quantifying effects of parameter values on model predictions [10]. Specifically, it is important to characterize the relative significance of various intracellular dynamic processes for modeling a growing cell. This can be done by an appropriate extension of MCA to the case of self-oscillations in autonomous systems . Let s(t, p) be a periodic solution with least period T(p), p is a vector of the system's parameters. A well known property of periodic solutions in autonomous systems is that a first-order sensitivity function, is generally unbounded when time tends to infinity [10,11]. Here k = 1,…, K and K is the number of parameters.
The key idea to understand and overcome this analytic difficultly can be seen from the differentiation of periodic condition s(t m (p), p) = s(t 0 , p) with respect to any scalar parameter p k , We find that while metabolic processes contribute to the fastest part of the 'cell clock' with division period T of 45 min, a slower part of the clock has to be transmitted between cell generations.

Conclusion
We are currently developing general MCA and BVP approaches to study the robustness of hybrid whole-cell models when parameters are allowed to vary. This includes relating variations in growth conditions to changes in the number of replicating chromosomes, sensitivity analysis of the DNA replication, identification of independent measurements to fit important parameters, etc. We hope that these approaches will also help us construct large-scale genome specific modules for E.coli and other genomically related gram-negative organisms (e.g. Shewanella oneidensis and Zymomonas mobilis).