Modeling the temporal dynamics of master regulators and CtrA proteolysis in Caulobacter crescentus cell cycle.

The cell cycle of Caulobacter crescentus involves the polar morphogenesis and an asymmetric cell division driven by precise interactions and regulations of proteins, which makes Caulobacter an ideal model organism for investigating bacterial cell development and differentiation. The abundance of molecular data accumulated on Caulobacter motivates system biologists to analyze the complex regulatory network of cell cycle via quantitative modeling. In this paper, We propose a comprehensive model to accurately characterize the underlying mechanisms of cell cycle regulation based on the study of: a) chromosome replication and methylation; b) interactive pathways of five master regulatory proteins including DnaA, GcrA, CcrM, CtrA, and SciP, as well as novel consideration of their corresponding mRNAs; c) cell cycle-dependent proteolysis of CtrA through hierarchical protease complexes. The temporal dynamics of our simulation results are able to closely replicate an extensive set of experimental observations and capture the main phenotype of seven mutant strains of Caulobacter crescentus. Collectively, the proposed model can be used to predict phenotypes of other mutant cases, especially for nonviable strains which are hard to cultivate and observe. Moreover, the module of cyclic proteolysis is an efficient tool to study the metabolism of proteins with similar mechanisms.

Introduction

Caulobacter crescentus (C. crescentus) is a model organism for exploring cell development and cell cycle regulation in prokaryotes. C. crescentus undergoes an asymmetric cell division producing two distinct progenies: a sessile stalked cell equipped with a stalk and a motile swarmer cell equipped with a flagellum (Fig 1). While the stalked cell immediately initiates chromosome replication and enters the next cell cycle, the swarmer cell searches for suitable environments and differentiates into a stalked cell (sw-to-st transition) before entering the cell cycle replication [1]. The dimorphic lifestyle makes C. crescentus feasible to survive in oligotrophic waters.

Fig 1

The asymmetric cell cycle of C. crescentus including G1, S, and G2 phases.

C. crescentus cell grows in G1, replicates DNA in S phase, and prepares for cell division in G2 phase. The predivisional (PD) cell divides asymmetrically into two different progenies: motile swarmer (SW) cell and non-motile stalked cell (ST). The dynamics of CtrA, GcrA, DnaA, CcrM, SciP, RcdA, and CpdR is indicated by color during each stage of the cell cycle.

The timed asymmetric cell progression of C. crescentus is highly regulated by a cell cycle-dependent regulatory network including four master regulators–DnaA, GcrA, CtrA, and CcrM [2, 3]. DnaA, GcrA, and CtrA are transcriptional factors that control over 200 cell cycle-regulated genes in C. crescentus. These proteins form a loop to control each other. DnaA activates gcrA expression, GcrA regulates the expression of ctrA and dnaA, and CtrA in turn influences the transcription of dnaA [4-6]. Furthermore, DnaA is a conserved DNA replication initiator, activating replication by binding directly with the chromosome origin (Cori) [7]. In addition, there are five binding sites for CtrA on Cori, where replication initiation is suppressed when being bound by the phosphorylated form of CtrA (CtrA∼P) [4]. CcrM, a conserved methyltransferase, is turned on at the completion of DNA replication to fully methylate the motif GANTC, which is carried by promoters of ctrA, dnaA, and ccrM (see Fig 2) [8, 9]. The short window of CcrM allows the maintenance of hemimethylated chromosome during replication, ensuring the robustness of cell cycle development. Moreover, CcrM has been reported to influence the expressions of more than 10% genes [8, 10]. Among these CcrM-regulated genes, more than 100 genes are likely influenced by a GANTC motif-dependent pathways, while the mechanisms of the rest genes are not clear [8]. Here, we take the CcrM-dependent methylation of GANTC motif into the regulatory network although it is dispensable for the replication control. Additionally, SciP is an antagonist of CtrA which is instrumental in cell cycle regulations but receives little attention. SciP spatiotemporally represses the transcription of CtrA-induced genes because most of these genes contain a SciP binding site upstream of a CtrA binding site in their promoters [2].

Fig 2

Methylation site locations of different genes on C. crescentus chromosome.

The elliptical curve represents the DNA fork in replication. Cori is the origin of DNA replication and Ter is the termination site. The CcrM methylation site is located upstream of the dnaA, ccrM and ctrA genes, represented as rectangles.

A wealth of experimental data for cell cycle-regulated genes and proteins in C. crescentus have been accumulated in last decade [11, 12]. System biologists have proposed different quantitative models to analyze underlying mechanisms and pathways of cell cycle regulation. For example, Li et al. [4, 13] quantitatively modeled the interactions between CtrA, DnaA, GcrA, and CcrM and studied the simulated behaviors of some mutants. Murray et al. [10] proposed a simplified model incorporating CtrA, CckA, and GcrA to capture the cell cycle features of C. crescentus and predict the behaviors of ΔgcrA cells. However, the proteolysis of CtrA is not explicitly modeled. Li et al. [4, 13] borrowed DivK while Murray et al. [10] used CckA to describe the proteolysis of CtrA; but both DivK and CckA are indirect factors influencing the proteolysis of CtrA through phosphorelay pathways [14]. There are a series of models working on the spatial regulatory networks in C. crescentus. Li et al. investigated the spatial regulations focusing on CtrA in a stochastic model [15], which preliminarily revealed roles of spatial phosphorylation on the asymmetric cell cycle in C. crescentus. Further, Chen et al. [16] and Xu et al. [17] proposed spatial models for the scaffolding protein PopZ in C. crescentus, which complemented Li et al.’s model [15] about the initial localization factors. Although previous mathematical models revealed some mechanisms of Caulobacter system, a comprehensive model for core regulators of the cell development, as well as a quantitative comparison between simulations and observations, have yet been explicitly investigated. Previous models didn’t consider the mRNA abundance and transcription process based on master regulators. Additionally, there is no mathematical model describing the cyclic proteolysis of master regulator CtrA, which plays important roles for cell development especially for the sw-to-st transition [18].

In this paper, we focus on five core components–DnaA, GcrA, CtrA, CcrM, and SciP that control over 90% of cell cycle-regulated genes, and propose a mathematical model that considers the regulation of DNA replication and methylation, as well as the gene-protein and protein-protein interactions. Since CtrA is essential in the cell cycle regulation and its proteolysis is distinctively and spatiotemporally regulated, we construct a hierarchical ClpXP complex network for its proteolysis, which is then integrated into the entire model. The simulated dynamics of mRNA and proteins are consistent with experimental observations. The ClpXP complex model can be used as a quantitative analysis tool to simulate other cyclic proteolysis in C. crescentus, such as the proteolysis of ShkA and TacA [19, 20].

