The general workflow of constructing an initial cell state is shown in Figure 1. Starting from the coordinates of the 503 ribosomes and cell boundary of a small cell with

$\sim$ 400-nm diameter from cryo-ET (Figure 1A), the DNA is folded around ribosomes as a circular self-avoiding polymer on a lattice in such a way that the sequence order and gene positions are maintained (Figures 1B and 1C) (

). Experimental 3C maps showed no significant features of persistent supercoiled domains or loops, so the chromosome is assumed to be in a relaxed state (

). Each replicate cell uses the same ribosome coordinates from cryo-ET, but a different chromosome configuration unique among the ensemble. The top of the membrane is cut away in Figure 1C to reveal the ribosomes and DNA. According to the 3C-seq maps and the proteomics of nucleoid-associated proteins (NAPs) (

;

), the DNA configurations are assumed to be relaxed with no supercoiling so that the genes are easily accessible. Figure 1D shows 120 degradosome complexes in red, 66 SecY proteins in blue, and 831 PtsG proteins in green as three examples. The remainder of the proteome consisting of over 77,000 proteins, 200 mRNA, and 5,800 tRNA are then randomly distributed throughout the cytoplasm and membrane, resulting in the crowded environment (Figure 1E).

Our spatial model includes a total of 7,765 unique molecules and intermediates and over 7,200 reactions including binding reactions such as RNAP binding to a gene start site. Where possible, kinetic parameters are obtained from single-molecule experiments, such as the smFRET experiment for the formation of the DnaA filament along the AT-rich single-stranded DNA (Figure 1F) near the origin. Otherwise, as described in STAR Methods (Metabolic rates and parameterization), kinetic parameters are developed from a targeted survey of the primary literature or kinetic databases (Figure 1G) as discussed above. Simulations of a spatially resolved cell are computationally expensive and require GPU (Figure 1I) acceleration to make them possible on a human timescale (

). The GPUs used for spatial simulations included NVIDIA Titan V and NVIDIA Tesla Volta V100 GPUs, which took 10 h and 8 h to simulate 20 min of cell time, respectively. Because of this computational expense, we simulated the first 20 min of the cell cycle for only 8 cells, a limited time frame during which we assume no substantial growth or DNA replication has occurred. The simulations provide insight into the numbers of active degradosomes, RNAP, and ribosomes (Figures 1J–1L). The early increases are due to the initial conditions of the spatial simulations. The cell is initialized with no active complexes (no RNAP are on the chromosome and no mRNA are on ribosomes and degradosomes), and the transient behavior reflects the time required for RNAP to diffuse to genes and mRNA to be translated or degraded.

calculated that anywhere from 15.5% to 36.2% of RNAP are active at any one time in E. coli depending on the doubling time, with slower-growing cells having a smaller fraction of active RNAP. The spatial model predicts that Syn3A will have an average of 63 of its 187 RNAP active (34%) early in the cell cycle, falling within the calculated range. For fast-growing E. coli, approximately 80% of the ribosomes are active (

;

), but for slow-growing E. coli, this number drops between 20%–50% depending on the growth rate (

). The spatial model predicts that, on average, 220 of the 503 ribosomes, roughly 45% of ribosomes, are active.

We calculated the mRNA half-lives, the number of times a gene is transcribed, and the number of times an mRNA is translated in its lifetime for all the 452 protein-coding genes (Figures 1M–1O). The average and median half-lives are in reasonable agreement with the 2-min average half-life experimentally measured in Mycoplasma gallisepticum (

). The broad distribution of half-lives, including the long tail out to 15 min, has been observed in a genome-wide study of mRNA half-lives in B. subtilis,

. Each gene is transcribed at least once within the ensemble of simulations, but not in every simulation. The number of times each gene is transcribed reflects both its length and, more importantly, its promoter strength, which is weighted relative to proteomics counts. Lastly, the genome-wide average translations per mRNA is four times, but several factors impact this number including gene length and how many times a mRNA can be read by a polysome using a polysome spacing of 120 nt estimated from a distribution of polysome sizes in E. coli (

).

