A Quantitative Study of the Division Cycle of Caulobacter crescentus Stalked Cells

Progression of a cell through the division cycle is tightly controlled at different steps to ensure the integrity of genome replication and partitioning to daughter cells. From published experimental evidence, we propose a molecular mechanism for control of the cell division cycle in Caulobacter crescentus. The mechanism, which is based on the synthesis and degradation of three “master regulator” proteins (CtrA, GcrA, and DnaA), is converted into a quantitative model, in order to study the temporal dynamics of these and other cell cycle proteins. The model accounts for important details of the physiology, biochemistry, and genetics of cell cycle control in stalked C. crescentus cell. It reproduces protein time courses in wild-type cells, mimics correctly the phenotypes of many mutant strains, and predicts the phenotypes of currently uncharacterized mutants. Since many of the proteins involved in regulating the cell cycle of C. crescentus are conserved among many genera of α-proteobacteria, the proposed mechanism may be applicable to other species of importance in agriculture and medicine.


Introduction
The events of the cell division cycle must be carried out in a coordinated fashion. Coordination is maintained by underlying molecular regulatory systems of great complexity. Intensive studies of these protein interaction networks by mathematical modeling have assisted our understanding of cell cycle regulation in yeasts [1][2][3][4][5], frog eggs [6][7][8], and even mammalian cells [9,10]. In addition to reproducing large amounts of experimental data, these models have made successful predictions and guided further experimental studies [11][12][13][14].
Although progress in understanding cell cycle regulation in bacteria has lagged behind eukaryotes, the recent discovery of master regulatory proteins, CtrA and GcrA, in Caulobacter crescentus [15,16] allowed us to propose a closed loop of signaling events controlling the cell cycle in this bacterium [17]. Central to this proposal is a CtrA-based bistable switch [17] that can be driven to the ON state by GcrA and to the OFF state by cell division. This simple model, though providing some insight into the logic of cell cycle regulation in Caulobacter, neglected many important aspects of the control system. In this paper, we add more genetic and molecular details to the Brazhnik-Tyson model in order to build a computational model of sufficient accuracy to account for the phenotypes of many mutant strains, both well-characterized and yet-to-be-studied strains. We have incorporated a third important regulatory protein, DnaA [18,19], and the effects of DNA methylation on gene expression. In particular, as the DNA replication fork progresses along the bacterial chromosome, it may turn on the expression of an inactive (fully methylated) gene by creating a pair of hemimethylated genes (old strand methylated, new strand unmethylated). These genes are returned to their fully methylated state by CcrM, a methyltransferase whose synthesis is induced by CtrA late in the cycle. (Some genes are active in the methylated state and inactive in the hemimethylated state.) Genes that are regulated in this fashion can be turned on/off in a strict temporal order, according to their location on the chromosome. The role of DNA methylation in governing progression through the Caulobacter cell cycle has recently been emphasized by Collier et al. [20].
At this stage of the model, the regulation of CtrA proteolysis has been incorporated in a simplistic way, concentrating on the phosphorylation of DivK in the stalked cell compartment at cell division and ignoring (for now) the roles of other proteins, such as RcdA, CpdR, and ClpXP [21]. In addition, we do not attempt to model the control of CtrA activity by phosphorylation [22,23], nor do we take into account explicitly the spatial localization of proteins [24,25]. These aspects of the control system are reserved for a later version of the model.
A Consensus Picture of Cell Cycle Controls in C. crescentus C. crescentus is a dimorphic bacterium that inhabits freshwater, seawater, and soils, where it plays an important role in global carbon cycling by mineralizing dissolved organic materials [26]. C. crescentus normally undergoes an asymmetric cell division resulting in two different progeny cells ( Figure  1): a motile, flagellated swarmer cell and a sessile stalked cell [22,23,27]. The nascent stalked cell then enters immediately into a new round of cell division and produces, about 90-120 min later, a new swarmer cell. The nascent swarmer cell swims

Author Summary
The cell cycle is the sequence of events by which a growing cell replicates all its components and divides them more or less evenly between two daughter cells. The timing and spatial organization of these events are controlled by gene-protein interaction networks of great complexity. A challenge for computational biology is to build realistic, accurate, predictive mathematical models of these control systems in a variety of organisms, both eukaryotes and prokaryotes. To this end, we present a model of a portion of the molecular network controlling DNA synthesis, cell cycle-related gene expression, DNA methylation, and cell division in stalked cells of the aproteobacterium Caulobacter crescentus. The model is formulated in terms of nonlinear ordinary differential equations for the major cell cycle regulatory proteins in Caulobacter: CtrA, GcrA, DnaA, CcrM, and DivK. Kinetic rate constants are estimated, and the model is tested against available experimental observations on wild-type and mutant cells. The model is viewed as a starting point for more comprehensive models of the future that will account, in addition, for the spatial asymmetry of Caulobacter reproduction (swarmer cells as well as stalked cells), the correlation of cell growth and division, and cell cycle checkpoints.
reactions are regulated to control the fraction of active CtrA is poorly understood.
GcrA is an activator of components of the replisome and of the segregation machinery [15], and also regulates genes such as ctrA, pleC, and podJ [15,19]. GcrA protein concentration varies through the cell division cycle, peaking early in the cycle in stalked cells and reaching its minimum in a swarmer cell, after cell division. The DNA replication-initiating protein, DnaA, is required for gcrA expression [18]. In addition, transcription of gcrA is repressed by the CtrA protein [15].

DNA Replication
DNA replication proceeds in three phases: initiation, elongation, and termination. The origin of DNA replication (C ori ) in C. crescentus has one potential binding site for DnaA, a protein involved in initiating DNA synthesis [51]. The DnaA binding site partially overlaps with five CtrA binding sites in C ori [33][34][35]. CtrA represses initiation of DNA replication [30]. Thus, DNA replication is only initiated when DnaA level is high and CtrA level is low. In addition, DNA replication cannot be re-initiated until the origin stie has been fully methylated [52,53]. These conditions prevail during the swarmer-to-stalked cell transition, and just after division in the stalked cell compartment [34]. Once initiated, DNA synthesis continues bidirectionally along the circular chromosome, with an average speed of ;20.5 kb/min in minimal broth, finishing in the late predivisional cell [54]. Elongation of newly replicating DNA strands requires a complex machinery, many components of which are under GcrA control [15].