Materials and methods

Model description

The regulatory network of bacterial cell cycle includes a series of complex mechanisms, such as genetic regulations, degradations, phosphorylation, dephosphorylation, and so on. The details of the complex regulatory network will be described in the following.

Module 1: The core regulatory network of cell cycle

The master regulatory network of C. crescentus, as summarized in Fig 3, is composed of DnaA (dnaA), GcrA (gcrA), CtrA (ctrA), CcrM (ccrM), and SciP (sciP). Specifically, DnaA promotes the expression of gcrA, while GcrA inhibits DnaA and activates one of the promoters (P1) of CtrA [7]. Conversely, CtrA∼P suppresses the initiation of DNA replication [21], activates the transcription of dnaA [4], inhibits the activity of P1, and activates itself through another promoter (P2) [3]. The accumulation of CtrA promotes the expressions of ccrM and sciP, where CcrM controls the methylation state of P1 of ctrA [22]. SciP downregulates CtrA and CcrM [2]. The regulatory network of the five master proteins and mRNAs governs cell cycle-regulated genes, thereby driving the cell cycle progression [11].

Fig 3

The master regulatory network of C. crescentus.

Solid lines represent activation/inhibition influences of master regulators (DnaA, GcrA, CcrM, CtrA, SciP) with arrow/bar, respectively. The dashed lines represent the methylation effects on dnaA, ctrA, ccrM genes from CcrM.

In normal cell cycle progression, active CtrA (phosphorylated form, blue color in Fig 1) is cleared during the sw-to-st transition; CtrA concentrations are generally low in stalked cells when the Z-ring is closed [23]. The activity of CtrA is controlled by synthesis, degradation, and phosphorylation, the latter of which is driven by the CckA-dependent pathway [18, 24] (Fig 3). As CckA∼P is the only known phosphoryl donor of CtrA [14], we involve the CckA-dependent phosphotransfer into our model. CtrA proteolysis depends on a particular protease complex comprising the protease ClpXP and four additional adaptors–CpdR, RcdA, PopA, and c-di-GMP (cdG) [18]. While the protease ClpXP presents throughout the entire cell cycle, RcdA and CpdR co-localize at the stalked pole during sw-to-st transition and stay in the predivisional cell’s stalked compartment (gold and black circles in Fig 1) [25]. The phosphorylation of CpdR is also regulated by CckA, thus CckA regulates the activity of CtrA through both phosphorylation and proteolysis. Additionally, CtrA controls the expression of RcdA and CpdR in C. crescentus.

Module 2: Cell cycle-dependent proteolysis of CtrA

The stability and activity of proteins strictly regulate cell cycle processes. Accordingly, proteolysis plays a significant role in cell development and response to internal/external stimuli [18, 26]. ClpXP, a highly conserved protease, is responsible for the proteolysis of a wide range of proteins including CtrA in C. crescentus [18]. Many substrates of ClpXP are cell cycle-regulated. Although ClpXP levels do not change significantly throughout the cell cycle, it requires additional cell cycle-dependent adaptors to cyclically degrade proteins [18]. Substrates of ClpXP-based proteolysis require different classes of protease complex assemblies [27]. Substrates, such as PdeA, only require ClpXP primed by unphosphorylated CpdR; we name this type of substrates as the first class substrate. Similarly, the second class substrates require primed ClpXP additionally with RcdA assembled; and the third class substrates, such as CtrA, require binding between PopA and c-di-GMP connected with the second class protease complex (see Fig 4).

Fig 4

Hierarchical proteolysis of the first (eg. PdeA), second (eg. TacA), and third (eg. CtrA) substrate.

The degradation of different substrates is dependent on the degree of adaptor assembly. Priming of the protease ClpXP by unphosphorylated CpdR results in PdeA decay, which recruits additional adaptor RcdA to degrade TacA. RcdA tethers cdG-bound PopA with primed ClpXP, which is responsible for the proteolysis of CtrA.

To simulate the protease complex for CtrA degradation, we use ‘Complex 1’, ‘Complex 2’, and ‘Complex 3’ to name the protease complexes that degrade the first, second, and third class substrates, respectively (see Fig 5). Unphosphorylated CpdR primes ClpXP to function as the first class protease complex (Complex 1) which degrades CpdR in turn [14]. (Table 1, Eq. 20). Primed ClpXP (Complex 1) recruits RcdA (Complex 2) to deliver the second class substrates to the protease ClpXP [27] (Table 1, Eq. 23). Additionally, the RcdA proteolysis has been shown to be catalyzed by Complex 1 [20]. Besides CpdR and RcdA, the third class proteolysis requires PopA bound with cdG, where cdG-bound PopA directly interacts with RcdA and CtrA ensuring the specific degradation of CtrA [18]. The diguanylate cyclase PleD and phosphodiesterase PdeA are included in our system as the main synthetase and hydrolase of cdG, respectively, where the expression of pleD and pdeA is controlled by CtrA∼P [28]. PdeA is proteolyzed by Complex 1, shown in Fig 5. As PopA bears a GGDEF domain and two receiver domains akin to PleD, we assume PopA functions as a dimer; thus, PopA dimer binds with two cdG molecules in the same way as PleD does [29, 30]. Since cdG levels in C. crescentus are less than 0.3 μM [31], which is much lower than most protein levels, we use cdG to represent the PopA:2c-di-GMP binding species in this model (Fig 5). Moreover, the phosphorylation of CpdR is controlled by the kinase CckA, similarly with CtrA [14]. cdG binds to CckA to inhibit its kinase activity [32], which means cdG participates in the degradation and dephosphorylation of CtrA. CckA and cdG connect the master regulators network and ClpXP-based proteolysis system through CtrA (Fig 5).

Fig 5

Hierarchical diagram of protease complexes.

Solid lines with arrow denote metabolisms; solid lines with filled circles denote binding processes; dashed lines with arrow denote activation effects. Complex 1 decays the first class of substrates; Complex 2 degrades the second class of substrates; Complex 3 degrades the third class of substrates.

Table 1

Equations of replication and methylation, transcription, translation, and proteolysis.

	Equations of DNA
