Computational Modeling of the Chlamydial Developmental Cycle Reveals a Potential Role for Asymmetric Division

ABSTRACT Chlamydia trachomatis is an obligate intracellular bacterium that progresses through an essential multicell form developmental cycle. Infection of the host is initiated by the elementary body (EB). Once in the host, the EB cell differentiates into the noninfectious, but replication-competent, reticulate body, or RB. After multiple rounds of replication, RBs undergo secondary differentiation, eventually producing newly infectious EBs. Here, we generated paired cell-type promoter reporter constructs and determined the kinetics of the activities of the euo, hctA, and hctB promoters. The paired constructs revealed that the developmental cycle produces at least three phenotypically distinct cell types, the RB (euoprom+), intermediate body (IB; hctAprom+), and EB (hctBprom+). The kinetic data from the three dual-promoter constructs were used to generate two computational agent-based models to reproduce the chlamydial developmental cycle. Both models simulated EB germination, RB amplification, IB formation, and EB production but differed in the mechanism that generated the IB. The direct conversion and the asymmetric production models predicted different behaviors for the RB population, which were experimentally testable. In agreement with the asymmetric production model, RBs acted as stem cells after the initial amplification stage, producing one IB and self-renewing after every division. We also demonstrated that IBs are a transient cell population, maturing directly into EBs after formation without the need for cell division. The culmination of these results suggests that the developmental cycle can be described by a four-stage model, EB germination, RB amplification/maturation, IB production, and EB formation. IMPORTANCE Chlamydia trachomatis is an obligate intracellular bacterial pathogen responsible for both ocular and sexually transmitted infections. All Chlamydiae are reliant on a complex developmental cycle, consisting of both infectious and noninfectious cell forms. The EB cell form initiates infection, whereas the RB cell replicates. The infectious cycle requires both cell types, as RB replication increases the cell population while EB formation disseminates the infection to new hosts. The mechanisms of RB-to-EB development are largely unknown. Here, we developed unique dual promoter reporters and used live-cell imaging and confocal microscopy to visualize the cycle at the single-cell and kinetic levels. These data were used to develop and test two agent-based models, simulating either direct conversion of RBs to EBs or production of EBs via asymmetric RB division. Our results suggest that RBs mature into a stem cell-like population producing intermediate cell forms through asymmetric division, followed by maturation of the intermediate cell type into the infectious EB. Ultimately, a more complete mechanistic understanding of the developmental cycle will lead to novel therapeutics targeting cell type development to eliminate chlamydial dissemination.

2019, the CDC reported 1.8 million C. trachomatis infections within the United States alone, with the most recent reports indicating that rates increased by 10.0% in women and 32.1% in men from 2015 to 2019 (1,2). This increase in infection rates has been reported across all racial/ethnic groups and affects all age groups (2).
Chlamydial growth and development have classically been characterized as a biphasic cycle, consisting of two primary cell forms, the elementary and reticulate body (3). These cell forms maintain a division of labor throughout the infectious cycle and are essential for chlamydial proliferation. The elementary body (EB) is the infectious cell form and initiates host cell invasion by pathogen-mediated endocytosis (4). The EB cell form is nonreplicative, and the chromosome is tightly compacted by nucleoid-associated proteins (5,6). Upon entry into the host, the EB undergoes large transcriptional and phenotypic changes, maturing into the reticulate body (RB) in a process that takes up to 12 h for serovar L2 (7,8). The RB is replication competent but noninfectious and must redifferentiate back into the EB to disseminate the infection to new host cells (3,9).
Electron micrographs have also shown the presence of a transitory cell form, termed the intermediate body (IB). IBs are present beginning between 20 and 24 h postinfection (hpi) for serovar L2 and are characterized by a semicondensed nucleoid similar to the EB, but they are significantly larger, which gives them a target-like appearance (7,10). Due to the presence of the IB and its appearance as a transitory form, it is currently hypothesized that a subset of RBs undergoes large morphological changes to convert directly into EBs.
We previously reported the development of a live-cell reporter system to follow the chlamydial cycle in real time at the single-inclusion level. A suite of different gene promoters was designed to drive the expression of fluorescent proteins in order to follow RB growth and EB development. These kinetic data suggested that the promoters fell into three temporal categories exemplified by the activity of the euo, hctA, and hctB promoters (10).
In this study, paired promoter reporter constructs were developed to determine the temporal and spatial relationships between the activities of the euo, hctA, and hctB promoters at the kinetic and single-cell levels. Based on the expression kinetics of these reporters and their kinetic relationships to each other, computational agent-based models were created to best represent the developmental cycle. Two models were developed to explain the data, an RB-to-IB direct conversion model and an asymmetric division/production model. The outputs of simulations from these models were compared to experimental data to determine which mechanism was best supported. Our model and data suggest a novel RB amplification/maturation step where the RB initially divides symmetrically to produce two RBs and increase RB numbers, followed by maturation of the RB to an asymmetrically dividing cell that produces one IB while regenerating the RB. Our data also support the direct maturation of the IB cell type into the EB without the need for cell division.

