Recent technological advances have made feasible studies of biological systems at the single-cell level ^[1]–^[4]. However, our current understanding of single-cell biochemistry and physiology has been largely inferred from averaged population measurements that often mask individual cell dynamics and lead to a distorted picture of cell behavior. Such cell population data can be difficult to reconcile with single-cell models, such as those that attempt to describe cell-cycle-dependent gene expression kinetics ^[5]–^[7]. In particular, mathematical models of single cells that rely on population data for constraints on biochemical parameters may arrive at incorrect conclusions.

Among the properties hidden by population averaging is cell-to-cell variability, such as that found in gene expression and protein production ^[8]–^[11]. We refer to the natural variation found between cells at the same position in their cell cycles as synchronous variability. A population experiment in which synchronous variability is the only source of variability can at most yield the average of the observable of interest (e.g., gene expression levels). However, in addition to the inherent synchronous variability, typical time-series experiments on cells contain a significant asynchronous variability: even if cells have been physically or chemically synchronized, individual cells within a synchronized population exist at variable points in their respective cell cycles. As a result, the extraction of ‘true’ temporal data from such populations is difficult, since contributions from cells in different stages of the cell cycle are averaged.

From a mathematical perspective, population asynchrony may be modeled as a kernel function that maps the average of an observable in the absence of asynchronous variability to the value measured at the population level. Population asynchrony has been modeled in yeast as both a time-dependent ^[12],^[13] and time-independent ^[14] source of variability. With an accurate asynchrony model, extracting the average of an observable becomes an inverse problem for which established regularization methods can be used. These computational methods can effectively remove from population data artifacts that are due solely to asynchrony, or uncover features that are masked by population averaging ^[13]–^[15]. The resulting data is thus better suited for comparison with single-cell models and parameter estimation.

Population asynchrony characterization is most easily done with a synchronizable system such as the dimorphic bacterium Caulobacter crescentus. Caulobacter begins its cycle as a motile ‘swarmer’ (SW) cell and differentiates to a non-motile ‘stalked’ (ST) cell just prior to the initiation of DNA replication. The SW stage is thus analogous to the G1 phase of the eukaryotic cell cycle, and the ST stage is analogous to the S and G2 phases ^[16]. At the SW-to-ST transition, the flagellum is released and a narrow cylindrical extension of the cell envelope (the ‘stalk’) is grown in its place. A new flagellar assembly is constructed at the pole opposite the stalk as the cell cycle progresses, and on cell division, a new motile, chemotactic SW cell is spawned. The remaining ST cell immediately commences another round of DNA replication and division while the SW cell begins the full cell cycle (Fig. 1). Centrifugation of a mixed culture of Caulobacter in Ludox or Percoll separates SW cells from all other cell types, so that nearly pure cultures of SW cells can be easily obtained ^[17],^[18]. However, even a perfectly pure culture of SW cells includes a mixture of new and old SW cells, and variance in the cell cycle times of individual cells within this synchronized population leads to a further increase in the heterogeneity of the population as time-series experiments progress. Additional heterogeneity is introduced following cell division, as each dividing cell results in both a SW and ST cell. Thus, even a perfectly synchronized population develops a significant and time-dependent population asynchrony.

We propose a simple model for the time-dependent distribution of Caulobacter cell types in a population during synchronized growth. Our model accurately matches observed distributions of synchronized Caulobacter cells during a time-series experiment, and may be extended to any organism for which the synchrony state can be characterized—particularly those that undergo asymmetric division. We then combine a generalization of deconvolution with our Caulobacter distribution model to extract the “single-cell”-like synchronous average of gene expression profiles from published cell cycle microarray data. The resulting expression profiles more accurately predict the cell-cycle position and size of gene expression peaks, display new features not evident in the original microarray data set, and demonstrate robustness to uncertainty in model parameters. This represents a new advance in the study of cell-cycle dependent gene expression in Caulobacter. The deconvolution method presented herein can be generally applied to characterize time-dependent processes in a variety of biological model systems.