(1)	dInidt=k_s,Ini·([DnaA]Θ_DnaA)^4Ja,Ini4+([CtrA∼P]Θ_ctrA)^4+([DnaA]Θ_DnaA)^4·(1+1Ji,Ini4+(h_CoriΘ_Cori)^4)
(2)	dElongdt=k_elong·Elong^4Elong^4+Pelong4;(3)dZringdt=k_s,Zring
(4)	dh_Coridt=−k_m,Cori·[CcrM]^4Jm,Cori4+[CcrM]^4·h_Cori
(5)	dh_ccrMdt=−k_m,ccrM·[CcrM]^4Jm,ccrM4+[CcrM]^4·h_ccrM
(6)	dh_ctrAdt=−k_m,ctrA·[CcrM]^4Jm,ctrA4+[CcrM]^4·h_ctrA
	Equations of mRNAs
(7)	dI_ccrMdt=k_s,I_ccrM·([CtrA∼P]^2Ja,ccrM-CtrA2+[CtrA∼P]^2·Ji,ccrM-SciP2Ji,ccrM-SciP2+[SciP]^2)·h_ccrM−k_d,I_ccrM·I_ccrM
(8)	dccrMdt=k_s,ccrM·I_ccrM−k_d,ccrM·ccrM
(9)	ddnaAdt=k_s,dnaA·(Ji,dnaA-GcrA2Ji,dnaA-GcrA2+[GcrA]^2)·(2−h_Cori)−k_d,dnaA·dnaA
(10)	dgcrAdt=k_s,gcrA·([DnaA]^2Ja,gcrA-DnaA2+[DnaA]^2·Ji,gcrA-CtrA2Ji,gcrA-CtrA2+[CtrA∼P]^2)−k_d,gcrA·gcrA
(11)	dsciPdt=k_s,sciP·[CtrA∼P]^2Ja,sciP-CtrA2+[CtrA∼P]^2−k_d,sciP·sciP
(12)	dctrAdt=k_s1,ctrA·([GcrA]^2Ja,ctrA-GcrA2+[GcrA]^2·Ji,ctrA-CtrA4Ji,ctrA-CtrA4+[CtrA∼P]^4·Ji,ctrA-SciP4Ji,ctrA-SciP4+[SciP]^4)·h_ctrA+k_s2,ctrA·[CtrA∼P]^2Ja,ctrA-CtrA2+[CtrA∼P]^2−k_d,ctrA·ctrA
	Equations of regulatory proteins
(13)	d[CcrM]dt=k_s,CcrM·ccrM−k_d,CcrM·[CcrM]
(14)	d[DnaA]dt=k_s,DnaA·dnaA−k_d,DnaA·[DnaA]
(15)	d[GcrA]dt=k_s,GcrA·gcrA−k_d,GcrA·[GcrA]
(16)	d[SciP]dt=k_s,SciP·sciP−k_d,SciP·[SciP]
(17)	d[CtrA]dt=k_s,CtrA·ctrA−(k_d,CtrA+k_d,CtrA-ClpXP·[Complex3]^2Jd,CtrA-CplXP2+[Complex3]^2)·[CtrA]−k_phoCtrA·[CckA~P]·[CtrA]+k_dephoCtrA·[CtrA~P]
(18)	d[CtrA~P]dt=−(k_d,CtrA+k_d,CtrA-ClpXP·[Complex3]^2Jd,CtrA-ClpXP2+[Complex3]^2)·[CtrA~P]+k_phoCtrA·[CckA~P]·[CtrA]−k_dephoCtrA·[CtrA~P]
	Equations of protease complexes
(19)	d[CckA~P]dt=k_phoCckA·(CckAT−[CckA∼P])−k_dephoCckA·(1+α_cdG·[cdG])·[CckA~P]
(20)	d[Complex1]dt=k1+·[ClpXP]·[CpdR]−k1−·[Complex1]−k2+·[Complex1]·[RcdA]+k2−·[Complex2]
(21)	d[CpdR]dt=k_s,CpdR·[CtrA∼P]^2Ja,CpdR-CtrA2+[CtrA∼P]^2−k_d,CpdR·[CpdR]·[Complex1]J_d,CpdR+[Complex1]+k1−·[Complex1]−k1+·[ClpXP]·[CpdR]+k_dephos,CpdR·[CpdR∼P]−k_phos,CpdR·[CckA~P]·[CpdR]
(22)	d[CpdR∼P]dt=−k_d,CpdR·[CpdR~P]·[Complex1]J_d,CpdR+[Complex1]+k_phos,CpdR·[CckA~P]·[CpdR]−k_dephos,CpdR·[CpdR~P]
(23)	d[Complex2]dt=k2+·[Complex1]·[RcdA]−k2−·[Complex2]+k3−·[Complex3]−k3+·[c-di-GMP]^2·[Complex2]
(24)	d[RcdA]dt=k_s,RcdA·[CtrA~P]^2Ja,RcdA-CtrA2+[CtrA∼P]^2−k_d,RcdA·[RcdA]·[Complex1]J_d,RcdA+[Complex1]
(25)	d[Complex3]dt=k3+·[c-di-GMP]^2·[Complex2]−k3−·[Complex3]
(26)	d[PleD]dt=k_s,PleD·[CtrA~P]^2Ja,PleD-CtrA2+[CtrA~P]^2−k_d,PleD·[PleD]−k_phosPleD·[PleD]+k_dephoPleD·[PleC]·[PleD~P]
(27)	d[PleD~P]dt=k_phosPleD·[PleD]−k_dephoPleD·[PleC]·[PleD~P]
(28)	d[PdeA]dt=k_s,PdeA·[CtrA~P]^2Ja,PdeA-CtrA2+[CtrA~P]^2−k_d,PdeA·[PdeA]·[Complex1]J_d,PdeA+[Complex1]
(29)	d[cdG]dt=k_s,cdG·(1+α_PleD·[PleD])·Ji,cdG-cdG2Ji,cdG-cdG2+[cdG]^2−k_d,cdG·(1+α_PdeA·[PdeA])·[cdG]+k3−·[Complex3]−k3+·[cdG]·[Complex2]

Only phosphorylated form of PleD is active to catalyze the synthesis of cdG [32]. As the phosphorylation of PleD is controlled by more than three enzymes, including PleC, DivJ, CckN, and at least one unknown kinase [33, 34], it is complicated to thoroughly involve phosphorylation pathway of PleD. We initially assumed that phosphorylated PleD has a similar trend over cell cycle with total PleD and used total PleD as the synthetase of cdG; but cdG simulation in predivisional cell was super high, inconsistent with experiments, although both PleD and PdeA fit data well. We hypothesize PleD∼P is relatively low in predivisional cell due to the regulation of its main phosphatase PleC and kinase DivJ. To verify our hypothesis, we quantify western blots of DivJ and PleC over cell cycle using ImageJ [35, 36] (Fig 6). Experimental data indicates that DivJ almost does not change during the cell cycle and PleC is high in predivisional cell. Therefore, it is reasonable that PleD∼P decreases in predivisional cell because of high phosphatase activity of PleC. We fit PleC data points with trigonometric functions: 80.09 × sin(0.013t + 1.74) + 78.77 × sin(0.013t + 4.85) (Fig 6A). The function of PleC is then introduced into our model to represent the PleC level regulating the phosphorylation of PleD.