The only sugar source in the defined medium is glucose, so the revised map for central metabolism (Figure 2) starts with the phosphorelay relay system (

;

,

), which is responsible for the uptake and phosphorylation of glucose to glucose-6-phosphate (g6p) and is shown in the inset. Each phosphate exchange reaction of the phosphorelay is simulated independently, and the overall kinetics for the phosphorelay predict that Syn3A takes up 15,000 glucose molecules per second for a cell with a radius of 200 nm. Syn3A does not have the proteins to perform oxidative phosphorylation, so all ATP in Syn3A is generated by the central metabolism (

). Pyruvate kinase (pyk/0221) (PYK) converts ADP to ATP and pep to pyruvate, competing for pep molecules with the glucose transport reaction, so the fluxes between these two reactions need to be carefully balanced throughout the cell cycle. Because the fructose-1,6-bisphosphate aldolase (fbaA/0131) (FBA) reaction splits fdp into two molecules, the rate of lower glycolysis is twice that of upper glycolysis. Therefore, a maximum of 45,000 ATP can be generated per second assuming no other NTPs are being formed: 30,000 ATP per second can be generated by the phosphoglycerate kinase (pgk/0606) (PGK) and 15,000 by pyruvate kinase (PYK) where pep is split between the glucose uptake reactions and PYKs. A much smaller amount of ATP is generated through acetate kinase (ackA/0230) (ACKr) upon the secretion of acetate. In our model, the cell was able to survive off of the ATP generated by central metabolism. Detailed discussions of the amino acid, cofactor, nucleotide, and lipid metabolic subsystems are included in STAR Methods.

The kinetic model of the genetic information processing reactions including DNA replication initiation and elongation, mRNA degradation, transcription, and translation are based on the kinetic model by

with a few modifications. Each of the elongation reactions are treated using a polymerization rate dependent on the respective monomer concentrations (dNTPs, NTPs, aa:tRNAs). For full details on these rate forms, see STAR Methods (Genetic information processes).

Unlike the spatial RDME-CME-ODE model, cells are simulated for whole cell cycles in the well-stirred CME-ODE model. To highlight the key features defining the cell cycle in Syn3A, we examine in Figure 3 the time dependence for the initiation of DNA replication, chromosome duplication, and the doubling of both the cell volume and surface area in the well-stirred simulations. In our previous work (

), only one replication event was allowed to occur per cell cycle, but here we complexify this model and allow for multiple replication initiation events based solely on the kinetics. We use the same equation for the rate of DnaA(III) filament formation on ssDNA for each independent origin

Both

$k_{on}$ and

$k_{off}$ for the binding of DnaA(III) were measured using smFRET (

). We allow replication initiation events on the daughter chromosomes in our model after they separate following the first replication cycle. The initial cell volume is used for the first replication initiation event, and the doubled cell volume is used for the daughter chromosomes late in the cell cycle. Based on these kinetics, if the number of DnaA doubles faster than the volume, it is likely that another replication initiation event will occur and the cell will have more than two chromosomes at the end of its cell cycle.

Although other bacteria have been shown to undergo multiple replication events per cell cycle (

;

), it had not yet been directly observed in Syn3A. From quantitative polymerase chain reaction (qPCR) experiments, we have determined the relative quantities of origins, quarter positions, and termini in Syn3A cells in both exponential and stationary phases (Figure 3D). The results are presented with all the quantities scaled to the number of termini. In the exponential phase cells, there are more than three times the number of origins than termini on average. This indicates that in many cells, after the first replication initiation event occurs, another replication initiation event will occur on the same chromosome before the first replication cycle completes. In our current model, a maximum ratio of 2 would occur, as we do not allow multiple initiation events to occur until the first replication cycle completes.

Critical to determining the length of the cell cycle are the times required to double the volume and surface area of the cell. The kinetics for cell growth are characterized by the formation of lipids and insertion of membrane proteins to determine the cell surface area and volume in the model. The cell volume is calculated from the surface area assuming the cell maintains spherical morphology until the onset of division. The volume doubles in the simulations anywhere from 56 to 72 min with an average of 64 min (Figure 3E). Without explicit kinetics for cell division by FtsZ/FtsA filaments, division is assumed to begin when the volume has doubled, during which the volume stays constant and the surface area continues to grow until two separate cell are formed. The surface area will double anywhere from 88 to 112 min with an average doubling time of 97 min (Figure 3F). Syn3A has an experimentally measured doubling time of 105 min in rich growth medium (

). We report a simulated doubling time (Figures 3E and 3F) based only on healthy cells, whereas the experimentally measured doubling time includes a whole population, which may include unhealthy cells, reducing the average doubling time of the colony.