While population-level experimental techniques typically allow for high-throughput and fast data collection, they are unable to capture many of the details present at the level of single cells. This is an unavoidable consequence of population averaging; population-based data are in fact transforms of organism- and condition-specific population asynchrony kernels with single-cell data. Thus, an assumption of equivalence of population and single-cell data is an assumption of a non-physical delta function integral kernel. Recognizing this, cell distribution models have been proposed with the aim of extracting more detailed information from biological time-series data. Perhaps the simplest improvement on the delta function model is a fixed kernel such as a Gaussian. Further improvements have been made by allowing for a Gaussian kernel whose width increases with time (e.g., ^[13]). However, a normal distribution of this kind is not sufficient to describe the complex cell-phase distribution of organisms that undergo asymmetric division, and attempts to deconvolve single-cell expression for such organisms will lead to unreliable results. As a result, we have developed an intuitive mathematical model of the cell-type (or, alternatively, the cell-phase) distribution of asymmetrically-dividing cells as a function of time following synchronization, using Caulobacter as a specific example. Our model takes into account the initial population asynchrony and, similar to the yeast cell cycle phase probability density model presented in Orlando et al. ^[12], captures the phase variability resulting from asymmetric cell division and differences in cell cycle times. An appealing aspect of our model is its simplicity; a knowledge of three easily-measured parameters—namely the mean SW-to-ST transition phase (or equivalent), division time COV, and SW/ST cell total volume fraction (or equivalent)—and the initial synchronization state (i.e., the cell-type distribution at the outset of an experiment) are all are that is required to describe the time-dependent cell-type distribution.

The aforementioned parameters and initial synchronization state are specific to a given model system and experimental condition. For a synchronized population of Caulobacter under normal growth conditions, we use a mean SW-to-ST transition phase of ∼0.25, division time COV of 0.13, cell volume partitioned 40% SW to 60% ST, and a simulated initial cell cycle phase distribution that accurately models the real synchronization process. But Caulobacter is not the only synchronizable model system to which our cell-type distribution model can be applied. Indeed, synchronizable model systems are found across the tree of life, including E. coli^[42], S. cerevisiae^[43], and mammalian cells ^[44]. A 1957 review by Campbell describes synchronization methods for 11 microbial species ^[45]. For the symmetrically dividing E. coli, the equivalent of the SW-to-ST transition phase would be set to zero, and the two daughter cells would (on average) have the same volume. In the case of S. cerevisiae, the SW-to-ST transition phase equivalent is equal to the average fraction of the cell cycle that the budded daughter cell remains in the early G1 stage ^[46], with the average size of the budded cell being smaller than that of its mother ^[47]. The division time COVs for a number of commonly studied systems have already been published (a compilation of these values can be found in ^[1]). Initial cell distributions for many of these organisms have to be determined.

We note that we have assumed a perfect Caulobacter synchrony, i.e., exactly 100% of the cells at the beginning of the experiment are SW cells. In real cell synchrony experiments, SW fractions are close but not necessarily equal to 100% (see, e.g., ^[22]). However, minor differences in the purity of a synchronized population are not expected to significantly alter our results. That our cell-phase distribution model is consistent with experimental observations of the time-dependent state of a Caulobacter population (Fig. 4) supports this assumption.

Along with characterization of cell distribution, there has been considerable interest in recent years in extracting “single-cell”-like information from population data using deconvolution-type algorithms ^[13]–^[15],^[48],^[49]. Although all algorithms of this kind are somewhat limited in the level of detail they can provide about biological systems—at best, only synchronous average information, and not the full stochastic variability between cells at identical phases, can be determined—they have been highly effective at uncovering features not visible in the population data. The model-based deconvolution method presented here is an extension to these previous methods and a powerful tool for the analysis of biological data, requiring no more information than the parameters described previously, and is applicable to any time-series data set for which the state of the synchrony is known or can be predicted. In particular, our method can be applied to time-series gene expression data to identify additional cell cycle-regulated genes not previously discovered and to complete meta-analyses across multiple platforms (i.e. competitive hybridization oligo arrays or non-competitive hybridization arrays such as Affymetrix). Although the differences in the data obtained from different platforms may require modifications to the kernel function, the method itself is independent of the experimental and biological details; indeed, the method supports arbitrary kernel functions.