Fig 6

Quantification of experimental data of PleC from [35] and DivJ from [36].

(A) Black curve is fitted from experimental data by MATLAB. The fitting function is 80.09 × sin(0.013t + 1.74) + 78.77 × sin(0.013t + 4.85). (B) Experimental data of DivJ indicates DivJ levels sharply drop during the sw-to-st transition and then almost do not change.

Module 3: Chromosome replication and methylation

We build the module for DNA replication following the recognized principle in Li et al’s work [4], which consists of initiation, elongation, and termination phases. During the sw-to-st transition, C. crescentus requires high levels of DnaA and low levels of CtrA to initiate DNA replication [13]. As DNA synthesis proceeds, the fully methylated chromosome becomes hemimethylated due to the semiconservative replication. Replication will not be initiated again until CcrM re-methylates Cori once more [13]. Additionally, the master regulator genes ctrA, dnaA, and ccrM have CcrM-targeted sequence GANTC in their promoters (see Fig 2). Therefore, the methylation state of these genes are likely influenced by CcrM abundance and the progression of replication. Taken together, the initiation of DNA replication occurs when CtrA concentration is low, DnaA concentration is high, and Cori is fully methylated. Once initiated, DNA replication continues in a bidirectional manner along circular chromosomes and terminates in the late predivisional cell [37]. Finally, the newly replicated chromosomes are separated into two daughter cells with the Z-ring constriction.

We use variables Ini and Elong to model the initiation and elongation phase of DNA replication, respectively, where Elong was built by Li et al [4] (Table 1, Eqs. 1-2). Ini = Pelong signals the end of initiation and the beginning of the elongation phase, where Pelong = 0.05 [13]. The initial value of Elong is 0.1 (2 × P) because the chromosome replication of C. crescentus is bidirectional. DNA replication is terminated when Elong = 1 and we reset Elong = 0 once replication is terminated. h, indicating the probability of hemimethylation [4], is introduced in this study to describe the methylation influences on transcript rate (see Table 1, Eqs. 7-12). As the position of dnaA is very close to Cori (see Fig 2), ‘(2 − h)’ is used to represent the methylation effect of dnaA [4] and dnaA transcription rate reduces to half when it is hemimethylated [5]. I is introduced for a time delay. The chromosomes are separated with Z-ring constriction; however, the Z-ring event is not modeled in this study. Experiments indicate the S-phase period of Caulobacter is approximately 90 min. Here, we introduce a variable Zring to control the timing of Z-ring constriction and cell division. The increase rate of Zring is set as a particular constant to control the time for Zring rising from 0 to 1; and we use the time event Zring = 1 (Table 2) to signal the separation of chromosomes, where the count of chromosomes goes from 2 to 1 [6]. Throughout the execution of our simulation, several events representing cellular phenomena, including time points of replication initiation and chromosome segregation, can be triggered given particular conditions (summarized in Table 2).

Table 2

Event list.

Event description	Condition	Change(s)
DNA replication initiates	Ini = Pelong	Ini = 0, Elong = 0.05, DNA = 1.05, ks,Zring = 0.011, hCori = 1
replication fork passes ccrM locus	Elong = 0.2	hccrM = 1
replication fork passes ctrA locus	Elong = 0.375	hctrA = 1
DNA elongation terminates	Elong = 1	Elong = 0
Z-ring constriction	Zring = 1	Zring = 0, ks,Zring = 0, DNA = 1

Model derivation

Some proteins are not uniformly distributed in Caulobacter cells (see Fig 1). As we focus on temporal dynamics of regulators and their contributions to cell development, we ignore the non-uniform distributions and assume the whole cell is well-mixed at this stage. We use the law of mass action to describe the general synthesis/degradation and binding/unbinding processes, while protein effects–activation and inhibition–are described by Hill functions. To be more specific, is converted as where X represents protein, x is the mRNA of X, ks,X is the rate constant of synthesis, and kd,X indicates the rate constant of degradation.

is converted as A binds to B to produce C, where k+ and k− represent binding and unbinding rates, respectively. Activation and inhibition effects are described by Hill functions as follow: where Ha(X) indicates activation, and Hi(X) indicates inhibition. Variables n, J represent the corresponding Hill coefficient and the microscopic dissociation constant, respectively.

Model parameters

Experimental data

To compare our simulations with experimental observations from different publications, we first normalize experimental data to [0, 1] as follows: where x indicates the original data point; z is the scaled normalized value of experiments. Second, considering the relative abundance of different species in experiments [38], we set different targeted ranges for different species in the model. For example, the abundances of DnaA and CcrM are relatively low while those of CtrA and SciP are relatively high in experiments [38] and our simulations. For the figures in the Result section, the normalized experimental data are scaled to the range of our simulations to evaluate the temporal dynamics.

Parameter description

All 86 parameters used in this study are summarized in Table 3. Among them, seven parameters are obtained from previous experimental or modeling publications (see parameters marked with * in Table 3).

Table 3

Parameter values.

Parameters marked with * are obtained from publications.

Parameters of DNA

Rate constants, units = min-1

ks,Ini = 3.104e-4, kelong*=6.53e-3 [4], Pelong* = 0.05 [4],

km,Cori = 1.5637, km,ccrM = 2.2763, km,ctrA = 1.4645

Binding constants (dimensionless)

Ja,Ini = 1, Ji,Ini = 1.4565, Jm,Cori = 0.95, Jm,ccrM = 0.95, Jm,ctrA = 0.95

Scaling variables (dimensionless)

ΘCtrA = 6.0, ΘDnaA = 0.5, ΘCori = 0.308

Parameters of mRNAs

Rate constants, units = min−1

k_s,I_ccrM=0.1105, k_d,I_ccrM=0.0696, ks,ccrM = 0.2557, kd,ccrM = 0.1005,

ks,dnaA = 0.199, kd,dnaA = 0.0693, ks,gcrA = 5.4235, kd,gcrA = 0.7342,