To determine the cell surface area, each lipid and membrane protein has an assigned surface area contribution (STAR Methods, Lipid metabolism). The contributions of proteins and lipids to the surface area are separated in Figure 3F where the two contributions are seen to maintain a rough 55% to 45% surface area contribution ratio (

). Because the translation and insertion of membrane proteins are both totally stochastic in the simulations, there is more variation in the surface area contribution from membrane proteins. The variation in lipids increases with simulation time as lipid synthesis genes are stochastically expressed in each individual cell.

The liponucleotide synthase cdsA/0304 catalyzing reaction DASYN (Figure S2) that adds CDP to a phospholipid precursor (phosphatidic acid, PA) was identified as a “choke point” in the phoshpolipid biosynthesis, in agreement with a recent whole-cell model of E. coli (

). Additionally, acyltransferase plsY/0117, which catalyzes the conjugation of fatty acid and glycerol moieties at the membrane was found to limit the production of the downstream intermediate PA. We attribute both of these effects to their low counts in the reported proteomics of Syn3A (

) and adjusted the counts to values similar to those observed in other bacterial species (Table S2 and STAR Methods, Lipid metabolism). Their low counts are likely due to the fact that both are multiple domain membrane proteins, which are known to be underreported by proteomics when only a trypsin digest is used (D. Gonzalez, personal communication).

Due to reduction in its minimal genome, Syn3A has few remaining transcription/translation/transport regulatory proteins and must adjust the fluxes through the cellular subsystems to maintain stable growth. The simplified map of the reaction network with fluxes from a representative cell from the well-stirred model early in its cell cycle demonstrates the balance among use of ATP and GTP, nucleotide metabolism, and glycolysis (Figure 4A). The glycolysis pathway and nucleotide metabolism are connected through the PYK reactions converting all (d)NDPs to (d)NTPs, which results in the shared usage of pep with glucose uptake. Syn3A exclusively makes pep by the action of enolase at the end of glyolysis, generating two pep per glucose taken up. Because the glucokinase was removed in genome reduction (

), we assume that the only way Syn3A can phosphorylate glucose is by PtsG in the phosphorelay. According to our kinetic model, if the cell runs out of pep, there is no way to continue glycolysis or phosphorylate NDPs except ADP, which can be converted to ATP by the reversible ATP synthase. On average, if PYK reactions use more than half of the pep formed, less glucose will gradually be taken up and the cell will run out of pep and cease to take up glucose. Roughly 16% of the cells in a total ensemble of 207 cells in the well-stirred simulations experienced pep’s shortage, leaving 174 cells that could successfully complete a cell cycle.

Previous studies have indicated the possibility of PYK, PFK, and GAPD being rate-limiting reactions of glycolysis (

;

); however, fructose-1,6-bisphosphate aldolase (FBA) appears to control the overall flux of glycolysis in Syn3A according to our simulations, which agrees with the findings of

). FBA has a low experimental proteomics count of 227 relative to the other glycolytic enzymes in Syn3A, having counts of 400 or greater. Relative to other bacteria, Syn3A appears to have a lower FBA count (Table S1), so in parameterizing the simulations, the count of FBA was scaled to 775 based on the enzyme’s concentration in E. coli so that kinetic parameters would better match the known equilibrium constant for the reaction. The FBA reaction is almost always at its maximum rate and is the slowest step in upper glycolysis in our models. The rates of the reactions in lower glycolysis are dictated by how much flux goes through the FBA reaction.

The balance of metabolic reactions also affects the rates of genetic information reactions through the monomer-dependent polymerization reaction rates discussed in STAR Methods (Genetic information processing). The time-dependent rates of each elongation reaction (DNA replication, transcription, and translation) scaled to its maximum rate for dnaA/0001 expression are shown in Figure 4B. The single cell in orange has only one complete DNA replication event occur in its cell cycle, and the green cell has a second. Their replication rates are at the maximum until the cell runs low on dNTPs (

$<$ 0.01 mM) (Figure S4), in this case dATP (data not shown). The slowdown gives the cell time to import more deoxynucleosides and generate more dNTPs. The rate of replication will fluctuate from minute to minute as long as the cell runs low on a particular dNTP. In general, DNA replication rate is the most frequently affected of the genetic information processes with its average (black) going as low as 75% of the maximum.

Translation is infrequently altered because the amount of charged tRNA (aa:tRNA) depends directly on the uptake of amino acids, which are present in millimolar concentrations in the medium. So even though cells will sometimes run low on amino acids, a brief slowdown is enough for them to import more amino acids and recover their charged tRNA levels (

$>{10}^2$ of each aa:tRNA).