Even with a detailed and accurate kernel and an accepted deconvolution-type algorithm, the precise shape of a deconvolved function is in general highly sensitive to the value of the regularization parameter ( in this work; see Eq. (12)). To objectively address this problem, we employ a cross-validation routine that provides a sensible and well-established criterion for determining the appropriate amount of regularization. Our use of cross-validation in deconvolution of time-series gene expression data thus represents an improvement over methods that use arbitrary regularization based only on visual inspection of the estimated profiles.

By construction, the model-based deconvolution method presented in this paper mitigates the effects of synchronization loss in expression experiments. However, as with all time series experiments, the estimates remain dependent on the sample rate of the data. If the sample rate is insufficiently high to capture salient gene activity, important events in the expression profile may be missed. In principal, lower sampling rates may be accommodated by increasing the number of assumptions made about the expression profile to be estimated. In this paper, smoothness (Eq. (12)), positivity (Eq. (13)), and continuity (Eq. (14)) were all used to decrease the effective degrees of freedom and supply a maximal, yet realistic, amount of a priori information. The cubic splines support a broad class of potential expression functions, however more restrictive models could be used to supply stronger assumptions and support lower sampling rates—at the cost of potentially being overly restrictive and not capturing the true gene expression profile. See, e.g., ^[50] for further consideration of sample rates in temporal data.

The synchronous average expression profiles extracted using our generalized deconvolution algorithm are, with the effects of population asynchrony removed, a much-improved reflection of biological reality. We demonstrated this with Caulobacter, calculating deconvolved expression profiles for 10 genes previously found to be cell cycle-regulated and essential for cell viability or polar cell development (Fig. 5). As mentioned in Results, the deconvolved expession profiles generally have their peaks shifted to later times relative to the population data. This is to be expected, since even a perfectly-synchronized population at the outset of an experiment contains both young SW cells ( ) and old SW cells (and all cells in between). Many of the genes analyzed here also show a narrowing of their expression peak(s) following deconvolution, although this is not universally true. The expression profile of divJ, for example, is shifted to later times but not otherwise fundamentally changed; the peak, located just after the SW-to-ST transition in the deconvolved profile, is as broad as in the population measurement. Thus, expression peak narrowing is not an artifact of the deconvolution method, but rather a property of an individual gene's expression profile. Here we highlight some of our Caulobacter-specific results that also demonstrate the power of combining an organism-specific kernel with a generalized deconvolution routine:

ctrA. As the master regulator of the Caulobacter cell cycle ^[51], ctrA is arguably the most well-characterized of Caulobacter genes. It has been shown that ctrA expression is controlled by two promoters (P1 and P2) that are differentially-regulated by phosphorylated CtrA (CtrA∼P): the weaker P1 is negatively-controlled by CtrA∼P and the stronger P2 is positively-controlled (Fig. 7A). The P1 promoter is activated in the early ST cell, immediately following replication of the chromosomal ctrA locus. Activation of the weak P1 promoter leads to an increase in the CtrA∼P concentration, which then activates the stronger P2 and represses P1 ^[52]. The differential regulation can be seen in Fig. 7B, left panel (data reproduced from Reisenauer and Shapiro ^[53]). Although these details are not visible in the population-level microarray data, they are revealed in the deconvolved expression profile (Fig. 7B, middle and right panels). For example, in the deconvolved profile, ctrA expression remains flat until DNA replication is initiated at the SW-to-ST transition. Perhaps most interestingly, the initial expression ‘shoulder’ is consistent with expression from P1, and the main peak beginning around the phase of cell compartmentalization (transition from EPDLPD), is consistent with expression from P2. The shape of the deconvolved ctrA profile is thus validated by our previous knowledge of the mechanism of ctrA regulation.

ftsZ. The tubulin homolog FtsZ is essential for bacterial cell division. It has been shown that transcription of ftsZ is repressed in SW cells and activated only when the DNA replication begins ^[20]. However, this regulation is not clear from the microarray data alone. Specifically, the raw microarray data show no delay in ftsZ transcription from the time the experiment begins (Fig. 8, left panel). In contrast, the deconvolved expression profile reveals the delay in transcription initiation until the beginning of the ST stage (Fig. 8, right panel), consistent with our understanding of ftsZ regulation.