ks1,ctrA = 1.0035, ks2,ctrA = 0.0937, kd,ctrA = 0.0983, ks,sciP = 0.583, kd,sciP = 0.0523,

Binding constants (dimensionless)

Ja,ccrM-CtrA = 5, Ji,ccrM-SciP = 6, Ji,dnaA-GcrA = 3, Ja,gcrA-DnaA = 1.25, Ji,gcrA-CtrA = 5,

Ja,ctrA-CtrA = 5, Ja,ctrA-GcrA = 3, Ji,ctrA-CtrA = 8, Ji,ctrA-SciP = 8, Ja,sciP-CtrA = 5

Parameters of master regulators

Rate constants, units = min−1

ks,DnaA = 0.0787, kd,DnaA*=0.07 [54], ks,GcrA = 0.032, kd,GcrA*=0.022 [4],

ks,CcrM = 0.0834, kd,CcrM*=0.07 [55], ks,SciP = 0.1294, kd,SciP = 0.0673,

ks,CtrA = 0.0404, kd,CtrA*=0.002 [4], kd,CtrA-ClpXP = 0.053

Phosphorylation constant, units = min−1

kphos,CtrA = 4.2919, kdephos,CtrA = 0.113, kphos,CckA = 1.027, kdephos,CckA = 0.9242

Binding constants (dimensionless)

Jd,CtrA-ClpXP = 4

Parameters of protease complexes

Rate constants, units = min−1

k1+=0.6072, k2+=1.4375, k3+=170.4913, k1−=3.3013, k2−=0.8164, k3−=2.3133,

kphos,PleD = 0.046, kdephos,PleD = 0.0414, ks,CpdR = 1.2227, kd,CpdR = 1.6152,

kpho,CpdR = 1.1239, kdepho,CpdR = 1.3854, ks,RcdA = 0.1642, kd,RcdA = 0.2323,

ks,cdG = 0.0099, kd,cdG = 0.9893,

ks,PleD = 0.0956, kd,PleD = 0.1314, ks,PdeA = 0.012, kd,PdeA = 0.5161

Binding constants (dimensionless)

Ja,CpdR-CtrA = 15, Jd,CpdR = 6, Ja,RcdA-CtrA = 15, Jd,RcdA = 2

Ja,PdeA-CtrA = 5, Jd,PdeA = 5, Ja,PleD-CtrA = 2.5, Ji,cdG-cdG*=0.2 [51]

Constants (dimensionless)

CckAT = 0.3, [ClpXP] = 1, αPdeA = 7, αPleD = 1500, αcdG = 10

The rest of the parameters are split into two groups: 1) 47 parameters (summarized in Table 4) that characterize major functionality of mRNAs and proteins, such as synthesis and degradation, are chosen for optimization; 2) the remaining 32 parameters are set with fixed values, including most dissociation constants.

Table 4

Parameter optimization with lower and upper bounds and starting point.

Parameter	[L,U]	Starting	Parameter	[L,U]	Starting
k m,Cori	[0.35, 5.6]	1.4	k m,ccrM	[0.35, 5.6]	1.4
k m,ctrA	[0.35, 5.6]	1.4	k_s,I_ccrM	[0.025 0.4]	0.1
k_d,I_ccrM	[0.016675, 0.2668]	0.0667	k s,ccrM	[0.064, 1.024]	0.256
k d,ccrM	[0.02, 0.32]	0.08	k s,dnaA	[0.0605,0.968]	0.242
k d,dnaA	[0.015, 0.24]	0.06	k s,gcrA	[1.4, 22.4]	5.6
k d,gcrA	[0.15, 2.4]	0.6	k s,sciP	[0.125, 2]	0.5
k d,sciP	[0.01, 0.16]	0.04	k s1,ctrA	[0.2475, 3.96]	0.99
k s2,ctrA	[0.0225, 0.36]	0.09	k d,ctrA	[0.02075, 0.332]	0.083
k s,DnaA	[0.01625, 0.26]	0.065	k s,GcrA	[0.007, 0.112]	0.028
k s,CcrM	[0.02125, 0.34]	0.085	k s,SciP	[0.0295, 0.472]	0.1183
k d,SciP	[0.015, 0.24]	0.06	k s,CtrA	[0.0108, 0.1728]	0.0432
k d,CtrA-ClpXP	[0.015, 0.24]	0.06	k1+	[0.15, 2.4]	0.6
k1−	[0.75, 12]	3	k s,CpdR	[0.175, 2.8]	0.7
k d,CpdR	[0.375, 6]	1.5	k dephos,CpdR	[0.25, 4]	1
k phos,CpdR	[0.25, 4]	1	k2+	[0.275, 4.4]	1.1
k2−	[0.25, 4]	1	k s,RcdA	[0.0375, 0.6]	0.15
k d,RcdA	[0.05, 0.8]	0.2	k3+	[35, 560]	140
k3−	[0.5, 8]	2	k s,cdG	[0.0025, 0.04]	0.01
k d,cdG	[0.25, 4]	1	k s,PleD	[0.025, 0.4]	0.1
k d,PleD	[0.0375, 0.6]	0.15	k phosPleD	[0.01, 0.16]	0.04
k dephosPleD	[0.01, 0.16]	0.04	k s,PdeA	[0.0025, 0.04]	0.01
k d,PdeA	[0.125, 2]	0.5	k phoCtrA	[1.25, 20]	5
k dephoCckA	[0.25, 4]	1	k dephoCtrA	[0.025, 0.4]	0.1
k phoCckA	[0.25, 4]	1

Multiobjective optimization

Let , p = 47 be the vector of parameters to be estimated in the caulobactor cell cycle model. For this optimization problem, we focus on two aspects: the quantitative difference between experimental data and simulated results, and the cell cycle time, both for wild type cells. The reasoning behind the two objectives is that the experimental data have inconsistent concentration levels between the beginning (t = 0 min) and ending (t = 150 min) states of the cell cycle, whereas our model must be consistent to ensure stable cell cycle regulation. This is also validated in our initial optimization test using a single objective function, where we observe minimizing the difference in species concentration results in high deviation in cell cycle time, and vice versa. Due to this conflict, we cannot use the common scalarization scheme to sum up the two objectives using weights, i.e., F(χ) = w1f1(χ) + w2f2(χ). Our parameter optimization is therefore defined as a multiobjective optimization problem (MOP). The two objective functions are: where x denotes the simulated concentration of species i at time j, y denotes the experimental data of species i at time j, and T is the simulated cell cycle time. Here, we have the experimental data for n = 15 species (see Table 5) and the number of data points m varies for different species. Note that we only use available observations of C. crescentus wild type (WT) cells for parameter fitting. The mutant cases of C. cell are used as our model validation. The optimization problem to be solved is where [L, U] is a search box in for model parameters. See the lower and upper bounds of parameters in Table 4.