DNA Methylation by CcrM
Several cell cycle-related genes (ctrA, gcrA, dnaA, ftsZ, and ccrM) have GANTC methylation sites in their promoters [19,31,40,52,53,55,56]. Hence, the expression of these genes may be sensitive to the methylation state of the promoter. DNA replication transforms a fully methylated gene (both strands methylated) into a pair of hemimethylated genes (only one strand methylated). At some later time, the unmethylated strands become methylated by the action of CcrM to return  [77]) Regulation of genes by CtrA is shown in red, by GcrA in blue, by DnaA in green, and by CcrM in cyan. The cyan stars indicate those genes whose transcription is regulated by DNA methylation. The CtrA-driven up-regulation of the dnaA gene (red line with ? mark, left) is suggested by microarray data [32]. DnaA self-regulation (blue line with ? mark, left) is proposed from the fact that the dnaA promoter has DnaA boxes [55]. doi:10.1371/journal.pcbi.0040009.g002 the genes to the fully methylated state [53]. These methylation transitions may affect the expression of cell cycle-related genes [53]. Methylation of C ori is also necessary for initiating a new round of DNA synthesis [34]. These methylation effects provide feedback from the progression of DNA replication to the cell cycle control system.
In C. crescentus and other a-proteobacteria, CcrM is the methyltransferase that accounts for methylation of newly synthesized DNA strands. ccrM transcription is activated by CtrA only from a hemimethylated chromosome for about 20 min, in a late predivisional cell (its expression peaks at ;105 min in the 150-min swarmer cell cycle) [57]. Lon protease is required for CcrM degradation [58]. The half-life of CcrM is less than 10 min in vivo [39]. Thus, CcrM activity is limited to a short portion of the predivisional cell phase, just before cell division.

The Septal or Z-Ring
The multicomponent Z-ring organelle, which forms and constricts at the mid-cell plane, plays an important role in compartmentation of the predivisional cell and its subsequent division [27]. Compartmentation lasts about 20 min [59]. After the late predivisional cell is divided into two progeny cells, the Z-ring is disassembled and degraded.
The Fts proteins (FtsZ, FtsQ, FtsA, and FtsW) have been identified as crucial elements of the Z-ring. ftsZ expression is positively and negatively regulated by CtrA [29,60], and it may also by regulated by DNA methylation since the ftsZ promoter has a methylation site [40,53]. The ftsQ gene is expressed only after CtrA-mediated activation in the late predivisional cell [41]. The FtsQ protein localizes predominantly to the mid-cell plane of the predivisional cell, consistently with the appearance of the Z-ring [61,62]. The FtsA protein exhibits the time course similar to FtsQ [61].
Polar Distribution of DivK and DivK;P divK transcription is activated by CtrA in late predivisional cells, which results in a slight elevation of DivK protein, otherwise present throughout the cell cycle at a nearly constant level [42,63]. The total amount of DivK;P, the form that promotes CtrA degradation, does not appear to undergo dramatic changes during the cell cycle. It is 50% 6 20% lower in swarmer cells than in predivisional cells [63]. However, DivK and DivK;P are dynamically localized during the cell division cycle [63][64][65][66][67][68]. Membrane-bound proteins DivJ and PleC, which localize at stalked and flagellated cell poles, respectively, regulate this process [64,65] by having opposite effects on DivK phosphorylation. DivJ is a kinase that continuously phosphorylates DivK at the stalked cell pole, and PleC promotes the continuous dephosphorylation of DivK;P at the flagellated cell pole [64,67]. Hence, opposing gradients of DivK and DivK;P are established between the two cell poles. Full constriction of the Z-ring disrupts the diffusion of DivK between the two poles [59,64]. As a result, DivK;P accumulates in the nascent stalked cell compartment and unphosphorylated DivK accumulates in the nascent swarmer cell compartment. High DivK;P promotes CtrA degradation in the stalked cell compartment [42,43], whereas high CtrA is maintained in the swarmer cell compartment [16]. The nonuniform distribution of DivK and DivK;P, and their corresponding effects on CtrA degradation, contribute largely to the different developmental programs of swarmer and stalked cells in C. crescentus. In addition, recent investigations indicate that CtrA phosphorylation is also at least partially under the control of DivK;P (as mentioned above), which shows that DivK;P not only controls the stability of CtrA, but also its activity [44,45].