RESULTS
Development of dual fluorescent cell reporter constructs to determine celltype gene expression during the developmental cycle. We previously reported the use of promoter reporter constructs that were expressed at different stages of the chlamydial developmental cycle (10). Fluorescence from the euoprom-Clover construct was first detected at ;14 h postinfection (hpi) and demonstrated a short exponential phase, followed by an expression plateau at ;24 hpi (10). Imaging of the hctAprom-Clover construct revealed that the Clover signal could first be detected at ;18 hpi and continued to increase until cell lysis (10). Additionally, we noted that hctBprom-Clover expression was initiated at ;22 hpi, ;4 h later than hctAprom, and also increased linearly until cell lysis (10).
Here, we extended these studies by combining these promoters in pairs on the same plasmid. In addition, the LVA degradation tag was used to target fluorescence proteins for degradation to increase temporal resolution. To visualize expression kinetics and expression relationships between these cell-type promoter reporters at the single-inclusion and single-cell level, we constructed three dual fluorescent developmental gene expression reporter strains, hctAprom-euoprom, hctBprom-hctAprom, and hctBprom-euoprom. For the hctAprom-euoprom (AMELVA) construct, the euo promoter was used to drive the expression of the green fluorescent protein variant, mNeonGreen (mNG) (11) fused in-frame to the LVA protein degradation tag, which reduced the fluorescent half-life to ;30 min (12), and the hctA promoter was used to drive the expression of the red fluorescent protein mKate2 (13) (see Fig. S1A in the supplemental material). The hctBprom-euoprom (BMELVA) dual-reporter construct was created by replacing hctAprom in AMELVA with hctBprom to drive mKate2 expression (Fig. S1B). Last, to create the hctBprom-hctAprom (BMALVA) dual reporter, the hctA promoter was used to drive mNG(LVA), and the hctB promoter was used to drive mKate2 expression (Fig. S1C). Each construct was transformed into Chlamydia trachomatis L2, creating the L2-AMELVA, L2-BMELVA, and L2-BMALVA reporter strains (Table 1) and compared to our previously reported dual-promoter strain L2-hctBprom_mKate2-euoprom_Clover (L2-BMEC) (10). To determine the relationships between the promoter reporters, we investigated the expression of each paired reporter at the single chlamydial cell level by confocal microscopy (Fig. 1A), and the single inclusion kinetic level by livecell imaging (Fig. 1B). For confocal imaging, host cells were infected with each strain for 24 h followed by fixation. All L2-BMEC chlamydial cells showed some level of green fluorescence, suggesting that all cells expressed from the euoprom-Clover construct at some point during the developmental cycle (Fig. 1A, L2-BMEC). However, confocal microscopy of the new strains using the short-half-life mNG(LVA) fluorescent reporter revealed that the three promoters were active exclusively from one another (Fig. 1A). The euoprom-mNG(LVA) signal from the L2-AMELVA strain was present in large RB-sized cells, while hctAprom-mKate2 was active in a subset of large and small cells but in a population distinct from euoprom 1 cells (Fig. 1A, L2-AMELVA). The euoprom-mNG(LVA) signal from L2-BMELVA was also present in large cells only, while the hctBprom-mKate2 signal was present in small cells and in an entirely distinct population (Fig. 1A, BMELVA). In the L2-BMALVAinfected cells, hctAprom-mNG(LVA) expression was visible in large and small cells, while the hctBprom-mKate2 signal was detected in small cells. Again, the two promoter reporters were active in distinct, nonoverlapping populations (Fig. 1A, L2-BMALVA). Occasionally, we saw paired hctAprom-expressing cells, but we do not know the origin or fate of these cells, as these are static images of a dynamic environment.
We quantified the expression of the paired fluorescent reporters at the single-cell level in four individual inclusions for each strain. The TrackMate plugin in Fiji (14) was used to identify green fluorescent protein (GFP)-and red fluorescent protein (RFP)-expressing cells and record the fluorescent signal for each reporter from entire confocal z stacks encompassing each inclusion. The paired fluorescent signals for each cell were plotted, the RFP on the y axis and the GFP signal on the x axis (Fig. 1A). For the L2-BMEC strain expressing the stable GFP protein, Clover, driven by the euo promoter, both green and red cells had Clover signal, demonstrating that both cell types expressed Clover from the euo promoter at some point during the developmental cycle of that cell. In contrast, the red-and green-expressing Chlamydia cells of the L2-BMELVA strain, which expresses a short-half-life fluorescent protein from the euo promoter, are distinct populations with minimal overlap. The same was seen for the L2-AMELVA and L2-BMALVA infections; the red-expressing cells had minimal green signal (AMELVA), and the green-expressing cells had minimal red signal (BMALVA). Together, these data demonstrate that the three promoters are active in distinct subpopulations of chlamydial cells.
In addition to the confocal images, the kinetics of each promoter was determined  Modeling of the Chlamydial Developmental Cycle mSystems at the single-inclusion level using live-cell microscopy. Host cells were infected with each strain and imaged for both GFP and RFP from 10 hpi until 50 hpi at 30-min intervals. Euoprom-mNG(LVA) signal from L2-AMELVA was first detected at ;14 hpi and increased exponentially until ;26 hpi, after which time the signal reached a plateau that was maintained for the duration of the infection (Fig. 1B, L2-AMELVA). The hctAprom(mKate2) signal was first detectable at ;18 hpi, with an exponential increase in expression until 28 hpi, followed by a linear increase until the end of the experiment at 50 hpi (Fig. 1B, L2-AMELVA). Like L2-AMELVA, the euoprom-mNG(LVA) signal for L2-BMELVA followed the same kinetics, with an early exponential increase followed by a signal plateau (Fig. 1B, L2-BMELVA). The hctBprom signal in these inclusions became detectable at ;24 hpi and increased exponentially until ;34 hpi. After this brief exponential phase, hctBprom(mKate2) signal increased at a linear rate until the end of the experiment (Fig. 1B, L2-BMELVA). Livecell kinetics of L2-BMALVA showed that hctAprom-mNG(LVA) activity initiated around 18 hpi; however, expression approached steady-state kinetics and did not accumulate (Fig. 1B, L2-BMALVA). HctBprom(mKate2) demonstrated the same kinetics as the L2-BMELVA strain (Fig. 1B). The short half-life of the LVA-tagged mNG allowed for spatial and kinetic resolution of three cell types. These results demonstrate that chlamydial cells follow a temporal gene expression program resulting in three distinct cell forms. The RB cell activity expresses from the euoprom, while the EB cell expresses from the hctBprom. The observation that hctAprom was temporarily active after euoprom and before hctBprom and the observation that hctAprom is active in a distinct cell population from euoprom-positive (euoprom 1 ) cells and hctBprom 1 cells suggest that hctAprom is active in the IB cell population. Overall, these data suggest that the developmental cycle can be represented by three phenotypically distinct cell types, RB cells (euoprom 1 ), IB cells (hctAprom 1 ), and EB cells (hctBprom 1 ). Live-cell imaging also revealed the interrelated kinetics of the activity of these promoters. The kinetic data from each individual inclusion of both L2-AMELVA ( Fig. 2A) and L2-BMELVA ( Fig. 2B) revealed that there was significant heterogeneity in the maximal expression level of the euoprom plateau. Interestingly, the variation in euoprom expression from the paired expression constructs correlated with the rates of hctAprom (L2-AMELVA) and hctBprom (L2-BMELVA) signal accumulation, i.e., inclusions exhibiting high euoprom plateaus had steeper slopes for hctAprom and hctBprom signal accumulation, while lower euoprom signal correlated with a lower rate of hctAprom and hctBprom signal accumulation ( Fig. 2A and B). When the hctAprom and hctBprom signals were normalized to the euoprom signal plateau for each inclusion, the variation in the slopes of hctAprom(mKate2) and hctBprom(mKate2) accumulation was dramatically reduced ( Fig. 2A and B). This variability in the maximum expression levels of euoprom suggested that each inclusion contained differing numbers of RB (euoprom 1 ) cells during the plateau phase, which, in turn, led to the various IB (hctAprom) and EB (hctBprom) accumulation rates.
To test whether the differing plateau signal of euoprom expression was due to differing RB numbers, cells were infected with L2-BMELVA at a multiplicity of infection (MOI) of ;0.1 and fixed and stained with DAPI (49,6-diamidino-2-phenylindole) every 2 h from 14 to 48 hpi. Infected cells were imaged by three-dimensional (3D) confocal microscopy, and cell type quantification was carried out using an automated cell counting workflow using the open-source software Fiji and the TrackMate plugin to count individual Chlamydia based on fluorescent reporter intensity (14). These experiments revealed that RB (euoprom 1 ) cells increased in number from 14 hpi to 26 hpi, reaching an average of ;30/inclusion (Fig. 2C). After this time point, the average number of RBs was unchanged. However, it was clear that  the maximum number of RBs in each inclusion varied significantly, with some inclusions having as few as one euoprom 1 cell, while others had as many as 59 RBs/inclusion during this plateau phase (26 to 48 hpi) (Fig. 2C). A similar time course was carried out for cells infected with L2-BMALVA, with similar results. IBs (hctAprom 1 ) cells increased in number from 17 hpi until reaching a maximum at 32 hpi, after which the average of IBs/inclusion remained steady (Fig. 2D). Again, similar to the euoprom 1 cells, the number of hctAprom 1 cells was significantly different on a per-inclusion basis, with some inclusions having as few as a single hctAprom 1 cell, while others had as many as 65 during the plateau phase (Fig. 2D). These data suggest that the chlamydial developmental cycle produces significant heterogeneity between inclusions, but despite this heterogeneity, each inclusion produces similar kinetic relationships between cell types.
Modeling the chlamydial developmental cycle. Dissecting the mechanisms that control the developmental cycle in Chlamydia has, in part, been inhibited by the reliance on population-based studies. Using the individual inclusion kinetic data and individual cell expression data generated from L2-AMELVA, L2-BMELVA, and L2-BMALVA, we divided the cycle into discrete steps, EB germination, RB amplification, IB formation, and EB maturation. To explore potential mechanisms involved in these steps, two agent-based models (ABMs) were developed to describe the developmental cycle using the Python-based bacterial growth simulation platform, Cellmodeller (15). Euo, HctA, and HctB protein expression was simulated in each cell type over time for both models (Text S1). Additionally, germination was set at 10 hpi, as the first replication event has been previously reported at this time point (10,16), and it agrees well with the initiation of euo promoter expression ( Fig. 1 and 2) (10). The two models differed in the mechanism controlling RB amplification and IB formation. The asymmetric production model used an RB maturation mechanism where early RBs (designated the RB R ) replicate to produce two identical RB R daughter cells, resulting in RB R number amplification. This step is followed by an RB R maturation phase where the RB R matures over time into a cell form that undergoes asymmetric cell division (designated the RB E ), producing one RB E daughter cell and one IB daughter cell. The direct conversion model used a stochastic direct conversion mechanism where early RBs replicate, resulting in RB number amplification followed by an increase in the chance that an RB transitions into the IB state. This stochastic chance of conversion increases over time before reaching a maintenance state where RB amplification and conversion rates are matched (Text S1). Both models reproduced the developmental kinetics that were observed using L2-AMELVA, L2-BMELVA, and L2-BMALVA ( Fig. 1 and 2; Fig. 3A). However, the stochastic model needed to be constrained to match the experimental data. When conversion became greater than RB replication, RB numbers dropped to extinction. Conversely, when the RB replication rate remained higher than the conversion rate, RBs quickly outnumbered EBs (Fig. 3B).
Although both models could produce similar kinetics at the population level, simulations of individual inclusions demonstrated large kinetic differences. For the asymmetric production model, simulated inclusions with high RB numbers had corresponding high EB production rates, while inclusions with low RB numbers had corresponding low EB production rates. The EB production rate for each inclusion, regardless of RB numbers, was linear, while the RB population numbers remained unchanged over time (Fig. 3C, asymmetric production). The kinetics for the direct conversion model had similar trends; however, there were obvious runs of over-and underamplification/conversion for both RB and EB numbers in the individual inclusion simulations, demonstrating that the stochastic mechanism can reliably reproduce the observed data only on the population level. These simulations suggest that the asymmetric production model is a better match for the observed developmental kinetics.  Steady state versus stem cell population. The asymmetric production model predicts that after amplification, the RBs act as a stem cell population, producing IBs while selfrenewing after every division (Movie S1). In contrast, the direct conversion model predicts that RBs convert directly into IBs and are subsequently replaced by a dividing RB population, producing a steady state of RBs (Movie S2). To determine which of these phenotypes the RB is exhibiting, live-cell imaging at high resolution (40Â objective) of cells infected with L2-BMELVA was used to follow RB behavior. Cells infected with L2-BMELVA were imaged for euoprom and hctBprom expression every 15 min starting at 24 hpi until 60 hpi. We imaged from 24 hpi until 60 hpi, as this covers the end of the RB amplification stage until the end of the cycle. Imaging revealed that the number of RBs (euoprom 1 ) per inclusion remained roughly the same throughout the experiment ( Fig. 4A and Movie S3). The RBs in each inclusion were easily tracked from one time point to the next and did not disappear when new RBs appeared. We also observed two binary divisions (RB amplification events) occurring between 26 hpi and 31 hpi (Fig. 4A, cells 3a and b and 5a and b, and Movie S3). These newly formed RBs also remained trackable within the inclusion for the remainder of the experiment. These observations are consistent with stem cell-like behavior of the RBs after amplification.
Both models predicted that the IB cell type is maintained at steady state, i.e., IBs mature into EBs, and new IBs replace the maturing cells. To investigate the dynamics of the IB  Together, these data strongly support the asymmetric production model that predicts that RB numbers are amplified between 12 hpi and 28 hpi followed by a stem cell-like behavior, producing one IB at division. In contrast, the IBs behave as a steady-state population where, once formed, the IBs are converted directly into EBs.
The role of cell division in the RB population. Simulations of the two competing ABMs predicted different developmental outcomes if cell division was inhibited. The asymmetric production model predicted that RBs produce IBs only at division. Therefore, an immediate block in the formation of new IBs would occur, but RB numbers would remain unaffected if replication was inhibited at 30 hpi (Fig. 5A). In contrast, the direct conversion model predicted that RB numbers would drop over time if cell division was inhibited at 30 hpi, eventually disappearing as the RBs converted into IBs but were not replaced by further RB replication (Fig. 5B). To test these predictions experimentally, two cell replication inhibitors were used, penicillin (Pen) and ciprofloxacin (Cip). Chlamydia does not contain a peptidoglycan cell wall and instead uses peptidoglycan in septum formation; therefore, Pen treatment of Chlamydia inhibits cell division (17). Ciprofloxacin prevents bacterial DNA replication by inhibiting topoisomerases and DNA gyrase (18). Cells were infected with L2-BMELVA, and .20 inclusions per treatment were imaged using a 60Â 1.4 numerical aperture (NA) objective. Individual inclusions were imaged at multiple Z planes to visualize and quantify all the RBs in the inclusion. Images were collected from the live cultures at antibiotic treatment (30 hpi) and 10 h later (40 hpi) (Fig. 5C and D). The number of euoprom 1 cells (RBs) was quantified on a per-inclusion basis. The same inclusions were quantified at each time point. Consistent with the confocal time-series experiment, there was a large variation in the number of RBs in individual inclusions, ranging from a single RB to greater than 50 RBs (Fig. 5C). However, the number of RBs per inclusion remained essentially constant between the 30-hpi and 40-hpi time points (Fig. 5C). The mean ratio of RBs at 30 hpi and 40 hpi per inclusion was 1.27 6 0.28 for untreated, 1.18 6 0.24 for Pen treated, and 1.51 6 0.35 for Cip treated (Fig. 5D). There was an increase in RB numbers in each inclusion when treated with Cip, but never more than double. We speculate that although Cip treatment inhibited initiation of DNA replication, it did not significantly impact the continuation of DNA replication and allowed cell division to finish in RBs that had already initiated DNA synthesis. This is in agreement with our published data demonstrating that RBs are undergoing continuous DNA replication, as indicated by a replication index of 1.5 (19). This index is a ratio of sequencing coverage near the origin of replication versus coverage near the terminus and represents the growth rate of the population. An RB replication index of 1.5 demonstrates the presence of partially replicated chromosomes (19).
To confirm that chlamydial DNA replication was inhibited by Cip, digital droplet PCR (ddPCR) was performed on L2-BMELVA-infected samples treated with either Cip, Pen, or mock at 30 hpi. Host monolayers were harvested every 4 h from 26 to 54 hpi. As previously reported, genome copy number continued to increase in the Pen-treated samples (10,20). There was, however, a large reduction in genome copy accumulation after Cip treatment compared to the mock and Pen-treated samples (Fig. S2B).
The role of cell division in IB population and EB production. Simulations of the two competing ABMs also predicted different developmental outcomes in the IB population if cell division was inhibited. The asymmetric production model produces IBs only at division; therefore, blocking cell division at 30 hpi predicted an immediate block in the formation of new IBs, leading to an immediate halt in the increase in IB gene expression (Fig. 6A, asymmetric production). In contrast, after replication was inhibited, the direct conversion model predicted that IB gene expression would continue to increase over a 12-h period as the RBs converted into IBs, at which point IB gene expression would halt, as IBs could not be replenished by RB cell division (Fig. 6A, direct conversion).
To test these predictions experimentally, cell replication was inhibited with Pen, Cip, or genetically by overexpression of FtsI. Penicillin-binding proteins (PBPs) are a suite of enzymes involved in peptidoglycan synthesis and cell division (21). Both PBP2 and PBP3 are (D,D)-transpeptidases that mediate the formation of peptidoglycan. In Escherichia coli, PBP2 is required for cell shape and stability, whereas PBP3 is involved in cell division (22). When overexpressed in E. coli, PBP2 and PBP3 have been shown to induce lysis in dividing cells (23,24). However, as Chlamydia does not contain a peptidoglycan cell wall, PBP2 and PBP3 have both been co-opted for use in septation (25). Therefore, we overexpressed FtsI (PBP3) to target dividing chlamydial cells. The open reading frame (ORF) of chlamydial ftsI was cloned into our translational expression system and tagged with a C-terminal 3ÂFLAG epitope (12). The euoprom-mNG(LVA)_hctBprom-mKate2 promoter-reporter cassette was cloned into the E-ftsI-3XFLAG vector and transformed into Chlamydia trachomatis L2, creating the strain L2-E-ftsI-BMELVA. Ectopic expression of full-length FtsI-FLAG in Chlamydia was verified by Western blotting (Fig. S2A). Ectopic FtsI expression resulted in a block in chlamydial replication, as demonstrated by inhibition of chromosome replication measured using digital droplet PCR (Fig. S2B). Additionally, we examined the impact of FtsI overexpression on the RB population using immunofluorescence confocal imaging of host cells infected with L2-E-ftsI-BMELVA (Fig. S2C). Infected cells were treated with vehicle only or induced for expression at 20 hpi. Inclusions from the 20-hpi and 30-hpi uninduced samples were filled with bacteria and contained multiple bright euoprom-mNG(LVA) 1 cells, indicating the presence of live RBs capable of active transcription and translation (Fig. S2C). Conversely, inclusions from the 30-hpi FtsI-induced samples were relatively empty compared to the 20-hpi and 30-hpi uninduced samples. A subset of cells in the induced sample stained positively for FLAG (i.e., indicating FtsI induction); these cells were misshapen and contained significantly less euoprom-mNG(LVA) fluorescence than uninduced cells, indicating a lack of active expression and suggesting RB lysis had occurred (Fig. S2C). This construct contained a FLAG tag at the C terminus of FtsI, potentially altering its function. We have, however, used C-terminal FLAG tags in other chlamydial protein fusion constructs and not observed any impact on the developmental cycle (12,19,26).
An additional dual promoter reporter strain was created using hctAprom-mEos3.2 paired with hctBprom-mKate2 (L2-BMAMEO). The mEos3.2 GFP variant accumulates over time, as it does not have the LVA degradation tag and will therefore be present in all cells that have expressed from the hctAprom (IBs and EBs). This dual-color promoter reporter cassette was also cloned into the E-ftsI-3XFLAG plasmid and transformed into Chlamydia, creating L2-E-ftsI-BMAMEO. To differentiate between asymmetric production and direct conversion, cells were infected with L2-BMAMEO or L2-E-ftsI-BMAMEO and imaged from 10 hpi until 50 hpi (Fig. 6A, experimental data). Pen and Cip were added at 30 hpi to the L2-BMAMEO cultures, and FtsI was induced by adding theophylline (Tph) at 30 hpi to the L2-E-ftsI-BMAMEO-infected cells to inhibit chlamydial replication. Live-cell imaging of the hctAprom-mEos3.2 produced kinetics consistent with the asymmetric production model simulations. The increase of IB production signal from hctAprom was almost immediately inhibited by all three treatments (Fig. 6A, experimental data). As we have previously documented, the hctAprom reactivated in the Pen-treated aberrant cells after about an ;8-h delay (10).
The two models also predicted differences in EB production when cell division was inhibited. The asymmetric production model predicted that EB production would halt ;6 h after all treatments indicated (Fig. 6B, asymmetric production). The stochastic direct conversion model simulations predicted that inhibition of cell division would result in a slowing followed by a halt in EB production ;18 h posttreatment (Fig. 6B, direct conversion). To directly test these predictions, cells were infected with L2-BMAMEO or L2-E-ftsI-BMAMEO. Pen or Cip were added to L2-BMAMEO cultures, and Tph was added to the L2-E-ftsI-BMAMEO-infected cells at 30 hpi. EBs were collected every 4 h from 10 to 50 hpi and used to infect fresh host cells for EB quantification. Our data show that EB production was inhibited ;8 h post-cell division inhibition with all three inhibitors (Fig. 6B, experimental data). These data, taken together, support a division-dependent asymmetric IB production model and not stochastic direct conversion of an RB to an IB.
We also asked whether inhibiting cell division through FtsI ectopic expression affected IB-to-EB formation. Cells were infected with L2-E-ftsI-BMELVA and induced for FtsI expression at 20 hpi. The infected cells were fixed and stained with DAPI, and confocal images were taken for DAPI, GFP, and RFP signals at 20 and 30 hpi. The inclusions in the L2-E-ftsI-BMELVAinfected cells at 20 hpi contained euoprom 1 cells (RBs) and DAPI-only stained cells (IBs) and little to no hctBprom 1 cells (EBs) (Fig. S3D; Fig. S4A). The uninduced L2-E-ftsI-BMELVA inclusions at 30 hpi contained euoprom 1 cells (RBs), DAPI-only positive cells (IBs), and a significant number of hctBprom 1 cells (EBs) (Fig. S4B). However, the Tph-induced L2-E-ftsI-BMELVA inclusions at 30 hpi had misshapen large chlamydial cells with just a trace of euoprom 1 signal (RBs), fewer DAPI-positive cells (IBs), and a significant number of hctBprom 1 cells (EBs) (Fig. 7B; Fig. S3D). Taken together, these data support the hypothesis that the IB matures directly into the EB cell form without undergoing cell division.
IBs convert directly to EBs. To verify that the IB directly matures into the EB without cell division (hctAprom 1 to hctAprom 1 or hctBprom 1 ), cells were infected with L2-BMAMEO to visualize the hctAprom 1 and hctBprom 1 cells, and cell division was inhibited with Cip at 18 hpi. Infected cells were then fixed and stained with DAPI at 22 hpi and 34 hpi. The inclusions were imaged for DNA (DAPI), GFP (hctAprom), and RFP (hctBprom) (Fig. 8A and B). The expression levels of GFP and RFP were determined for .500 individual chlamydial cells in three representative inclusions by identifying cells using the DAPI channel and measuring expression levels using the TrackMate plugin in Fiji ( Fig. 8A and B). Confocal microscopy of untreated cells at 22 hpi revealed that there were a number of DAPI-only (RBs) and hctAprom 1 cells (IBs) with very few hctAprom 1 and hctBprom 1 cells (IB!EBs) (Fig. 8A, UNT). At 34 hpi, there were again populations of DAPI-only (RBs) and hctAprom 1 cells (IBs); however, many of the hctAprom 1 cells were also positive for hctBprom (IB!EB) (Fig. 8B, UNT). Similar trends were seen when cell division was inhibited by the addition of Cip at 18 hpi. The chlamydial population at 22 hpi consisted of primarily DAPI-only (RBs) and hctAprom 1 cells (IBs) (Fig. 8A, CIP), while many of the hctAprom 1 cells were also hctBprom 1 (IB!EB) by 34 hpi (Fig. 8B, CIP). This increase in the number of double-positive cells in cell division-inhibited Chlamydia suggests that the hctAprom 1 cells were activating hctBprom over time.
We also asked whether the IB could directly become an EB after replication was inhibited by ectopic expression of FtsI. Cells were infected with L2-E-ftsI-BMAMEO, and FtsI expression was induced at 18 hpi. At 22 hpi, there were a number of DAPI-only (RBs) and hctAprom 1 cells (IBs) in the induced population. Anti-FLAG staining revealed that many of the DAPI-only cells were expressing FtsI-FLAG (Fig. 8A, FtsI). The inclusions at this time point contained very few hctBprom 1 cells (EBs). At 34 hpi, there was a significant increase in the number of hctAprom/hctBprom double-positive chlamydial cells (IB!EB), suggesting  the hctAprom 1 cells activated hctBprom without cell division (Fig. 8B, FtsI). The uninduced samples showed a similar pattern. These data together demonstrate that IBs (hctAprom 1 ) mature directly into EBs (hctBprom 1 ) and that cell division is not required for this progression.