Table 5

Sources for experimental data used to evaluate our models.

Species	Data source	Species	Data source
ccrM	[12]	CcrM	[47, 48]
dnaA	[12]	DnaA	[7, 47]
gcrA	[12]	GcrA	[2, 7]
sciP	[12]	SciP	[2]
ctrA	[12]	CtrA	[47, 48]
CpdR	[25]	RcdA	[44]
PleD	[45]	PdeA	[45]
cdG	[31]

We apply two MOP algorithms to our optimization problem for comparison: one is the widely used nondominated sorting genetic algorithm (NSGA-II) [39]; the other is the more recent VTMOP [40] based on VTdirect [41] and QNSTOP [42, 43] that uses response surface and trust region methodologies, and an adaptive weighting scheme. Initial values in Table 6 are the levels of corresponding variables used as the beginning state of each simulated cell cycle. Fig 7 shows the combined Pareto front from both methods after multiple runs with different optimization settings. Observe that the Pareto front is a nonconvex curve, showing that the multiobjective optimization problem is very difficult. The best parameter estimates are found by VTMOP and listed in Table 3, with f1 = 1.57 and f2 = 0.02. The sensitivity of parameters is 18% for experimental data fitting (f1) and 72% for cell cycle time (f2) if we perturb the parameters of three Pareto points (marked as black circle in Fig 7) by 10%. Note that the sensitivity of the second objective is large when f2 is very close to zero, thus points near zero are not included in the analysis. The root mean square error (goodness of fit) is RMS(f1) ≈ 1.51, RMS(f2) ≈ 0.84 for the three selected Pareto points.

Table 6

Initial values of model variables.

DNA variables	Initial values	Master regulator vars.	Initial values
Ini	0.0383	CcrM	0.435
Elong	0	DnaA	2.638
DNA	1	GcrA	3.841
Count	1	SciP	12.485
Zring	0	CtrA	1.973
hCori	0	CtrA∼P	3.960
hccrM	0
hctrA	0
mRNA variables	Initial values	Protease complex vars.	Initial values
ccrM	0.173	Complex 1	0.211
dnaA	3.154	CpdR	1.045
gcrA	4.525	CpdR∼P	0.042
sciP	6.335	Complex 2	0.187
ctrA	0.658	RcdA	0.789
		Complex 3	3.407
		CckAP	0.042
		cdG	0.511
		PleD	0.526
		PleD∼P	0.663
		PdeA	0.228

Fig 7

Pareto front returned by NSGA-II and VTMOP.

f1(χ) and f1(χ) are the two objective function values. The red point is selected as the best parameter set for our model.

Results

Integrating the hierarchical proteolysis (Fig 5) into the master regulatory network (Fig 3), we propose a model to capture the temporal dynamics of cell cycle regulators and glean insights about bacterial protein proteolysis systems. Non-uniform distributions of molecules in space are ignored at this stage.

Our model accurately describes gene transcription patterns and temporal dynamics of key regulators during the replication cycle of Caulobacter wild type cells

Chromosome replication and methylation

We follow the DNA replication process as the rationale to formulate a set of ordinary differential equations (ODEs) modeling initiation, elongation, and termination of DNA replication as well as methylation states, as shown in Table 1 (Eqs. 1-6). The initiation of DNA replication requires a fully methylated state (both strands methylated, h* = 0), while semiconservative replication creates two hemimethylated copies of genes. As such, the variables h* in our model spike when the corresponding gene is being replicated (Fig 8). Later in the cell cycle, the hemimethylated copies (h* = 1) are re-methylated by CcrM, returning to the fully methylated state. Therefore, h* then plunge as the newly created, hemimethylated copies become fully methylated by CcrM. The CcrM-dependent methylation in the control system ensures DNA replication initiates once per cell cycle.

Fig 8

(A) Simulated chromosome/DNA replication and elongation process; (B) The probability of loci (Cori, ccrM, ctrA) being hemimethylated in a single cell cycle.

The proteolysis of CtrA is controlled by hierarchical protease complexes

In addition to replication and transcription, we investigate the proteolysis regulation of CtrA, the essential component of cell cycle control system, and explore the contribution of the conserved proteolysis module. Based on the hierarchical diagram of protease complexes (Fig 5), we use ODEs to simulate the temporal dynamics of three classes of protease complexes (Eqs. 20-29 of Table 1). Since there is no experimental data of protease complexes, we evaluate our simulations using western blots of CpdR, RcdA, PleD, PdeA, and cdG [14, 44, 45] (see Table 5), where numerical values are extracted by ImageJ or GetData, shown as the red circles in Fig 9. Those proteins are essential components of ClpXP-dependent proteolysis system.

Fig 9

(A-E) The dynamics of total CpdR, RcdA, cdG, PdeA, total PleD, and PleD∼P in simulation with the corresponding experimental data. Experimental data of CpdR is from Iniesta et al. [25], RcdA is from McGrath et al. [44], cdG is from Abel et al. [31], and PdeA as well as total PleD are from Abel et al. [45].

Our simulated CpdR, PleD and PdeA match well the experimental dynamics (see Fig 9). The general trend of modeled RcdA and cdG agrees with experiments, whereas cdG peaks a little bit late compared with experimental data. The discrepancy may derive from other regulatory enzymes of cdG or PleD which are not involved in our current model. As most proteins involved in protease complexes are modeled reasonably, We use the hierarchical model to simulate the cyclic proteolysis of CtrA (Eq. 17 of Table 1). In addition to degradation regulation, the hierarchical model influences the phosphorylation of CtrA via cdG and CckA, while phosphorylated CtrA in turn impacts the expression of components involved in degradation module, including cpdR, rcdA, and pleD.

Temporal dynamics of mRNA and master regulators