divK and ccrM. DivK is an essential single-domain response regulator that is transcriptionally-activated by CtrA∼P and plays a role in the cell cycle-regulated proteolysis of CtrA ^[54]. The essential ccrM DNA methyltransferase gene ^[55] has an expression profile similar to that of divK. In both cases, deconvolution reveals that expression begins in the EPD cell, and that the change from zero to maximal expression happens over a much shorter time (i.e., the response is more switch-like) than is evident from the population data.

cckA. One of the more interesting results is the predicted transcription profile of cckA, which encodes an essential histidine kinase responsible for CtrA phosphorylation ^[56]. The population-level microarray measurements show a single expression peak approximately half-way through the cell cycle, while the deconvolved profile shows two peaks: one beginning at the SW-to-ST transition and another peaking in the EPD cell. Although this result has not been previously reported, it does suggest the interesting possibility that cckA is under the control of additional and unknown layers of transcriptional regulation during the cell cycle.

These deconvolution results appear to be relatively insensitive to changes in model parameters. Of the parameters used in the cell cycle phase distribution model, the mean SW-to-ST transition phase is the one that is known with the least certainty. However, we found that precise knowledge of the mean transition phase under a given condition is not absolutely necessary for extraction of average single-cell data with our deconvolution algorithm. Even a substantial change in the assumed SW-to-ST transition phase had only a small effect on the deconvolved profiles. With respect to the single-cell volume model employed in the deconvolution algorithm, even the extreme and false assumption of fixed cell volume had an insignificant effect on the shape of the deconvolved expression profile.

One Caulobacter-specific result that merits further discussion is the SW-to-ST transition phase. Although accepted as around 0.25, or even up to 0.33, for standard growth in a rolling tube or shaken flask ^[17],^[22],^[39],^[40], it can change under other conditions. We present data showing that the mean transition phase is reduced to 0.15 in a microfluidic environment in which the cells are rapidly growing. We recognize that a possible explanation for this low value may be that the timing of the SW-to-ST transition in our microfluidic growth experiments is skewed by a division control system in which ST cells that have just transitioned from the SW stage divide on a different time scale than ST cells that follow from cell division. However, we are not aware of any data that would suggest that this is the case. Indeed, the morphology of ST cells after the transition from SW cells appears to be the same as the morphology of ST cells after division, and a single mean SW-to-ST transition phase in our model is consistent with experimental observations (Fig. 4). Furthermore, given that a population of Caulobacter cells starved for carbon or nitrogen tend to arrest during the SW phase ^[57],^[58], it is likely that the SW-to-ST transition phase can both increase and decrease, and be well above 0.33 under less-favorable environmental conditions. That the timing of this cell cycle ‘checkpoint’ may vary with growth conditions is a fascinating result that deserves more detailed study.

To our knowledge, our deconvolution method is the first to specifically deal with the unique analytical challenges posed by dimorphic organisms. Although this method can be applied to any time-series measurement made on a cellular population, we have demonstrated its utility with an analysis of cell-cycle regulated gene expression in Caulobacter. Certainly, directly measuring the concentration of individual transcripts in real time in single cells remains the gold standard in quantifying the gene expression behavior of single cells; the insights provided by such real-time, single-cell studies of mRNA have been profound ^[59]–^[62]. Still, despite recent progress and a number of successes, the real-time measurement of mRNA in single cells remains a challenging problem. Our method allows for the simple analysis of mRNA concentrations measured with common laboratory tools and advances the performance of population-level methods closer to that of single-cell studies. Thus, combining high-throughput experimental expression data with novel computational algorithms can provide new and exciting insights into the function of cellular systems.

The authors have declared that no competing interests exist.

DS-G is supported by the MBI, which receives major funding from the National Science Foundation Division of Mathematical Sciences and is supported by The Ohio State University. The Mathematical Biosciences Institute adheres to the AA/EOE guidelines. SC is supported by a Beckman Young Investigator Award from the Arnold and Mabel Beckman Foundation and by a grant from the Mallinckrodt Foundation. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The authors would like to thank Malgorzata Rowicka-Kudlicka, Andrzej Kudlicki, and Patrick McGrath for helpful discussions, and Alison Hottes for valuable comments on the manuscript.