DISCUSSION
The developmental cycle of Chlamydia is central to its ability to cause disease. The cycle produces phenotypically distinct cell types with functional specificity. This includes the RB cell, which replicates, leading to organism proliferation, and the EB cell type, which mediates entry into new cells, disseminating the infection to new hosts. The developmental cycle has conventionally been broken down into two stages, RB replication and EB conversion. While the broad strokes of the overall cycle are described, the molecular details that regulate this process are poorly understood. Our data suggest that the current understanding of the cycle is an oversimplification. We propose that the developmental cycle can be divided into multiple stages, including EB germination, RB amplification and maturation, IB production, and EB formation.
Our data show that the chlamydial developmental cycle produces at least three different phenotypic cell types, RB (euoprom 1 ), IB (hctAprom 1 ), and EB (hctBprom 1 ). The data also demonstrate that RB proliferation initiates at ;10 to 12 hpi and reaches a plateau by 24 to 26 hpi, after which RB numbers remain virtually unchanged until the end of the cycle. Remarkably, at the individual inclusion level, the number of RBs at plateau was highly variable, with some inclusions containing just a few RBs, while others had as many as 60.
To understand the potential mechanisms that underlie these observations, we developed computational agent-based models (ABMs). Two mechanistic models were developed that could recapitulate our experimental data. Both models estimate the germination time based on published time to first RB replication (16) and time to euoprom expression (this study). The models differ primarily in the IB production mechanism, which we hypothesize is the committed step to EB formation. Regulatory mechanisms, such as RB access to or competition for inclusion membrane contact (26), reduction in RB size (7), and responses to changes in nutrient availability (27), have been proposed to explain the regulation of EB formation. Although the triggering signal for differentiation differs in these models, all propose a stochastic direct RB-to-EB conversion mechanism. Therefore, the direct conversion model utilized a direct conversion mechanism to control RB amplification and IB production. Our second ABM, the asymmetric production model, used an RB maturation/asymmetric division mechanism to control RB amplification and IB production. In this model, we proposed two RB subtypes, the RB R , which, upon replication, produces two identical RB R daughter cells, and an RB E that, upon cell division, produces one IB and one RB E . The RB R subtype matures into the RB E subtype during the first 24 h of infection. For both models, the IB cell exits the cell cycle and matures directly into the EB. Although both models could reproduce the measured kinetics of the developmental cycle at the population level, only the maturation/asymmetric model could reliably reproduce the observed kinetics at the single-inclusion level. Even after constraining replication and conversion to avoid RB extinction or overpopulation, the stochastic direct conversion model resulted in runs of over-and under-IB production that were not readily evident in the measured data.
The asymmetric production model predicted that the mature RB cell (RB E ) acts as a stem cell, producing one IB upon division, while the other daughter cell remains an RB E . The direct conversion model predicted that RB numbers are at steady state where cell division creates new RBs, while other RBs convert into IBs. Using live-cell microscopy, we determined that individual RBs were consistently present and trackable from time point to time point over a 30-h time period, strongly suggesting the RB population is acting like stem cells and not a steady-state population. Conversely, both models predicted that the IB population is at steady state, with IBs forming and transitioning into EBs, while new IBs replace those that became EBs. The behavior of individual IBs revealed by live-cell imaging was consistent with a steady-state population, with IBs disappearing and new ones appearing over time.
The models also predicted that inhibition of cell replication would have differing effects on the RB, IB, and EB populations. The stochastic direct conversion model predicted that RB numbers would decline after cell division inhibition, as RBs that directly converted to IBs were not replaced by new RBs through cell division. In contrast, the asymmetric production model predicted that inhibition of cell division would have no impact on the RB population. We found that treatment with two different cell replication inhibitors (Pen and Cip) resulted in no decrease in the RB population over a 10-h period, again strongly supporting an asymmetric division model.
The kinetics of the IB and EB population was also predicted to differ in the two models after replication was inhibited. The direct conversion model predicted that inhibition of cell division would lead to a protracted inhibition of new IB formation, resulting in decreasing IB production over ;12 h, while EB formation would also slowly decline until completely inhibited at ;18 h after treatment. The asymmetric production model predicted that inhibiting cell division would inhibit new IB production nearly immediately and inhibit EB formation after the ;8-h IB-to-EB maturation time (10). Live-cell microscopy demonstrated that the kinetics of cell division inhibited Chlamydia, which clearly supported the asymmetric production model over the direct conversion model, as IB production halted nearly instantly, and EB production halted after ;8 h.
Both models predicted that IBs mature directly into EBs. This prediction is supported by the data showing that the IB cell type progressed to hctBprom expression after inhibiting cell division with Cip or Pen or through ectopic expression of FtsI. Additionally, we demonstrated that IBs mature directly from hctAprom single-positive to hctAprom/hctBprom double-positive cells, even after inhibiting cell division with Cip or overexpression of FtsI.
Cell size has been proposed as a mechanism controlling stochastic RB-to-IB conversion. The average chlamydial cell size becomes smaller and more heterogeneous over time (7). Our asymmetric production model also predicts that cell size will behave similarly, decreasing in size and increasing in heterogeneity over time. In our model, early in the developmental cycle, the RB R cell progresses through the cell cycle and divides to produce two RB R s, which immediately reenter the cell cycle, gaining mass before the next division. As the cycle progresses, RB R s mature into RB E s. The RB E continues in the cell cycle and, upon division, produces one RB E and one IB. After division, the RB E remains in the cell cycle, increasing in size; however, the IB does not reenter the cell cycle, does not begin to add mass, and is therefore smaller than the RB R/E cells. The IB then proceeds to reduce in volume as it matures to the EB cell over an ;8-h time frame (10). This process would produce IB cells of various sizes with a trend toward smaller and smaller cells in the inclusion, as early IBs increase in number compared to RBs.
Overall, the data support a developmental cycle that includes an asymmetrically dividing RB population. Asymmetry has been documented for both the EB and RB cells (28)(29)(30)(31). In addition, it has been shown that the division plane can form asymmetrically during RB division (28,32). Asymmetric cell division is a common mechanism to generate phenotypically distinct cell populations in bacteria. Many of these systems have evolved to create two-cell populations; one cell acts as the stem cell, while the other cell disseminates the bacterial colony to new environments. Both Caulobacter crescentus and some members of Chlamydia's nearest phylogenetic neighbors, the Planctomycetes, undergo a division cycle that includes a surface-attached mother cell that produces a planktonic swarmer cell upon division in order to extend the population to new ecological niches (33,34). The swarmer cell in the case of C. crescentus is nonreplicating and is out of the cell cycle (35). Our data support a similar role for the EB cell, as the EB disseminates the infection, does not replicate, and is out of the cell cycle (19).
The mother/swarmer cell developmental model fits our data for IB and EB production but does not explain RB expansion. We have modeled this process (RB expansion to RB and IB asymmetric division) as maturation over time from an RB that produces two RBs, thereby increasing the dividing cell population (RB R ) to an RB mother/stem cell (RB E ) that produces an IB through asymmetric division. Currently, the mechanism for RB R to RB E maturation is unknown but could be influenced by several factors, such as EB age at infection, differences in nutrient acquisition, and the effects of inclusion membrane interaction, or through a yetto-be-described stochastic maturation mechanism. Additionally, our data show that this expansion step produces significant heterogeneity in RB E numbers in individual inclusions. This differential amplification is potentially a novel evolutionary adaptation to balance RB maturation and early EB production, i.e., late maturation leads to late EB production but a high EB production rate, whereas early maturation leads to early EB production with a low EB production rate. We currently have no expression markers that show differences between the two theoretical RB subtypes, the RB R and the RB E . Our data and subsequent model suggest that this is the best explanation for our data and will be a continuing area of investigation.
The asymmetric model suggests that asymmetry is generated during cell division after a maturation process. We have not, however, detected differential expression between forming daughter cells during division using our currently available dual-color promoter reporters. The two forming daughter cells appear to express similar levels of fluorescence from the euoprom-mNG(LVA) reporter and also appear to be similar sizes. This has been reported by other studies that showed that although there is evidence that the division plane initiates asymmetrically, at cytokinesis, the two daughter cells appear to be of similar size (8,18,24,25). Our data suggest that differential gene regulation of the hctA promoter in the two daughter cells happens after the cells have separated, thereby increasing the difficulty in identifying the daughter cell pairs.
Our data support a four-stage model, EB germination, RB amplification and maturation, IB production, and EB formation (Fig. 9). A key aspect of this model is the stochastic amplification and maturation of RBs (RB R ) from an expanding population to a stem cell RB population (RB E ) that produces IBs through asymmetric replication reminiscent of stalk/swarmer cell dynamics (Fig. 9). This is followed by an EB formation stage where the IB undergoes a dramatic phenotypic change that includes the expression of the nucleoid-associated protein, HctA, and ultimately ends with the expression of a second nucleoid-associated protein, HctB. This EB formation stage takes ;8 to 10 h, ultimately resulting in the infectious EB (10) (Fig. 9). Our model suggests that the cycle starts with RB amplification and quickly converts to a mixed environment of RB R s producing RBs and RB E s producing IBs. The IB phenotype takes ;8 to 10 h to complete maturation into the infectious EB. Our data suggest that replication is significantly faster than IB-to-EB maturation; therefore, IBs accumulate and are a significant cell type in the inclusion for most of the cycle, outnumbering RBs as soon as ;18 hpi. It is currently unclear what role the IB cell type has in the cell biology of chlamydial infection.
Temporal gene expression during chlamydial infection has been extensively reported at the population level. These data have historically categorized developmental gene expression into early-, mid-, and late-expressed genes (8,36). However, these studies were all performed at a population level and normalized accordingly. Our model and supporting data suggest that this gene expression cascade is explained by a shift in cell types within the inclusion. Early in the cycle (;10 to 16 hpi), the inclusion is dominated by RB cells; by 18 hpi, there is a mix of RBs and IBs; and by 24 hpi, early IBs have matured into EBs, and IBs begin to outnumber the RBs. Additionally, from 30 hpi and on, EBs quickly accumulate, becoming the dominant cell form, overwhelming the gene expression signal of the other cell types. Our data suggest that to fully understand the developmental cycle, the kinetically dynamic mixed-cell environment needs to be taken into consideration.
Our experiments and models have focused on Chlamydia serovar L2 living in close to ideal growth conditions (10) However, the developmental cycle is shared by all of the Chlamydiaceae family. Our studies focused on serovar L2, as it has become established as a model organism for understanding the basic life cycle of these organisms, as it lends itself to genetic manipulation and infects cultured cells efficiently. Our studies are likely to result in a further understanding of the chlamydial developmental cycle across the family. Eventually, we would like to determine if what is reported here describes all the Chlamydia species and chlamydial-like organisms, but further studies are needed. Additionally, it is likely that Chlamydia reacts and adapts the cycle to nutrient and other stresses. The creation of an ABM that models the mechanisms of the developmental cycle under ideal conditions will provide a tool to visualize and understand the convolved data obtained from nutrientlimiting, pharmacological, genetic, and molecular experiments. Ultimately, a better mechanistic understanding of the developmental cycle will lead to novel therapeutics targeting development, as breaking the cycle will eliminate dissemination and chlamydial disease.
was utilized for quantification of chlamydial genomic copies. A 2Â ddPCR supermix for probes-no dUTP kit (Bio-Rad) and a custom copN-specific primer/probe set was used for DNA detection (Table S1).
Live-cell microscopy. Monolayers were seeded on a multiwell glass-bottom plate and infected with Ctr-L2 EBs. Infections were grown in an Oko-Touch CO 2 -heated stage incubator. Fluorescence images were acquired via epifluorescence microscopy using a Nikon Eclipse TE300 inverted microscope with a ScopeLED lamp at 470 nm and 595 nm and BrightLine band-pass filters at 514/30 nm and 590/20 nm. We used 20Â/0.4 NA dry, 40Â/0.6 NA dry, and 60Â/1.40 NA oil objective lenses. Differential interference contrast (DIC) was used to autofocus images. Image acquisition was performed using an Andor Zyla sCMOS camera in conjugation with mManager software (40). Images were taken in 30-min intervals unless otherwise stated. Imaging ranged from 10 to 60 h after Ctr-L2 infection, depending on the experiment. Multiple fields were imaged for each treatment, and the fluorescent intensity of individual inclusions was monitored using the TrackMate plugin in Fiji (14). Inclusion fluorescent intensities were averaged and graphed in Python as previously described (39).
Confocal microscopy. Cos-7 cells were seeded onto glass coverslips and infected with the appropriate Ctr-L2 strains. Samples were fixed at the designated times in 2% paraformaldehyde in phosphatebuffered saline (PBS) at room temperature overnight. Samples were then washed with filtered PBS and stained with DAPI to visualize DNA and monoclonal anti-FLAG M2 antibody (Sigma, Thermo Scientific) with Alexa 647 anti-mouse secondary antibody to visualize FtsI-3XFLAG expression. Coverslips were mounted onto a microscope slide using Mowiol (100 mg/mL 150 Mowiol 4-88, 25% glycerol, 0.1 M Tris, pH 8.5). Images were acquired using a Nikon spinning disk confocal inverted microscope with a 100Â/ 1.45 NA oil objective with a laser lamp at 405 nm, 490 nm, 568 nm, and 660 nm. Image acquisition was performed using an Andor Ixon electron-multiplying charge-coupled-device (EMCCD) camera and the Nikon Elements software. Collected images were taken as z stacks at 0.2-mm intervals, capturing the entirety of each inclusion. Multiple inclusions were imaged for each treatment and time point, and quantification of individual cells was performed using TrackMate. Chlamydial cell numbers were then analyzed in custom Python notebooks. Representative confocal micrographs are maximal-intensity projections of 3D data sets, and brightness and contrast were adjusted equally for comparisons.
Western blot analysis. To ensure that ectopically expressed FtsI was full size, infected monolayers were lysed in reducing lane marker sample buffer, and protein lysates were separated on a 6% SDS-PAGE gels and transferred to a nitrocellulose membrane for Western blot analysis of the FLAG-tagged protein. The membrane was blocked with PBS with 0.1% Tween 20 (PBS-T) and 5% nonfat milk prior to incubating in monoclonal anti-FLAG M2 antibody (1:40,000; Sigma, Thermo Scientific) overnight at 4°C followed by goat-anti mouse IgG-horseradish peroxidase (HRP) secondary antibody (Invitrogen) at room temperature for 2 h. The membrane was developed with the SuperSignal West Dura luminol and peroxide solution (Thermo Scientific) and imaged using an Amersham Imager 600.
Computational modeling. Modeling was done using the CellModeller platform (15). The model description scripts and model data analysis scripts are described in Text S1 and available on GitHub (https://github.com/SGrasshopper/Chlamydial-developmental-cycle).
Data availability. All data, bacterial strains, and methodologies are available upon request.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. VIDEO S1, MOV file, 6.8 MB.