We convert the regulatory network diagram in Fig 3 into ODEs shown in Table 1 (Eqs. 7-18) to simulate the temporal dynamics of five master regulators and their mRNA. The proposed hierarchical protease complexes are applied to simulate the cyclic degradation of CtrA. Fig 10A–10E and 10F–10J exhibit the comparisons between simulations (black curves) in our model and experimental data (red circles and blue triangles) for mRNA (dnaA, gcrA, ctrA, ccrM, and sciP) and protein (DnaA, GcrA, CtrA, CcrM, and SciP) levels, respectively. In general, our simulations fit the experimental observations well. As we capture the genetic information flowing from mRNA to proteins, the protein concentration curves generally resemble the corresponding mRNA abundance curves. dnaA transcription is reduced by hemi-methylation state, which in part explains the dip in our simulation of dnaA during DNA replication (t ∈ [30, 110] min in Fig 10A). Additionally, the expression of dnaA is activated by CtrA [5] and inhibited by GcrA [46]. Thus, the high levels of GcrA and low levels of CtrA during sw-to-stalk transition reinforce the decrease of dnaA expression (Fig 10A, 10G and 10H). When the replication fork passes ccrM and ctrA right before and after 50 min, h and h are switched from 0 to 1 (Fig 8), which explains the increase of CcrM (ccrM) and CtrA (ctrA) at the corresponding time (Fig 10B, 10D, 10I and 10H). Meanwhile, the high levels of activator GcrA and low levels of inhibitor SciP amplify the increase of ctrA. DnaA and CtrA collaborate to regulate the initiation of DNA replication: 1) during sw-to-st transition, initiator DnaA is high and suppressor CtrA is low, allowing the cell to initiate replication; 2) during DNA replication, DnaA keeps low and CtrA is high, avoiding another initiation of replication in the same cycle. Under the combined functions of DnaA and CtrA, the transcription of of gcrA increases in the beginning and decreases in the predivisional stage, which agrees with the observation of gcrA transcription (see Fig 10B). sciP expression is activated by CtrA, which is observed in our simulation as well (Fig 10E). Fig 11A shows the maximum levels of our simulated master regulators, in which the relative scales agree with experiments [38]. We summarize the simulated and observed abundance of five master regulators in a single cell cycle in a bar chart (Fig 11B), where our simulation shows similar translation patterns with experiments. Even though the experimental data comes from a variety of sources and experimental techniques, visual inspection suggests fair agreement between the timing of master regulator abundance in simulation and experimental data.

Fig 10

(A-E) Experimental mRNA concentration of dnaA, gcrA, ctrA, ccrM, sciP (curves) with corresponding simulated data (red circles, from Schrader et al. [12]), and (F-J) simulated protein concentration of DnaA, GcrA, total CtrA (CtrA∼P), CcrM, SciP (curves) with experimental data (circles or triangles) over a single cell cycle. For the sources of experimental data, DnaA data is from Shen et al. [47] and Collier et al. [7]; GcrA data is from Collier et al. [7] and Tan et al. [2]; CtrA data and CcrM data are both from Reisenauer et al. [48] and Shen et al. [47]; and SciP data is from Tan et al. [2].

Fig 11

(A) Relative maximum concentrations of master regulators across one swarmer cell cycle. (B) Abundance of five master regulators (CcrM, CtrA, DnaA, GcrA, and SciP) from simulated results and experimental data. Horizontal bars represent the time periods of protein abundance across the swarmer cell cycle. Blue bars indicate the time frame where simulated protein levels are above the mid-range concentrations and red bars are the corresponding experimental data from [24].

One objection worth noting is that some of our simulations deviate from the experimental data at the beginning or the end of the cell cycle. For example in Fig 10C, at the end of the cell cycle, the expression level of ctrA is considerably lower in our simulation than the experimental data suggest. Additionally, this type of discrepancy can be witnessed in the simulated GcrA, where the simulated level is lower than experimental observations after t = 100 min (Fig 10G). This disagreement stems from a limitation that the simulated endpoint has to be equal to the starting point (t = 0 min), because we do not model the asymmetrical heritage of two distinct daughter cells after progenies are completely separated. Although there are several mismatches between the simulation and experiments, our model fits most data points and captures key trends during cell cycle, such as the dynamics of regulators during sw-to-st transition. Our model exhibits that key regulators interact with each other through transcription, degradation, and phosphorylation regulations to determine the timing of cell differentiation and reproduction.

Hierarchical protease complexes contribute to the timed cell cycle progression

Our modeled hierarchical cyclic proteolysis module performs well in the simulation of CtrA. Here, we explore the contribution of this module for cell development. We replace the cyclic proteolysis by a constant for CtrA, CpdR, and RcdA, separately, setting Jd,CtrA−ClpXP, Jd,CpdR, or Jd,RcdA as 0. In the simulation of Jd,CtrA−ClpXP = 0, where the degradation rate of CtrA is constant, the system still oscillates during cell cycles whereas the amplitude of CtrA and SciP shows noteworthy reduces. The cycle time increases, resulting in delays of master regulators in simulation, including CtrA and CcrM (Fig 12A). With a constant degradation of CpdR or RcdA, simulations show severe defects, especially for the dynamics of CtrA. The oscillation of CtrA almost disappears and methylation states are abnormal under these conditions (Fig 12B and 12C). In summary, the cyclic proteolysis deriving from the hierarchical protease complexes shows significant impacts on the system. We further replace all cyclic complexes with constants, setting Jd,CtrA−ClpXP, Jd,CpdR, and Jd,RcdA as 0 simultaneously. The corresponding simulation is similar with the Jd,CtrA−ClpXP = 0 mutant, which shows delayed cell cycles and reduced amplitudes of several species (Fig 12D). Simulations of these cyclic proteolysis mutants suggest the cyclic proteolysis of CtrA is key to regulate the entire system, because both deletion (Jd,CtrA−ClpXP = 0) and changes (Jd,CpdR = 0 and Jd,RcdA = 0) of CtrA cyclic degradation would screw up the dynamics pattern of both master regulators and their mRNA. Moreover, the system is more sensitive without the cyclic proteolysis module. We increased the degradation rate of CtrA to 5-fold for system with and without the cyclic proteolysis module; wild type system still has an acceptable cell cycle while cyclic proteolysis mutant systems have severe deficiencies. Taken together, our model suggests the hierarchical cyclic proteolysis module contributes the timed cell cycle and robustness of the Caulobacter system.

Fig 12

(A-D) Simulated results of mutating the cyclic proteolysis of CtrA, CpdR, or/and RcdA. (A). Jd,CtrA−ClpXP = 0 indicates the cyclic proteolysis of CtrA is replaced by a constant. (B). Jd,CpdR = 0 indicates the cyclic proteolysis of CpdR is replaced by a constant. (C). Jd,RcdA = 0 indicates the cyclic proteolysis of RcdA is replaced by a constant. (D). The cyclic proteolysis of CtrA, CpdR, and RcdA are all mutated as constant degradation.