Results
Using Figure 2, we create a wiring diagram (Figure 3) of the molecular interactions that we deem to be most important for regulation of the division cycle in stalked cells. The diagram is then converted into a mathematical model (Table  1), as described in Materials and Methods. Full details of the model can be found on our Web site (http://mpf.biol.vt.edu/ research/caulobacter/pp/), including machine-readable versions of the model (for MATLAB and SBML) and an online simulator.

The Model Accurately Describes Protein Expression Patterns during the Division Cycle of Wild-Type Cells
To simulate the molecular regulation of a wild-type stalked-cell division cycle, we solve the equations in Table 1 subject to the parameter values and initial conditions in Tables 2 and 3. Figure 4 illustrates how scaled protein concentrations and other variables of the model change during repetitive cycling of a stalked cell. The duration of a wild-type stalked-cell division cycle in our simulations is 120 min (;90 min for S phase and ;30 min for G2/M phase), as typically observed in experiments [22,23,59].
The main physiological events of the division cycle can be traced back to characteristic signatures of protein expression, as described in the Introduction. The division cycle starts with initiation of DNA replication ( Figure 4A) from a fully methylated origin site by elevated DnaA, when CtrA is low and GcrA is sufficiently high (to induce production of required components of the replication machinery) ( Figure  4C and 4D). Immediately after DNA replication starts, C ori is hemimethylated.
As DNA synthesis progresses, certain genetic loci become hemimethylated in order along the chromosome ( Figure 4B). Consequently, the regulatory proteins are produced and reach their peak concentrations sequentially. By contrast, dnaA expression seems to be activated by full methylation [55], so its expression rate declines immediately after DNA replication starts. The effect of methylation on dnaA expression is minor compared to the regulatory signals coming from GcrA and CtrA. When the replication fork passes the ccrM locus, the gene becomes available for transcription, but is not immediately expressed, because CtrA level is low. In a predivisional cell, at approximately 35 min after start of DNA replication, the replication fork passes the ctrA gene ( Figure 4B), and its expression is immediately activated by GcrA ( Figure 4C) and then further up-regulated by CtrA itself. Later on, when CtrA level becomes high, expression of the ccrM gene and, later, hemimethylated fts genes (at ;65 min), are expressed by the activation from high-level CtrA ( Figure 4D).
High CtrA down-regulates gcrA expression. When DNA replication is finished, the new DNA strands are methylated by elevated CcrM in about 20 min. DNA methylation shuts down production of CtrA, CcrM, and Fts proteins. Meanwhile, elevated Fts proteins promote Z-ring formation and constriction ( Figure 4D), which separates the predivisional cell into two compartments, thereby restricting access of DivK and DivK;P to only one of the old poles of the cell. As a result, in the stalked cell compartment, most DivK is converted into DivK;P, accelerating CtrA proteolysis there ( Figure 4C). In a nascent stalked cell, low CtrA concentration releases gcrA expression, and GcrA protein level rises. Then, low CtrA, high GcrA, and high DnaA drive the nascent stalked cell into a new round of DNA synthesis from the fully methylated chromosome. These computed properties of the model agree reasonably well with what is known (or expected) about cell cycle progression in C. crescentus.
In Figure 5, we compare our simulation with experimental data. The data, collected from literature, were obtained by different research group with various experimental techniques. In most cases, experimental uncertainties of the data were not reported, but it is reasonable to assume that the error bounds are quite generous. Therefore, based on a visual comparison, we conclude that our model is in reasonable agreement with experimental observations.
The only serious objection that might be raised is to our simulation of DivK;P ( Figure 5C, green curve), which increases rapidly in the stalked-cell compartment after the Z-ring closes and DivK;P is cut off from its phosphatase at the swarmer cell pole. Jacobs et al. [62] reported roughly constant levels of DivK;P in predivisional stalked cells, i.e., until just before Z-ring constriction, and significant differences of DivK;P levels between stalked cells and swarmer cells. Our waveform for DivK;P is consistent with this report and predicts that there should be a distinct peak of DivK phosphorylation in the stalked cell compartment at the end of the division cycle. This peak seems to be an inevitable consequence of the current belief that, upon Z-ring constriction, DivK becomes dephosphorylated in the swarmer cell compartment and remains heavily phosphorylated in the stalked cell compartment.

The Model Agrees with the Phenotypes of Mutant Strains
The phenotypes of mutant cells provide crucial hints for deciphering the biochemical circuitry underlying any aspect of cell physiology. A mathematical model must be consistent with known phenotypes of relevant mutants. To make this test, we simulate cell cycle mutants of C. crescentus using exactly the same differential equations, parameter values, and initial conditions as for wild-type cells (Tables 1, 2, and 3), except for those modifications to parameters dictated by the nature of the mutation (Table 4). Our simulations of 16 classes of mutants are in agreement with experimentally observed phenotypes, as described here.
DctrA-P1. This mutant is obtained by deleting the normal ctrA gene and placing a copy of ctrA near the terminus of the chromosome, where it is almost fully methylated throughout the cell cycle [40]. In this case, expression from the P1 promoter is greatly reduced, but expression from the P2 promoter is still possible [40]. As a result, a large proportion of mutated cells are longer than wild-type cells, and the appearance of the CtrA protein in the mutant cells is delayed relative to wild-type cells [40]. The mutant cells are able to The double-stranded closed curve at the top represents DNA. C ori is the origin of DNA replication and Ter stands for the termination site. All proteins (ovals) are assumed to be produced and degraded at specific rates. Only the degradation of CtrA is shown on the diagram (four small circles indicate products of CtrA degradation), in order to indicate how this step is regulated by closure of the Z-ring (Zclosed) and subsequent phosphorylation of DivK. Dashed lines denote regulatory effects among the components. DNA methylation sites on genes are marked by cyan stars. doi:10.1371/journal.pcbi.0040009.g003 divide, but with some time delay compared to wild-type cells. In our simulation ( Figure S1), CtrA reappears with approximately 30 min delay, compared to wild-type cells. After that, the cell finishes the division cycle, as observed [40].
DctrA. Experiments show that the CtrA depletion strain, DctrA::spec þ P xylX -ctrA (the genomic copy of ctrA is disrupted, and a wild-type copy of the gene under control of a xyloseinducible promoter is integrated into the genome), is arrested and becomes filamentous without xylose inducer [46,47]. In our corresponding simulation ( Figure S2), insufficient CtrA fails to stimulate expression of ccrM and fts genes. Therefore, the DNA remains hemimethylated, and the Z-ring cannot constrict. Cell division, but not metabolism and growth, is blocked, so the cell is expected to become filamentous. ctrA401 ts . When the temperature-sensitive strain, ctrA401 ts , is shifted to the restrictive temperature (37 8C), ctrA transcription is reduced by approximately 50%. This mutant has been studied thoroughly [15,29,30,32,36,47]. In our simulation ( Figure S3), reduction of CtrA production to 70% of wildtype level leaves progression through the cell cycle unaffected. However, if ctrA expression lies between 70% and 35% of wild-type level, then cell division fails, but DNA synthesis continues repetitively, producing cells with multiple chromosomes. This curious phenotype was observed experimentally [30,36]. Lowering ctrA production rate (in our model) leads to decreasing expression (compared to wild-type) of ccrM and fts genes, with the latter becoming insufficient for Zring formation, as observed experimentally [29,32,41,47]. Since the Z-ring remains open, DivK;P does not accumulate  m;cori þ ½CcrM 4 ½h cori ; when ½Elong ¼ P elong ; ½h cori ¼ 1 k ¼ rate constants (min À1 ), J ¼ binding constants (dimensionless), h ¼ thresholds (dimensionless), P ¼ position of genes relative to origin site (dimensionless). a ¼ activation, d ¼ degradation, i ¼ inactivation, s ¼ synthesis, trans ¼ transformation (phosphorylation and dephosphorylation in our case).
[I] ¼ activity of an intermediate component or pathway that introduces a time delay between CtrA activation of ccrM transcription and the accumulation of CcrM protein.
[Fts] ¼ a lumped representation of Fts series proteins as required for Z-ring constriction.
[Ini] and [Elong] represent the initiation of DNA replication forks and the progression of forks along the DNA molecule, respectively.
[DNA] ¼ the total amount of DNA in the cell, scaled to the amount of a full chromosome. The total DNA starts to increase when initiation finishes and is updated with same speed as the elongation [Elong] grows. Both the elongation and the total amount of DNA increase with possible multiple chromosome replication in specific mutants. Count ¼ a number of chromosomes in the cell. Count is used here in case the multiple initiation happens from several chromosomes when cell does not divide (initiation then is assumed to happen simultaneously on all chromosomes). When initiation is successfully accomplished, ''Count'' is doubled.
[h cori ], [h ccrM ], [h ctrA ], and [h fts ] represent probabilities of hemimethylated states of C ori and dnaA, ccrM, ctrA, and fts genes, respectively, during the cell division cycle. All other variables represent scaled protein concentrations. Notes on equations: The rates of protein expression from ctrA promoter can be written as a product of the probability of transcription factors to be bound and the probability that the locus is hemimethylated. The same assumption is also adopted for dnaA, fts, and ccrM expression. Specifically, dnaA expression is assumed to be half-decreased when the gene is hemimethylated reflecting the experimental facts that dnaA has a peak expression only during the swarmer-to-stalked cell transition and that corresponding protein level during the stalked cell cycle is almost constant [55]. Furthermore, as dnaA is very close to C ori , and in order to avoid introducing an extra variable [h dnaA ], we use (2 À [h cori ]) term to describe the methylation effect on DnaA production rate. and CtrA stays elevated longer than in a wild-type cell, therefore lowering GcrA. The two latter effects result in DnaA being elevated so significantly that DNA replication is initiated in an undivided cell, despite an elevated level of CtrA and a depressed level of GcrA. When ctrA expression is reduced below 35% of wild-type level, the simulation is similar to DctrA ( Figure S2).
ctrA op and ctrAD3. A number of different mutations can cause increased levels of CtrA in cells: by constitutive expression of the gene, by producing a poorly degraded form of CtrA, by producing a constitutively active form of CtrA (not needing to be phosphorylated), or by combinations of these mutational strategies [30,38]. The phenotypes of these different mutant strains differ widely, from normal cell cycling to mixed arrest (in G1 and G2/M phases) to dominant G1 arrest. Although our model is consistent with some of these phenotypes, it is not consistent with them all because it does not take into account CtrA phosphorylation (see model assumptions in Materials and Methods).
Overproduction of CtrA does not interfere with normal cell cycling: the DctrA1::spec þ P xylX -ctrA mutant grows and divides normally when CtrA is expressed constitutively from a xylose-inducible gene on a high copy-number plasmid [38]. By contrast, in our model ( Figure S4), these cells arrest in G1 because they produce CtrA too early in the cell cycle. In reality, although G1 cells have an elevated level of CtrA, it is inactive (unphosphorylated), so these cells can proceed normally into S phase. This mutant emphasizes the importance of regulated CtrA phosphorylation. Since we do not yet account for CtrA activation and inactivation, this mutant is beyond the scope of our present model.
When the genomic copy of ctrA is missing the coding sequence for the last three amino acids (ctrAD3X), the encoded mutant CtrA protein is more stable [38]. When this gene is introduced on a high copy-number plasmid (ctrA þ [wild type] þ P xylX -ctrAD3) and the cells are grown on 0.2% xylose to overexpress the stable form of CtrA, then the mutant cells become filamentous, arresting either in G1 phase (unreplicated DNA) or in G2 phase (replicated DNA) [30]. Our model exhibits this behavior; whether cells arrest in G1 or G2 depends on cell cycle phase when CtrA production is induced ( Figure 6). If CtrA production is induced in early-tomid S phase, it does not interfere with DNA elongation, but does stimulate ccrM expression earlier than in wild-type cells, and the DNA becomes fully methylated promptly. Consequently, fts gene expression is extremely low, and cells become arrested with 2n chromosomes (Figure 6-I). On the other hand, if CtrA production is induced after fts genes have been turned on (Figure 6-II), then cell division proceeds on schedule, CtrA is not removed by proteolysis, and cells arrest  in G1, unable to initiate a new round of DNA replication, because CtrA is blocking C ori . DgcrA. In the GcrA depletion strain, DgcrA:: X þ P xylX -gcrA, the chromosomal gene is disrupted and the gcrA coding sequence (under control of a xylose-inducible promoter) is integrated into the chromosome. In this mutant, when xylose is removed from the medium, the CtrA protein and divK transcription levels are decreased and dnaA expression is increased [15]. In our simulations ( Figure S5), CtrA and DivK protein variations follow the experimentally observed trends after the first cycle. DnaA protein is elevated due to release from GcrA inhibition, as expected. In the simulation, GcrA depletion was introduced at t ¼ 120 min, when DNA synthesis has already been successfully initiated. Since in our model progression of DNA replication is independent of GcrA, it persists until the end of the already initiated cycle. In real cells, GcrA is involved in maintenance of the replication machinery. Therefore, whether DNA replication will finish after the deletion of the gcrA gene depends on the abundance of GcrA protein relative to its rate of degradation. We are unaware of experimental studies relevant to this issue.
DdnaA. When the dnaA:: X þ P xylX ::dnaA mutant is shifted from xylose to glucose medium in order to block production of DnaA, DNA synthesis is arrested, and cell division is blocked [33,51]. Consequently, the expression of many genes is blocked. But when the cell is shifted from glucose back to xylose, DNA replication resumes, and the cell returns to a typical wild-type morphology. In our model mimicking this mutant, initiation of a new round of DNA synthesis fails due to insufficient DnaA protein ( Figure S6). As a result, genes that are expressed only from hemimethylated DNA are never transcribed and the cell arrests in G1. In our simulation, the cell is viable when it is in the xylose medium, as dnaA is expressed constitutively (Figure S7, 120 , t , 300 min). If this simulated cell is shifted to glucose medium for a while ( Figure  S7, 300 , t , 370 min) and back to xylose medium again (t . 370 min), then it returns to normal cell cycling ( Figure S7), in agreement with experimental observations [18].
DftsQ and DftsZ. In the DxylX-ftsQ::pBGPxQ mutant, the chromosomal copy of ftsQ is deleted by recombination with a plasmid carrying the ftsQ gene under control of the P xylX promoter. Blocking FtsQ production, by growing the cells on glucose, results in filamentous cells [61]. In our model, ftsQ is assumed to be essential for Z-ring formation and/or constriction. In the simulation ( Figure S8), after fts expression is turned off, a cell can proceed through one round of DNA replication, but the Z-ring never fully constricts, and the cell never divides. Because growth continues (CtrA level is high), the cell is expected to become filamentous.
Similarly, when ftsZ is deleted (ftsZ163DC þ P xylX -ftsZ), cytokinesis is inhibited, and cells become filamentous [67]. Our simulation ( Figure S9), which assumes no Z-ring formation and constriction in the DftsZ mutant, produces results similar to those for DftsQ and is in agreement with observations. DdivK. Deletion of divK (DdivK::Spec R ) was reported to be lethal [68,69], although how cells die when DivK is depleted has yet to be determined. In our simulation (Figure 7), after blocking expression of divK, cell divisions continue for some time, as divK protein is slowly lost from the cell (half-life ¼ 350 min). Eventually, DNA replication and cell division stop with a high level of CtrA in the cell. In the terminal state, our simulated cell contains two DNA copies. The simulated cell eventually fails to divide because of insufficient levels of Fts proteins ( Figure 7D). Compared to wild-type cells, the CtrA level is elevated, which stimulates production of Fts proteins and CcrM. The latter quickly methylates DNA, leaving a narrower window for expression of Fts proteins ( Figure 7B). We were not able to find any experimental information about the DNA content of this mutant.
The model predicts that cells depleted of DivK may undergo four to five division cycles before dying, because (in the model) DivK protein is only slowly degraded, and its only job is to trigger degradation of CtrA. To the extent that  Figure 4A in [40], compares well with our simulation (red and blue curves). doi:10.1371/journal.pcbi.0040009.g005 other processes may render DivK nonfunctional and that DivK may have other essential roles in the cell [42,68], real cells may lose viability after DivK depletion more quickly than predicted here. divK341 cs . A point mutation (D90G) of divK creates a coldsensitive allele, divK341 cs , that (at the restrictive temperature, 25 8C) maintains a constant high level of CtrA and does not initiate DNA replication. Cells elongate, grow stalks, and become arrested in G1 phase with one chromosome [16,42,69]. In our simulation (Figure 8), we assume that, immediately after transfer to the restrictive temperature, DivK;P loses its capacity to stimulate CtrA degradation [42]. After the origin of replication is initiated, DNA synthesis and methylation proceed normally in this case. The Z-ring fully constricts, but CtrA stays elevated. After the first division, DNA synthesis and cell division cease (in our simulation) because CtrA cannot be degraded. (The model also predicts a delayed initiation of DNA replication at t . 400 min, not evident in Figure 8. This is an artifact of our oversimplified differential equation for the DNA initiation variable, [Ini].) If the mutation is introduced in our simulation before DNA replication starts, its start will be delayed because it takes longer for CtrA to be reduced (by the background degradation rate alone) to the replication-permissible level ( Figure S10). Importantly, in the terminal state, the cell will block in G2 with replicated DNA.
If divK341 cs is returned to the permissive temperature within about 70 min, cells recover to normal division cycles in our simulations ( Figure S11; mutant is put into nonpermissive temperature at t ¼ 120 min and returned to permissive temperature at t ¼ 190 min), as in experiments [42].
DivK D53A , DivK E9A , and DivK D10A . In this family of DivK amino acid replacement mutations, the phosphorylation of DivK is impaired, and DivK protein remains homogeneously dispersed through the cell [68]. In experiments, these mutant cells arrest in G1 phase and became filamentous [68]. These characteristics are observed in our simulations ( Figure S12), which are similar to the simulation results (Figure 8) for the divK341 cs mutant. divK op . In experiments, a divK-cfp fusion gene under control of the xylose-inducible promoter on a medium copy-number plasmid was introduced into wild-type cells [68]. When overexpressed, DivK was mislocalized, cell division was blocked, and cells became filamentous [68]. In our simulation, several-fold overeproduction of DivK causes arrest in G1 phase after one or more division cycles, depending on how much DivK is produced (Figure 9). Cells with lower divK overexpression undergo more cell cycles before arresting. These simulations are consistent with the observed phenotypes of divK op mutants. (Of course, mislocalization of DivK may have other effects on the cell cycle [68] that are not taken into account in our model.) In [49], DivK was overexpressed from a high-copy plasmid. In our simulation, when DivK is strongly overexpressed (e.g., 20-fold), the cell becomes arrested in G2 phase. The lower level of CtrA is able to stimulate enough synthesis of DnaA to initiate a new round of DNA replication but fails to increase Fts sufficiently to trigger cell division. Thus the 20-fold overexpression mutant has a terminal phenotype similar to ctrA ts ( Figure S3).
DccrM. In a CcrM-depleted strain (DccrM þ P xylX -ccrM, where the chromosomal ccrM locus is inactivated and CcrM is supplied from a xylose-inducible promoter on a low copynumber plasmid), methylation of the new DNA strand ceases, when cells are grown on glucose [39,40]. Cell growth and DNA replication cease in 6-8 h. Cell viability starts dropping after approximately 4 h. Our simulation of this mutant ( Figure 10) shows that, without sufficient CcrM, DNA remains hemimethylated; hence, initiation of a new round of DNA replication is repressed. In addition, ctrA and fts genes can now be periodically expressed in the absence of DNA replication, which is usually required to return these genetic loci to a hemimethylated state. Elevated, periodic variations of Fts proteins can induce premature constriction of the Zring, if other relevant conditions are also satisfied. That leads to cell division without DNA replication, suggesting a precipitous loss of viability of these mutant cells, as observed experimentally [39,40]. ccrM op . In the CcrM overproduction mutant, the wild-type ccrM gene is supplemented by a second chromosomal copy expressed constitutively from a P lac promoter. CcrM overproduction causes rapid methylation of newly synthesized strands of DNA in experiments; some cells accumulate multiple chromosomes because additional initiations of DNA synthesis occur before cell separation [52]. About half of the mutant cells are longer than wild-type, and cell division is morphologically aberrant [52]. In our simulation ( Figure  11), a constitutively high level of CcrM accelerates methylation of newly synthesized DNA; hence, ctrA has little chance to be transcribed, and consequently, the Z-ring never fully constricts. At the same time, conditions are right for repeated rounds of DNA synthesis. The number of excess chromosomes per cell depends on the extent of overexpression of ccrM in our simulations, as can be seen by comparing Figure 11 (2-fold overexpression) with Figure S13 (50% overexpression).
Lon null mutant. Lon protease is required for CcrM
DdnaC and dnaE ts . Elongation of DNA during S phase ceases in cells depleted of DnaC [70]. The dnaE ts mutant also blocks DNA elongation at the nonpermissive temperature (37 8C), and cells arrest with an undetectable level of CtrA [41]. To simulate these mutants, DNA elongation was interrupted at different times during the cell cycle ( Figure S15). If DNA elongation is blocked in early-to-mid S phase, then few fts genes can be successfully expressed, and the Z-ring stays unconstricted. The mutant cell arrests in an early S phase (GcrA high, CtrA low). If the knock-out is made late in S phase, then cell division occurs normally, provided that incompletely replicated chromosomes can be separated to progeny cells. The progeny cells would then arrest pre-S phase. In reality, incompletely replicated chromosomes may prevent cell division [71], in which case, cells would arrest in a late predivisional stage. We do not include an unreplicated-DNA checkpoint in our model.

The Model Predicts Phenotypes of Novel Mutants
DctrA-P2. Based on our simulations (Figure 12), the cell cycle in a ctrA-P2 deletion mutant should arrest in G2 (replicated chromosome) as long (filamentous) cells. CtrA never reaches a high concentration and hence fails to activate much expression of ccrM and fts genes. Without sufficient Fts proteins, the Z-ring does not constrict. At the same time, DNA replication can continue, once it has been initiated. However, new DNA strands cannot be fully methylated because of insufficient CcrM, and thus a new round of DNA replication cannot be initiated.
gcrA op . No mutant overproducing GcrA has been reported in the literature. Our model predicts that small overexpression of gcrA (;110% in our simulations) speeds up production of CtrA, which accelerates progression through the cell cycle ( Figure S16). Higher levels of GcrA (150%-200% overexpression in our simulations) cause a significant decrease of DnaA protein. As DnaA level drops, eventually DNA synthesis cannot be initiated, and the cell should arrest in G1 phase ( Figure 13). dnaA op . Phenotypes of dnaA overexpressing mutants have not been reported in the literature, to our knowledge. According to our simulations, modest overexpression of dnaA (110%) stimulates gcrA expression quicker than in wild-type cells. High levels of DnaA and GcrA, combined with low CtrA, accelerate initiation of DNA replication and speed up the cell cycle somewhat (;100 min in our simulation; Figure S17). If dnaA is overexpressed by 150% or more, the elevated level of DnaA protein causes rapid, repeated initiations of DNA replication forks, in our simulations (Figure 14), suggesting that at least several initiations of DNA replication may take place. Our model permits instantaneous initiation of DNA replication as soon as the initiation condition is satisfied, whereas in real cells, re-initiation must take some minimal time. Hence, our prediction is that several re-initiations of

Discussion
We propose (Figure 3) a realistic mechanism for regulating the cell division cycle of stalked cells of C. crescentus. The mechanism includes three master-regulatory proteins (GcrA, DnaA, and CtrA), a DNA methylase (CcrM), Z-ring components (Fts proteins), and an end-of-cycle protein (DivK) in its inactive and active (phosphorylated) forms. Cytokinesis is represented by a phenomenological variable that describes the extent of constriction of the Z-ring. DNA synthesis is described in terms of initiation, elongation, and termination. We assume that initiation of DNA replication requires high DnaA and GcrA, low CtrA, and full methylation of the origin  site, and that the rate of DNA elongation is independent of DnaA, GcrA, and CtrA, and is almost linear. Transcription of some genes occurs only from an unmethylated DNA sequence; hence, the expression of such genes depends on their location on the newly synthesized DNA strand. Compartmentation in the predivisional cell is assumed to result in localization of phosphorylated DivK to the stalked compartment of the dividing cell, promoting CtrA degradation there.
These assumptions are formulated as a mathematical model (Table 1) consisting of 16 nonlinear, ordinary differential equations for seven proteins, the state of the Z-ring, the progression of DNA synthesis, and the methylation state of five gene sites on the DNA. The rate equations entail 44 parameters (rate constants, binding constants, and thresh-  olds; Table 2) that need to be determined by fitting the model to specific experimental observations. For the present, parameter estimation is done by trial and error, so we can only claim that our model equations and parameter set are sufficient to account for many properties of cell cycle control in C. crescentus. Because we fit the model to many different mutant phenotypes, we have a wealth of data to fix the parameters and to provide meaningful confirmation of the mechanism. Table 2 is in no sense an optimal parameter set, nor can we quantify how robust the system is, although our experience suggests that the model is quite hardy.
Our present model is based heavily on an earlier conjecture [17] that the C. crescentus cell cycle is controlled by a bistable switch, created by positive feedback in the molecular circuitry of the ctrA gene. In that conjecture, the switch is flipped from the off-state (CtrA low) to the on-state (CtrA high) by GcrA  accumulation as cells enter S phase, and then switched back to the off-state by DivK activation (phosphorylation) as cells divide (the CtrA-DivK negative feedback loop). The original model did not account for the ways in which gene expression is linked to DNA methylation, thereby anchoring the protein interaction network to the progression of DNA replication forks. By incorporating DNA synthesis and methylation into the Brazhnik-Tyson model, the present model provides a more satisfactory account of cell cycle regulation in C. crescentus, and it can be tested by comparison to a broad spectrum of mutant phenotypes. Because the new model successfully reproduces the behavior of wild-type and mutant cells in many quantitative details, we conclude that our present understanding of the control system ( Figure 3 and Table 1), properly interpreted, is accurate and adequate. On the other hand, the proposed mechanism must be considered as an evolving hypothesis that will be continually examined, revised, and improved as new observations tell us more about the control system. Some obvious improvements to the model include refined criteria for DNA initiation, regulated phosphorylation of CtrA, spatial localization of proteins, inclusion of a swarmer cell compartment, and an account of the swarmer-to-stalk cell transition.
Finally, most of division-control proteins (such as CtrA, DivK, CcrM, FtsZ, and FtsQ) are conserved among aproteobacteria [72], suggesting that the computational model proposed here for C. crescentus may prove applicable to other types of a-proteobacteria, including symbiotic nitrogenfixing genera (Rhizobia) and pathogenic genera (Brucella spp., Coxiella spp, etc.)

Materials and Methods
Scope of the model. To understand the molecular logic of cell cycle regulation in C. crescentus, we constructed a mathematical model of the temporal dynamics of the regulatory genes and proteins. Following standard rules of chemical kinetics, we converted the wiring diagram in Figure 3 into a set of rate equations describing the temporal dynamics of the model. Justification of our approach is described in detail in [17].
Two phenomenological variables, Z (the state of closure of the septal Z-ring) and I (introducing a delay between activation of ccrM transcription and later activation of CcrM protein production).
The progression of DNA replication (including initiation, elongation, and termination) and its methylation (including probabilities of hemimethylation of ccrM, ctrA, dnaA, and fts genes, and of the replication origin site, C ori ).
Accordingly, our mathematical model consists of 16 nonlinear differential equations presented in Table 1, including 28 kinetic constants (k's), 11 binding constants (J's), and five thresholds (h's). Our choice of parameter values is given in Table 2.
Assumptions of the model. A common trend in developing complex models in molecular cell biology is to start from a simple coarse-grained (''phenomenological'') model and then refine and expand it step by step (as data become available) into an increasingly more comprehensive model. (A good example is the progression of models of the budding yeast cell cycle [2,4,73].) We have taken this approach in our study of the C. crescentus cell cycle. We have limited the scope of our model so that it can be based largely on experimental observations, is not overwhelmed with assumptions, and is able to make predictions. Obviously, at any stage of modeling there will be facts that have not yet been incorporated and thus are out of the scope of the model. Our modeling assumptions are described here.
First, we propose to model, at this stage, only the average behavior of cells and do not address naturally occurring fluctuations in cell cycle progression.
Second, the rise of DivK;P in stalked compartments after constriction of the Z-ring is a necessary, but not sufficient, condition for CtrA degradation. In our coarse-grained model of CtrA proteolysis, we use DivK;P as a signal for starting rapid degradation of CtrA. In other words, DivK;P determines when the degradation of CtrA is turned on, but the how (the machinery that degrades CtrA, involving RcdA, CpdR, and ClpXP) is assumed to be there when needed and is not modeled at present.
Third, CtrA is activated by phosphorylation (by kinases CckA and DivL), and a complete model of the Caulobacter cell cycle should take this into account. Unfortunately, little is known about the phosphorylation and dephosphorylation of CtrA and how these processes are temporally regulated. During the division cycle of wild-type cells, the levels of CtrA and CtrA;P rise and fall together [22,46], so we need not distinguish between the two forms. Therefore, in the current model, we keep track of CtrA synthesis and degradation only, assuming that CtrA;P is a fixed fraction of total CtrA. This assumption, though a great oversimplification, is harmless enough for most of the mutants we consider in this paper. But it seems to cause serious problems for exactly those mutants (ctrA op , ctrAD3, ctrAD51E, and ctrAD51ED3 in wild-type background) that interfere with normal synthesis, degradation, or activation of CtrA [34].
Later versions of the model will have to include CtrA;P as a variable, when we have a better of idea of the mechanisms controlling CtrA phosphorylation.
It is known that DivK;P promotes the proteolysis of CtrA;P [42] and negatively regulates CckA activity, thereby reducing phosphorylation of CtrA [49,74]. Hence, DivK;P works to eliminate CtrA;P activity by two independent pathways. We lump these two effects together as a single DivK;P promoted reaction for removing active CtrA.
Fourth, the dnaA locus is very close to the origin site (C ori ) [28]. Within its promoter, potential CtrA and DnaA boxes and methylation sites exist for regulating its expression [20,34,52,55]. GcrA is a repressor for dnaA expression [15], and CtrA seems to be an activator [32]. However, DnaA protein concentration varies very little during the Caulobacter cell cycle [55]. Although we include the regulatory signals in the model, they do not much affect the dynamics of a stalked cell because DnaA level is nearly constant throughout the cell cycle due to DnaA's long half-life.
Fifth, initiation of DNA replication is triggered by the combined conditions of low CtrA, high DnaA, and fully methylated DNA origin site. In addition, initiation requires sufficient replication machinery, which is correlated to a high level of GcrA. We combine these regulatory effects into a single term. We assume that once initiation of DNA replication is successful, DNA elongation starts immediately. Elongation of new DNA strands is linear in time until it finishes, based on experimental data indicating that the speed of DNA replication in C. crescentus is almost constant [54].
Sixth, full constriction of the Z-ring requires accumulation and activation of a number of proteins, including FtsZ, FtsQ, FtsA, and FtsW, some of which are stimulated by CtrA. To simplify the model, we use Fts as a combined component to relay the signal from CtrA to Z-ring constriction. The transition from Z-ring open (¼ 1) to fully constricted (¼ 0) is modeled as a Goldbeter-Koshland ultrasensitive switch [75].
Seventh, we include the effects of DNA methylation on gene expression in our model because these effects mediate important feedback loops between DNA synthesis and the master regulatory proteins, and because DNA methylation can be a useful target for new drug development. In our model, the genes ccrM, dnaA, ctrA, and fts as well as the origin of DNA replication are regulated by methylation.
Methylation plays a minor role in the regulation of GcrA production [19], so we disregard it in our model. We allow a modest contribution of DNA methylation to regulating the production of DnaA. ccrM gene expression is significantly affected by its methylation state [40,57]. The activity of ctrA-P1 is known to depend on hemimethylation [36], and the activity of ctrA-P2 seems to depend in some other way on DNA replication [37]. For simplicity, we assume that both ctrA promoters are turned on by hemimethylation of the gene.
Among fts genes, the ftsZ promoter has a methylation site [40,53], but the ftsQ promoter does not [41]. Scanning the ftsQ gene for the consensus sequence GANTC using the Regulatory Sequence Analysis Tools (http://rsat.ulb.ac.be/rsat/), we found a GAGTC segment in the coding sequence, suggesting that the ftsQ gene might also be affected by methylation. Since our ''Fts'' variable is a combination of Fts proteins, we conclude that our fts gene should be regulated by methylation.
The effects of methylation on gene promoters and C ori are described by probabilities to be methylated or hemimethylated during the cell cycle. The probabilities (h .. variables) are in turn controlled by the progression of DNA replication and by the activity of CcrM [52,53].
Eighth, ccrM transcription is tightly regulated by CtrA protein, but accumulation of CcrM protein shows a noticeable delay from the transcriptional activation of its gene [37], resulting in delayed activation of DNA methylation [57]. This delay is mimicked in our model by an intermediate variable I in the CtrA-to-CcrM pathway.
Ninth, we recognize the importance of spatial controls in the Caulobacter cell cycle. However, at this stage, we are trying to model the stalked cell cycle as far as possible without explicitly tracking the spatial localization of regulatory proteins. That would require a more sophisticated mathematical framework and is planned for the next stage of the model. As the result of this simplification, our model makes no distinction between the stalked and swarmer parts of the predivisional cell. Right after compartmentation and before cytokinesis, we keep track of proteins in the stalked cell compartment only. At this stage, the distinction between swarmer and stalked cells is made by the phosphorylation state of DivK (being completely phosphorylated in the stalked compartment).
Tenth, we assume cells grow steadily in time, with a mass-doubling time of about 120 min and with the accumulated material shed at each division in the swarmer cell. In the present model, there is no coupling between cell growth and division, as in our models of eukaryotic cell proliferation [10]. Hence, there is no need for us to keep track of cell size, except to notice that if cell division is delayed or blocked, then the stalked cell will grow longer than normal and eventually be described as having a filamentous morphology.
Parameter values and initial conditions. Parameter values for our model ( Table 2) were determined from available experimental data, wherever possible. Rate constants of degradation were estimated from experimentally observed half-lives of proteins. Rate constants of protein synthesis were adjusted to fit variations of protein concentration observed in experiments. Parameter values of Z-ring dynamics were set to be consistent with observed durations of the open (;100 min) and constricted (;20 min) states of the Z-ring [59]. Rate constants of DivK phosphorylation and dephosphorylation were estimated from the difference of DivK;P concentration before and after Z-ring closing in predivisional cells [63]. Successful initiation of DNA replication depends on satisfying four requirements: low CtrA ([CtrA] , h CtrA ), high DnaA ([DnaA] . h DnaA ), high GcrA ([GcrA] . h GcrA ), and a fully methylated origin site (h cori , h Cori ). The thresholds were adjusted to position the onset of the S phase correctly in wild-type cells. Replication-fork progression (elongation) begins at each successful initiation ([Ini] ¼ 0.05) and stops when DNA replication is complete ([Elong] ¼ 1).The constant rate of elongation is consistent with an 80min delay for copying the chromosome. Due to the constant rate of DNA replication, those genes that must be hemimethylated in order to be transcribed will be expressed in a temporal sequence determined by their positions on the chromosome from the origin of replication [39,76]. To model this effect, the variable h gene is set to 1 (hemimethylated) when [Elong] ¼ distance of gene from C ori . Some time thereafter, when CcrM activity is high, the h gene decays exponentially back to 0 (fully methylated). Most Hill function exponents are assumed to be 2, with a higher value (n H ¼ 4) where sharper switching was required. Initial conditions (Table 3) were taken to represent the beginning of a stalked cell cycle in a wild-type cell.
Simulation of mutants. The phenotypes of relevant mutants were collected from the literature. To simulate each mutant, we use exactly the same equations (Table 1) and parameter values ( Table 2) except for values of those parameters directly affected by the mutation (Table 4). Mutations are introduced in our model after 120 min of simulation of the wild-type cell. For gene deletion, the rate of synthesis of the corresponding protein is set to zero. For gene overexpression, an additional constant rate of synthesis of the corresponding protein is introduced into the equations, because proteins are typically overexpressed from an extra copy of the gene under control of an inducible promoter. For heat-or cold-sensitive mutants, the relevant rate constant(s) retains its wild-type value at the permissive temperature and is set to zero at the restrictive temperature. For partial deletions, the relevant parameter value is assumed to lie between 0% and 100% of the wild-type value, according to the experimental characterization of the mutation.
Equations of the model were solved numerically with Matlab 2006a (The MathWorks). Machine-readable files for reproducing our simulations are made available in Text S1 and on our Web site (http://mpf.biol.vt.edu/research/caulobacter/pp/).