Complete in silico cell — integrating genome, transcriptome, proteome, and metabolome to simulate the entire cell cycle from first principles (à la Karr 2012 M. genitalium model)
The Karr et al. (2012) whole-cell model of Mycoplasma genitalium was a landmark achievement in systems biology: the first computational model that predicts cellular behavior from the complete genetic sequence. It integrates 28 sub-models spanning all major cellular processes — metabolism, gene expression, DNA replication, cell division — into a single executable simulation of one complete cell cycle. Understanding its architecture reveals how whole-cell simulation is built.
The WCM is organized as 28 sub-models, each representing one cellular process, integrated through a common state representation and a time-stepping simulation framework:
State variables (shared across all sub-models): • Metabolites: concentrations of 469 metabolites (mM, cytoplasm + membrane) • RNAs: copy numbers of 478 RNA species (mRNA, tRNA, rRNA, ncRNA) • Proteins: copy numbers of 381 protein species (monomers + complexes) • DNA: nucleotide sequence + modification state (methylation) + supercoiling state • Geometry: cell length, volume, surface area • Time: discrete time steps of 1 second
Key sub-models and their mathematical frameworks: 1. Metabolism (FBA): • Flux balance analysis: maximize biomass = min(nutrient uptake, max theoretical yield) • Stoichiometric matrix S (m×n): S·v = 0 (steady-state mass balance) • Bounds on fluxes: transport reactions constrained by transporter copy numbers × kcat
2. Transcription (stochastic): • Promoter states: RNAP bound/unbound (Markov chain) • Initiation: rate = k_init × [RNAP-σ factor] • Elongation: 40 nt/s (M. genitalium RNA polymerase) • Termination: intrinsic (stem-loop) or Rho-dependent
3. Translation (stochastic): • Ribosome binding to RBS (rate ∝ Kozak sequence strength) • Elongation: 5 aa/s per ribosome • Protein folding: immediate (simplified) or chaperone-dependent (HSP70/60 modeled)
4. DNA replication (deterministic ODE + stochastic): • DnaA accumulation → DnaA-ATP complex → oriC loading (threshold model) • Fork progression: 150 nt/s per fork × 2 forks (bidirectional) • Okazaki fragment synthesis: stochastic (primase binding, extension, ligation)
5. Cell division: • FtsZ polymerization kinetics (nucleation-elongation model) • Division checkpoint: chromosome fully replicated AND cell length threshold
Drug perturbation in the WCM: • Each drug maps to its molecular target in specific sub-model • Concentration-effect relationship: Hill equation E = Emax × C^n / (EC50^n + C^n) • Drug updates: RNAP elongation rate reduced → mRNA production cascade effect • Emergent cell-level phenotype: growth arrest, filamentation, lysis • Critical advantage: drug combinations — WCM predicts synergy/antagonism from first principles
Flux Balance Analysis (FBA) is the mathematical backbone of whole-cell metabolic simulation. By treating metabolism as a steady-state linear optimization problem — maximize growth rate subject to stoichiometric constraints and enzyme capacity bounds — FBA can predict which metabolic pathways are active, how fast the cell grows on any carbon source, and how antibiotic interference with specific enzymes ripples through the entire metabolic network.
Flux Balance Analysis formulation:
Optimization problem: maximize c^T · v (maximize biomass objective) subject to: S · v = 0 (stoichiometric mass balance, steady state) lb ≤ v ≤ ub (reaction flux bounds)
Where: • v = flux vector (mmol/gDW/h for each of n reactions) • S = stoichiometric matrix (m metabolites × n reactions) S[i,j] = stoichiometric coefficient of metabolite i in reaction j (negative=consumed, positive=produced) • c = objective coefficients (1 for biomass reaction, 0 for all others) • lb, ub = lower/upper flux bounds (reversibility, max enzyme capacity)
Biomass objective reaction: • Empirically determined: drain of all biosynthetic precursors at ratios matching measured cell composition • For M. genitalium: 0.35 DNA precursors + 0.25 RNA precursors + 0.45 amino acids + 0.08 lipids + ... • Growth rate μ = v_biomass [1/h]; doubling time = ln(2)/μ
Reaction constraints: • Irreversible reactions: lb=0 (transport, spontaneous thermodynamically driven) • Reversible reactions: lb=-1000, ub=+1000 (unconstrained) • Enzyme-capacity constrained: ub = E_i × kcat_i / MW (if proteomics data available) • Exchange reactions: uptake of nutrients constrained by transporter abundances and Km
GEM reconstruction pipeline: 1. Genome annotation: all metabolic genes identified (sequence homology + functional annotation) 2. Reaction assignment: each gene → enzyme → reaction (via KEGG/BiGG database) 3. Stoichiometry: balance each reaction for mass + charge 4. Gap filling: identify missing reactions to connect dead-end metabolites (automated gap-filling tools) 5. Validation: predict growth on minimal media → compare to experimental data
Enzyme-constrained FBA (ecFBA): • Adds protein cost: Σ_i (v_i / kcat_i × MW_i / E_total) ≤ 1 • Limits total enzyme usage to measured proteome mass • More accurate growth rate predictions; reveals enzyme bottlenecks • GECKO (Sánchez et al.): automatic ecFBA model construction from any GEM + proteomics
Dynamic FBA (dFBA): • Time-dependent: external metabolite concentrations change as fluxes consume/produce them • ODE for extracellular [S]: d[S]/dt = -v_uptake × X (where X = biomass density) • Predicts batch culture dynamics: lag phase, exponential growth, diauxic shifts (glucose → acetate) • Used in WCM to couple metabolism to cell growth
DNA replication is arguably the most complex molecular process in the cell, requiring the coordinated action of over 20 proteins to accurately copy 580,000 base pairs in M. genitalium (3 billion in human) with an error rate of 1 in 10^9. Simulating it at the mechanistic level reveals how antibiotics like ciprofloxacin and novobiocin trap DNA topoisomerases and why fluoroquinolone resistance mutations in gyrA are so predictable.
Replication cycle simulation:
1. Initiation (DnaA-driven): • DnaA-ATP accumulates through cell cycle (transcription rate × cell volume growth) • At threshold (Karr model: ~10 DnaA-ATP complexes): loads onto oriC (9 DnaA boxes) • DnaA-ATP cooperativity: Hill coefficient n=3 → sharp threshold switch • DnaBCD helicase loading → replication bubble opens → primase (DnaG) synthesizes RNA primer • Bidirectional: two forks depart oriC → terminate at ter sites (180° from oriC)
2. Fork progression simulation (stochastic): • Each fork tracked as position on chromosome (integer, 0–580,000) • Advance rate: 150 nt/s × (dNTP availability) × (1 - gyrase inhibition) • dNTP pools: coupled to metabolism (FBA provides dNTP synthesis fluxes) • Leading strand: continuous synthesis • Lagging strand: Okazaki fragments (stochastic primase binding every ~1–2kb) • Each Okazaki: primase binds (k_on = 0.1/s), extends 1–2kb, DNA Pol I removes primer, ligase seals
3. Topoisomerase function and drug mechanism: • DNA supercoiling: replication fork introduces positive supercoils ahead • DNA gyrase (Topo II): introduces negative supercoils (2 ATP per decatenation) • Without gyrase: Lk positive → fork stalls (DNA over-twisted) • Ciprofloxacin: binds gyrase-DNA cleavage complex → prevents re-ligation → DSB • Simulation: gyrase inhibited → fork rate drops proportionally → cell cycle arrest • Resistance: GyrA A83T (E. coli numbering) → reduced ciprofloxacin binding (Kd 100× higher)
4. SOS response (drug-induced): • DSBs activate RecA → filaments on ssDNA • RecA-ssDNA → LexA autocleavage → LexA repressor removed → SOS genes induced • SOS regulon: 40+ genes including umuC/umuD (error-prone Pol V), recA/recN • Pol V: copies past lesions (translesion synthesis) → increased mutation rate → resistance evolution • WCM tracks: each DSB event → RecA loading → LexA repressor level → SOS gene activation timecourse
Gene expression is the ultimate link between genotype and phenotype: the same genome can produce radically different protein compositions (and therefore behaviors) depending on which genes are transcribed, how efficiently mRNAs are translated, and how quickly proteins are degraded. Whole-cell simulation resolves stochastic gene expression at the single-molecule level — capturing the burst dynamics, translational noise, and feedback loops that determine cell-to-cell variability.
Gene expression simulation in the WCM uses a multi-scale stochastic framework:
1. Transcription model: Promoter occupancy (Markov chain): • State: free / RNAP-bound / elongating • k_on: σ-factor × RNAP concentration × promoter recognition rate • k_off: transcript abortion rate • k_elg: elongation initiation rate (rate-limiting for many stress-regulated genes) Elongation tracking: • RNAP position tracked (nt resolution) — steric exclusion between consecutive RNAPs • Pause sites: specific sequences → RNAP pause (power-law dwell time distribution) • Termination at terminator hairpin or Rho-dependent termination
Burst kinetics: • Prokaryotic transcription: genes fire in "bursts" (multiple mRNAs during one RNAP binding event) • Burst frequency: determined by promoter on-rate (k_on) • Burst size: geometric distribution with mean = k_elongation / (k_abort + k_terminate) • Protein variation: Fano factor (Var/mean) increases with burst size
2. Translation model: Polysome dynamics: • Multiple ribosomes on single mRNA (polysome) • Ribosome initiation rate: k_init × [30S complex] × SD sequence strength • Codon usage: rare codons → slow elongation → ribosome queue • Protein folding: cotranslational folding modeled by probability of aggregation Protein degradation: • ATP-dependent: Lon, ClpXP, ClpAP proteases • Rate: first-order ~ k_deg × [P] (half-life: stable proteins 1/μ = dilution; unstable: k_deg >>μ) • Ssrpeptides: tmRNA tagging of stalled ribosomes → ClpXP-directed degradation
3. Cell-to-cell variability: • Stochastic simulation (Gillespie algorithm) of each mRNA and protein • Division partitioning noise: each molecule randomly partitioned (binomial) to daughter cells • Result: even isogenic population shows 10–30% protein-level CV (coefficient of variation) • Drug tolerance: cells with low gyrase copies tolerate ciprofloxacin → "persister" cells
4. Protein complex assembly: • Ribosome (21 proteins + rRNA): co-transcriptional assembly modeled • RNA polymerase holoenzyme: core α₂ββ'ω + σ factor (stochastic σ factor competition) • Replisome: 12 subunits assembled at oriC → tracked as single complex
Cell division is the endpoint of the cell cycle and the process that determines bacterial fitness and antibiotic susceptibility. The FtsZ protein polymerizes into a Z-ring at the cell midpoint and drives membrane constriction — a process that has been the target of intense antibiotic development. Whole-cell simulation of FtsZ dynamics can predict the effect of FtsZ-targeting antibiotics and reveals why FtsZ inhibitors are promising against drug-resistant bacteria that have exhausted conventional antibiotic targets.
FtsZ division simulation model:
1. FtsZ polymerization kinetics: • GTP-FtsZ monomer → nucleation (rate: k_n × [FtsZ]^2) → elongation • Elongation: k_on × [FtsZ-GTP] (association); k_off (dissociation) • GTP hydrolysis: k_hyd → GDP-FtsZ → depolymerization • Net: treadmilling behavior (GTP end grows, GDP end shrinks) • Midcell localization: MinCDE oscillation system defines Z-ring position MinC/MinD oscillate pole-to-pole every ~40s → time-averaged inhibition highest at poles → FtsZ excluded from poles
2. Z-ring assembly and maturation: • Critical concentration for polymerization: ~1–2 μM FtsZ • FtsZ polymers condense to midcell ring structure • Divisome assembly: FtsZ → FtsA/ZipA (anchors) → FtsK → FtsQ/L/B/W/I/N (inner membrane proteins) • Constriction: FtsZ GTP hydrolysis → bending force on membrane • FtsZ bending model: protofilaments preferentially adopt curved conformation → inward pulling force
3. Division checkpoint: • DNA segregation check: FtsK helicase translocates DNA away from division plane • SOS check: sulA (recA-induced) inhibits FtsZ polymerization during DNA damage → delay division • Chromosome completion: last base pair replicated → FtsK releases→ division proceeds • Cell length threshold: sensor feeds back to FtsZ nucleation rate
4. Drug mechanism simulation (FtsZ inhibitors): • PC190723: binds cleft between C-terminal H7 α-helix and T7 loop of FtsZ • Mechanism: stabilizes polymers → reduces GTPase → no treadmilling → static ring → no constriction • WCM simulation: FtsZ treadmilling rate × (1-E_drug) → constriction velocity reduced → cell filaments • Berberine: inhibits FtsZ polymerization (competition with GTP) • Synergy with ciprofloxacin: ciprofloxacin activates SOS → FtsZ inhibited by SulA AND PC190723 → strong synergy • Resistance: FtsZ A316P mutation (conformational change) → PC190723 Kd 100× higher
5. Whole-cell model validation: • Karr 2012 WCM validation against 110 experimental observations: - Cell cycle duration within 2-fold of measured - Gene essentiality predictions: 92% accuracy vs. transposon screen - Growth rate on different media: within 30% - RNA/protein mass ratios: within 20% • Drug response prediction accuracy: 78% of drug MIC predictions within 4-fold of measured • Modern extensions: VIVARIUM (open-source WCM framework), now applied to E.coli (4,285 genes)
The most striking validation of the Karr whole-cell model was its prediction of 19 previously-unknown essential genes in M. genitalium — genes that FBA alone would classify as dispensable because they have no metabolic role, but which are essential through indirect regulatory connections only visible in an integrated whole-cell simulation. This demonstrates the unique power of WCM over any single-omics analysis: emergent essentiality from network connectivity.