Our model captures major phenotypes of mutant strains

To further test the validity of our model, we use the same equations and initial values to simulate seven different mutant strains (Fig 13). Among these mutant strains, cell cycle of ΔdnaA, where dnaA is knocked out (ks, = 0), is arrested. The other six mutant strains are all viable. Our mutant simulations correctly capture the viability of these seven mutant strains. To be more specific:

Fig 13

(A-G) Simulated results of mutant strains:ΔccrM, ΔgcrA, ΔdnaA, ctrAΔ3Ω, cdG0, ΔpdeA, and ΔpleD. In knock out mutant simulations, we set ks,i = 0, where i indicates corresponding species including ccrM, gcrA, dnaA, cdG, PdeA, and PleD. In the simulation of ctrAΔ3Ω, the cyclic proteolysis rate of CtrA is reduced to 10%.

ΔccrM: ccrM is verified to be dispensable for cell viability [10]. The doubling time is about 162±9 min, longer than that for WT. Our simulated ΔccrM (ks, = 0) has a 164 min cycle time, which fits the experimental observation well (Fig 13A). In our simulation, all h can not be returned to 0 because there is no CcrM re-methylating the chromosome. Additionally, experiments have suggested the cell cycle is also regulated by CcrM independent with GANTC motif. This study does not include the GANTC motif independent influence of CcrM, so the simulation of ΔccrM only shows the potential of deleting methylation of GANTC motif by CcrM.

ΔgcrA: In gcrA knocked out strain, the doubling time is 40% longer than for the WT [10]. Our simulated gcrA mutant (ks, = 0) has approximately 10 min longer cell cycle compared with the WT simulation (Fig 13B), which less than then the experimental observation. The gap is likely derived from the forced modeling of Z-ring constriction process, which is not explicitly modeled in this study.

ctrAΔΩ3: ctrAΔΩ3 contains modifications to the C-terminal amino acids of ctrA, resulting in a non-proteolizable CtrA allele [49]. Here, we decreases kd,CtrA-ClpXP to 10% of WT in simulation. In Fig 13D, The average CtrA levels increase in simulation with less fluctuation because of the non-proteolizable CtrA allele. Our simulation suggests the proteolysis of CtrA is important for its cell cycle-dependent regulation. h in simulated ctrAΔΩ3 can not decrease to 0, suggesting the levels of CcrM in ctrAΔΩ3 are not sufficient to completely re-methylate chromosome, while the lower levels of CcrM is caused by higher levels of the inhibitor CtrA.

cdG related mutant strains

cdG0 mutant strain has been verified to be viable, although it shows morphology defects [32]. In Fig 13E, our simulation of cdG0 strain (ks,cdG = 0) is viable and shows a horizontal shift which may result in morphology defects. The CtrA levels increase with less fluctuation which is caused by the deletion of cdG. pleD knocked out mutant (ks,PleD = 0) results in a lower cdG levels (Fig 13G), which shows a similar phenotype with the simulation of cdG0. pdeA mutant increases cdG levels in simulation (Fig 13F, ks,PdeA = 0). Both ΔpleD and ΔpdeA are viable in simulation, consistent with observations [32, 45]. Oscillations exist but shifts little bit in the simulations of these three cdG regulated mutants, as shown in Fig 13E–13G.

Discussion

The five major regulators–DnaA, GcrA, CcrM, CtrA, and SciP–work in tandem to drive the cell cycle progression of C. crescentus. Here, we investigated the interactions among master regulators to study the underlying mechanisms of DNA replication, methylation, transcription, and proteolysis of cyclic regulators. We applied the central dogma of molecular biology to simulate the temporal dynamics of mRNAs and proteins. Furthermore, we mathematically built a hierarchical model to simulate protease complexes and apply this model to CtrA degradation. Two MOP approaches (NSGA-II and VTMOP) have been applied to estimate parameters in this complicated system.

In C. crescentus, the protease ClpXP primed by one assistant adaptor recruits additional adaptors in sequence [27]. The hierarchical adaptor assembly determines the time and location of the proteolysis of hierarchical substrates. Our hierarchical model correctly captures the key dynamics of CpdR, PleD, and PdeA; it shows fair agreement with the trend of RcdA and cdG. Additionally, the protease model performs well in modeling the proteolysis of CtrA. Deleting the hierarchical protease module causes defects of cell cycle development and protein oscillations. Considering the fast formation of the protease Complex 3 (ClpXP bound with CpdR, RcdA, PopA, and cdG), we test quasi-steady-state assumption (QSSA) for Complex 3. QSSA shows similar simulated results in both wild type cells and mutant cases, suggesting QSSA might be a good approach in reducing model complexity of biological systems. As a wide range of proteins is degraded by ClpXP protease complex, our model provides a good quantitative tool to analyze the proteolysis of these proteins in C. crescentus, such as TacA and ShkA. As most components of the hierarchic protease complexes are conserved in bacterial species, our model has the potential for a wide range of applications. Moreover, cdG is a significant component of Complex 3 and participates in several essential pathways of cell cycle regulation in C. crescentus. For example, cdG binds to CckA and ShkA to induce phosphatase and kinase activity, respectively. While CckA controls the phosphorylation/dephosphorylation of several proteins, such as CtrA and CpdR, ShkA:cdG regulates the phosphorylation of TacA, which downregulates the stalked pole muramidase homolog SpmX and the stalk length regulator StaR [50]. Additionally, cdG has been verified to participate in the stress response, contributing to the survival of Caulobacter in oligotrophic environments [51]. Due to the importance of cdG, our protease complex model is potentially a valuable tool for understanding the regulatory network of C. crescentus.

With the advances in experimental technologies, mRNA and protein abundance of master regulators have been monitored and measured throughout the cell cycle. However, there is a limited comparison between experiments and simulations. Our results align very well with the experimental data. Satisfactory simulation results of our model, as indicated by visual inspection, suggest that the proposed regulatory network appropriately characterizes the Caulobacter cell cycle progression. This study also suggests the cell cycle dependent proteolysis of CtrA is significant for the cell cycle regulations and robustness. Our model can capture major features of seven mutant strains, which has the potential to predict phenotypes of nonviable mutant strains and functions of involved proteins. As most molecules involved in our model (CtrA, CcrM, GcrA, DnaA, etc.) are conserved among proteobacteria [18, 52, 53], this framework could be applied to the study of other proteobacteria. Last but not least, this work is a successful application of multiobjective optimization problem, showing that MOP is a promising approach for handling conflicting objectives in biological modeling.