A key feature of our model is the explicit tracking of every energy molecule used by activated reactions in both metabolism and genetic information processes. Already in 1973,

comprehensively calculated the amount of ATP required for the formation of a microbial cell based on the studies available on energetic costs at the time. He broke down the ATP demands of a cell into requirements for formation of polysaccharides, proteins, lipids, RNA, and DNA, uptake of amino acids, phosphate, and ions, and turnover of RNA. He notes that exactly accounting for transport reactions is difficult and little information on their ATP costs were available at the time, so only a few transport reactions are included in his calculations. More recently, a comprehensive review of the literature and calculations for both bacteria and eukaryotes (

assigned costs for the synthesis of macromolecules and determined the total ATP costs for DNA replication, transcription, and translation.

While the previous reviews calculated ATP costs for overall cell formation, we advance ATP cost calculations to single-cell, single-reaction resolution as a function of time for the cellular networks of Syn3A including both metabolism and genetic information processing, which is made possible by fully dynamical kinetic modeling. Figure 5A shows the ATP generated and used as functions of time for a representative cell from the well-stirred simulations. Note that these plots do not represent the number of phosphate bonds made and broken, but the number of ATP molecules being used. Another notable difference is that tRNA charging is only counted as one ATP, whereas it is typically counted as two phosphate bonds (

). The charging reactions still convert ATP to AMP and pyrophosphate in the simulations (Figure S3), but because the metabolic network explicitly accounts for the ATP cost of regenerating an AMP to an ADP through the ADK1 reaction (Figure S1). tRNA charging is counted as a single ATP cost in Figure 5. Additionally, because translation elongation uses GTP instead of ATP, it is not shown in the cost plot, but it is twice what we defined as the ATP cost of charging tRNAs from two GTP per amino acid: one during the loading of an aa:tRNA into the A site of the ribosome and a second during translocation to the next step on the mRNA (

).

Quantifying the exact cost of activated transport reactions has been a difficult challenge in both the energy calculations by

and

, as well as in other recent whole-cell models (

;

). Because each transport reaction is simulated independently, we know their exact ATP costs by recording the fluxes through each reaction. Unlike most organisms, which have synthesis pathways for most of its building blocks, Syn3A has been reduced to the point where it relies on having to transport them in. From crystal structures of related transporters, it is clear that the majority of them contain ATPase domains that require, on average, at least one ATP for every nutrient molecule imported (

). It is assumed here that one ATP is used for every molecule taken up. The cost of active transport is approximately twice the ATP cost of tRNA charging (roughly 5,000:2,500 ATP per second), making it one of the largest energetic costs in the minimal cell. This is the most exact calculation of the costs of transport presented to date and reveals how transport reactions are critical for a minimal cell, which relies almost entirely on nutrients from its environment.

In Mycoplasmas, the ATP synthase typically breaks down ATP and pumps out protons to maintain a basic intracellular environment (

,

). The reversible kinetics for ATP synthase depend on the ATP, ADP, and phosphate levels inside the cell. While ATP synthase accounts for roughly 10% of the overall ATP costs for the majority of the time, in response to low ATP levels in the cell, it can momentarily change direction.

Homeostasis is a property of a normal cell to maintain constant intracellular concentrations over a cell cycle suggesting that upon experiencing a perturbation whether from the environment or an intracellular reaction, it responds by adjusting its biochemical pathways to bring the concentrations back into an acceptable range for stable functioning of its networks (

). The prevailing wisdom is that some degree of regulation is required to control any large fluctuations. The time traces for a representative selection of metabolites and macromolecules from the well-stirred simulations are shown in Figure 6. The complete time traces of all chemical species and reaction fluxes are provided in Data S1. The dNTP concentrations decrease late in the cell cycle for cells where multiple replication events have occurred because they are being incorporated into new chromosomes. For dATP and dTTP, it appears that the same cells maintain average concentrations that are relatively constant over the whole cell cycle. The concentrations of dGTP and dCTP both continually increase over the cell cycle, which calls for investigating the possibility of any regulation on the uptake of their precursor deoxynucleosides or the thioredoxin reactions converting GDP and CDP to dGDP and dCDP in nucleotide metabolism. The concentrations of UTP and CTP are fairly constant throughout most of the cell cycle, likely because the only annotated way to make new CTP in Syn3A is by converting UTP to CTP through a CTP synthase reaction CTPS2 in nucleotide metabolism (Figure S1). ATP and GTP, on the other hand, continue to increase over the cell cycle in both cells with single and multiple replication events.

Even though homeostasis can be observed for a population, there can be significant variation among individual cells due to stochastic fluctuations in gene expression. Some of the largest variations in our simulations were for phosphate (PI), pyrophosphate (PPI), and fructose bisphosphate (fdp). Cells with multiple replication events have a wide range of phosphate levels at the end of the cell cycle. In contrast, cells with only a single replication event return to consistently lower phosphate levels after the first replication event is complete. PPI is given off from DNA replication reactions, so cells with multiple replication events will see a broader range of PPI and therefore PI levels. While such high concentrations have been reported in yeast (

), they may be inaccurate for Syn3A. Syn3A still has the gene phoU/0428 coding for a phosphate regulator, so it is a primary candidate to be included in a complexified model that includes regulation.

To gauge proteome-wide homeostasis, the scaled protein counts for all proteins reported in the proteomics data (

) with counts of 10 or greater excluding ribosomal proteins, a total of 350/452 proteins, at the end of a cell cycle are shown in the histogram in Figure 6. Ribosomal proteins are excluded from this plot because their counts do not reflect the 503 ribosomes observed in cryo-ET, with many having counts fewer than 300. For a protein to maintain a near-constant concentration, its count must double over a cell cycle as the volume doubles. On average, the overwhelming majority of proteins end the cell cycle with 1.75 to 2.25 times their initial protein counts. Outliers include proteins whose genes are longer than 4,000 nt on the low end and priB/0026 on the high end, which has a transcript only 441 nt long.

To gauge the quality of our parameterized well-stirred hybrid CME-ODE model to reproduce the results of the 3D spatially resolved hybrid RDME-gCME-ODE model with diffusion (see Video S1), we compare counts of mRNAs for genes involved in genetic information processes, nucletoide pools, protein distributions, and mRNA half-life distributions (Figure 7). The mRNA counts in the spatial model are higher on average likely due to the difference in transcription rates between the two models (see STAR Methods, Transcription). The spatial model has lower concentrations of nucleotides (Figures 7C and 7D) than the well-stirred model, which can be tied to two factors: first, more nucleotides are being incorporated into mRNA, and second, the current spatial model does not include DNA replication, which initiates, on average, around 10 min in the well-stirred model (Figure 3B). With the genes for the nucleoside transporters being close to the origin, their genes would be duplicated early in the cell cycle. Consequently, the spatial model should have fewer nucleoside transporters, which would result in slightly lower uptake rates of nucleosides. Scaled protein counts are compared in Figures 7E and 7F. The slightly higher counts produced in the well-stirred simulations reflect that partial replication events have taken place in a few of the cells within the first 20 min. Finally, mRNA half-lives are compared between the two methods (Figures 7G and 7H), where the well-stirred half-lives depend on the active degradosome statistics in the spatial model (STAR Methods, mRNA Degradation). The well-stirred model has longer half-lives on average than the half-lives in the spatial model, with average half-lives of 3.4 min and 1.97 min, respectively.

The RDME-CME-ODE simulations are currently limited by having the degradosome, RNAP, and ribosome complexes all initialized in inactive states. Starting from an inactive state results in the initial transience in Figures 1M–1O, where the first few minutes of simulation are dominated by the complexes reaching steady-state processing of mRNAs, which may be overshadowing other interesting phenomena. In the absence of experimental knowledge of the average active complexes in Syn3A, the initial transience emphasizes importance of diffusion and how the ensemble-averaged results of the spatial model could be used to parameterize probabilistic factors in the well-stirred CME-ODE kinetic model, which account for diffusion. Future simulations with averaged occupation states will address this limitation. Besides the lack of explicit regulatory factors discussed above, the reactions to modify nucleobases in DNA and rRNA have been neglected and the time-dependent assembly mechanisms of protein complexes and ribosomes have been absorbed into the overall rates. Regulation is important and will be included in future models. No kinetic model for formation of FtsZ/FtsA filament and septum formation prior to formation of the daughter cells is considered. In the future, these processes will be addressed as we modify the spatially resolved kinetic model to allow a change in cell morphology, DNA replication, and ribosome biogenesis similar to the rule-based model we created for DNA replication in a slow-growing and dividing E. coli cell (

).