Correlations of single-cell division times with and without periodic forcing

Periodic forcing of nonlinear oscillators leads to a large number of dynamic behaviors. The coupling of the cell-cycle to the circadian clock provides a biological realization of such forcing. Using high throughput single-cell microscopy, we have studied the correlations between cell cycle duration in discrete lineages of several different organisms including those with known coupling to a circadian clock and those without known coupling to a circadian clock. Correlations between cell cycles duration in discrete lineages observed in the organisms with a circadian clock cannot be explained by a simple statistical model but are consistent with predictions of a biologically plausible two dimensional nonlinear map. Surprisingly, the nonlinear map is equivalent to a classic nonlinear map called the fattened Arnold map. The model predicts that circadian coupling may increase cell to cell variability in a clonal population of cells. In agreement with this prediction, deletion of the circadian clock reduces variability. Our results show that simple correlations can identify systems under periodic forcing and that studies of nonlinear coupling of biological oscillators provide insight into basic cellular processes of growth.

Periodic forcing of nonlinear oscillators leads to a large number of dynamic behaviors. The coupling of the cell-cycle to the circadian clock provides a biological realization of such forcing. Using highthroughput single-cell microscopy, we have studied the correlations between cell cycle duration in discrete lineages of several different organisms including those with known coupling to a circadian clock and those without known coupling to a circadian clock. Correlations between cell cycles duration in discrete lineages observed in the organisms with a circadian clock cannot be explained by a simple statistical model but are consistent with predictions of a biologically plausible two dimensional nonlinear map. Surprisingly, the nonlinear map is equivalent to a classic nonlinear map called the fattened Arnold map. The model predicts that circadian coupling may increase cell to cell variability in a clonal population of cells. In agreement with this prediction, deletion of the circadian clock reduces variability. Our results show that simple correlations can identify systems under periodic forcing and that studies of nonlinear coupling of biological oscillators provide insight into basic cellular processes of growth.

I. INTRODUCTION
The process of cell division has fascinated scientists since the invention of the microscope. During the pro- * nathalie.balaban@mail.huji.ac.il cess of cell division in organisms that divide symmetrically, a cell generates two almost identical copies of itself. Many mechanisms act in concert to enhance the fidelity of replication and division, including proofreading, DNA repair enzymes, an elaborate partitioning apparatus and feedbacks. The statistics of the cell division process in single cells has been proposed to provide an unbiased way to uncover the type of feedback that control the process of cell division ( [1,2] and therein) and have motivated researchers to gather as much data as possible. Powell, one of the pioneers of single-cells measurements, was described as sitting in a 37 • C room for many hours watching bacteria divide and recording manually cell-division events [3]. Recent technological advances in microscopy and microfluidics [4,5] have boosted our ability to gather information over tens of thousands of cells and opened the door to quantitative analyses of the process of celldivision on lineages [6,7].
Cell division is a discrete process ( Fig. 1 and Movie 1) which occurs at each generation. The cellular components inherited from the previous division govern the state of the cell at birth. Therefore, it is appealing to describe this process with discrete maps. For maps, consecutive timing of key events such as the duration of sleep [8], phase of cardiac firing in a stimulation cycle [9], or cell cycle duration [10], depend in an iterative manner on the timing of previous events. In most studies involving maps, one assumes that the map is identical under subsequent iterations. However, assuming cell division can be described by a map, the map itself would be duplicated subjected to noise [10,11], generating not only one time-series but additional branches at each generation. The information contained in the lineages is much larger than on a single branch [12,13]. In the case of non linear maps, rich dynamics can emerge from the iterative process, but these dynamics have rarely been analyzed on lineages [11]. If the replication process is perfect, in principle no new information is provided by the branches. However, noise is always present and introduces new information on how the process propagates along the branches. Thus, the study of cell lineages introduces a new class of problems involving dynamics of nonlinear maps that are themselves duplicated with noise in subsequent generations.
Our goal is to formalize the expected statistics of inter and intra-generation correlations for discrete time-series and compare our results to experimental data measured on lineages of single-cells (Fig. 1). For this purpose, we focus on the measurement of the duration of the cellcycle, which is the time between two consecutive cell divisions, and its inheritance along lineages which generates the time-series T n and its branches. We expand on our previous work analyzing these correlations in mouse lymphocytes [10] and compare the data in different organisms, with a without putative periodic forcing. We analyze the correlations between cells in the lineages using a discrete map representing the periodic forcing of the cell-cycle by an external non-linear oscillator [10].

I.1. Periodic forcing
One of the classic problems in mathematics involves the effects of period forcing of a nonlinear oscillator [14]. The nonlinear oscillator can be represented by a differential equation containing a stable limit cycle. The periodic forcing is typically either a continuous periodic input or a pulsatile stimulus. In either case, one expects the appearance of certain universal features that can often arise independent of the detailed equations or the nature of the stimuli. If the intrinsic period of the driven oscillator is sufficiently close to the period of the forcing oscillator, then there will typically be entrainment or 1:1 phase locking where the two oscillators are synchronized. As one changes the relative frequencies of the oscillators and the strength of influence of the periodic forcing on the driven oscillator, then many different behaviors can arise. One such behavior is n : m phase locking in which there is a stable periodic rhythm in which there are n cycles of the forcing oscillator to m cycles of the driven oscillator. Two different types of aperiodic rhythms are possible. These can be distinguished by considering what happens to the trajectories starting from two different nearby phases of the driven oscillator as time proceeds. If the distance between the two trajectories grows with time, then the dynamics are chaotic and if the distance stays approximately the same, then the dynamics are quasiperiodic [15]. Moreover, many features of the locking zones have a universal structure. For low forcing amplitudes, there are typically zones of stable phase locking, called Arnold tongues, where the ratio m/n monotonically increases as the ratio of forcing oscillator period to the intrinsic period of the forced oscillator increases. These basic insights emerge from research stretching back to Poincaré with major insights from Arnold [16], Smale [17] and many others. Study of periodically forced nonlinear oscillators are not only of interest in themselves, but they can also help understand dynamics in a wide range of physical [14,18,19] and biological systems [9,[20][21][22].
One such problem involves the effects of an external rhythm on the cell cycle duration, T n , [10,22,23]. This external rhythm can be externally imposed, leading to new insight on the cell-cycle itself [24], or originate from a clock encoded within the cell, such as the circadian clock [25]. The circadian clock is an internal cellular oscillator which has an approximate period of 24 hours and which can influence cellular processes depending on the time of the day [25]. For example, in some organisms that have a circadian clock, there are time intervals in the day during which the progression through the cell cycle proceeds slowly. This phenomenon is called gating [26,27]. It is widely believed that in the presence of a gate, cell cycle states synchronize to the circadian signal [22]. Since biological systems typically have considerable variability, even in situations in which there is believed to be synchronization, there can be considerable variation in the timing of the events, and quantitative analysis and modeling is needed to test hypotheses.
Two types of theoretical models have been proposed to interpret the experimental results on entrainment of oscillations: nonlinear differential equations and nonlinear maps. Nonlinear differential equations are often developed specifically for particular systems with parameters generally determined by optimizing fits to data [22,28]. In situations in which the differential equations are based on realistic models, as for example in models of the interactions of the circadian clock and the cell division cycle, it is possible to determine some of the parameters based on different sets of experiments than those used to model entrainment data [23]. Maps constitute an alternative type of model.
In what follows we discuss the synchronization of the cell cycle to the circadian rhythm in the context of nonlinear dynamics models of periodic forcing. In order to simplify previous approaches and reduce the number of parameters, we apply a stochastic nonlinear map to study cell cycle time correlations in lineages from several different species. The analysis demonstrates the importance of deterministic non-linear factors in controlling the cell cycle and also suggests new directions for theoretical analyses. In Section II, we present the mathematical model of periodic forcing. Section III gives the experimental results, mathematical analysis of a null model and the fitting of the mathematical model of periodic forcing to the data. Section IV is the discussion of the results, and in Section V we present the experimental methods.

II. MATHEMATICAL MODEL
In recent work, some of us proposed a biologically plausible model for the interaction of the circadian rhythm with the cell-cycle, called the "kicked cell-cycle" [10].
This model is based on the assumption that the cycle duration of a daughter cell depends linearly on the cycle duration of the mother cell, as well as on the circadian time of the cell division. The basic idea is that the cell-cycle duration, T n , of a cell in generation n can be influenced both by the cell-cycle duration of its mother T n−1 and the phase at its birth of a forcing oscillator such as the circadian clock. Given the birth time of a cell in generation n as t n , the phase in the forcing rhythm is t n /T osc , where T osc is the period of the forcing oscillator (approx. 24 h for the circadian clock). The model for analysis is: where t n represents the birth time of a cell in generation n, T n represents the cell-cycle duration of a cell in generation n, τ 0 is the intrinsic cell-cycle period in the absence of circadian forcing, k is a parameter that controls the magnitude of the coupling between the cell-cycle and the circadian oscillation, α ∈ [−1, 1] is a parameter that allows a tuning for the influence of the mother cellcycle duration on the current cell-cycle duration and ξ ± n is white noise with ξ ± n = 0 and ξ ± n ξ ± m = ξ 2 δ n,m δ +,− , where · denotes average over realizations. The superscripts (±) designate the two sisters cells. For the situation in which ξ = 0, both sisters have the same cell-cycle duration. In this case, the model has a rich history. For α = 0 it is equivalent to the Arnold circle map [16], for α = 1 it is the standard map [29], and for general values of α it is the fattened Arnold circle map [30] also known as the dissipative kicked rotor [31,32]. Depending on the parameters, the dynamics of the model include fixed points (for example around τ 0 = T osc ), regions of periodic, quasi-periodic or chaotic behavior.
Despite the large amount of theoretical work on discrete maps that generate time-series, very little is known about the statistics of the time-series data that emerge on lineages created by the self-replicating process in the presence of noise. An early exception was a study of cellcycle variability by Grasman [11] that focused on modeling the correlation between the cell cycles of mother and daughter cells using the logistic map [33]. More recently, our group has proposed the kicked cell cycle map [10] to model lineage data. In the current paper we expand on our previous work by analyzing cycle duration correlations in several organisms and determining the values of the parameters in the kicked cell-cycle model giving a good fit to data in each organism.

III.1. Experimental correlations in lineages of cells
We measured or analyzed the cell-cycle duration in thousands of single cells and in various organisms: E.coli, Corynebacteria, Cyanobacteria, EMT6 human cells and L1210 mouse lymphocytes (Table I). While the first two organisms are not known to be controlled by an external oscillator, the circadian clock coupling to the Cyanobacteria cell-cycle has been extensively studied and a similar coupling has more recently been suggested to be active in mammalian cells [25]. Corynebacteria were chosen because of their division mode which occurs by snapping (see Movie 1) thus reducing the experimental noise in the determination of the cell division event.
Mother cells divide into two daughter cells, called sisters. The sisters again divide into four cells. Two daughter cells from different sisters are called cousins (Fig. 1). As shown in our previous work on the analysis of the L1210 data [10], meaningful information can be gained from measuring the correlations between different cells within the same lineage. In particular, we measured the Spearman correlations in cell-cycle duration between sister cells ρ s−s , between mother and daughter cells ρ m−d , and between cousin cells ρ c−c , (See Table I for measured correlations). The averages in the correlation coefficients are over different lineages. We found significant correlations between sister cells in all data sets in agreement with earlier results [34]. In most data sets, the ρ m−d was low or non-significant, suggesting that the memory of the cell-cycle duration was lost within one cell-cycle. The coefficient of variation of the cell-cycle duration varied from 10% to 40% in the various data sets, similarly to typical values in the literature. This variation is larger, for example, than the one expected from the measured noise between sisters, ξ, according to a simple model of inheritance, the bifurcation auto-regression (BAR) model [35]. Another departure from this model is the observation that cousins often had correlated cell-cycle duration which could not be attributable to microenviromental conditions [10]. Intuitively, the correlation between cousin cells despite the absence of correlation between mother and daughters is surprising. In order to formulate this intuition more rigorously and for the general case, we consider a process that proceeds in a tree-like fashion, as the division process does (Fig. 1A).

III.2. Expected correlation of cell-cycle duration in a lineage
Similarly to the assumptions of the BAR model, we assume that the fate of a daughter cell is determined by inherited factors from its mother cell and that there are no external influences. Likewise, we assume that the fate of two sister cells are directly related due to the resemblance in their composition at birth. Thus the correlations between cells are determined by the lineage relations.
In order to neutralize the effect of the fate of cell Z in Fig.1A, on the correlation of cells X and Y , we use the partial correlation function. The partial correlation ρ XY,Z between random variables X and Y, removing the effect of the variable Z, is: We now return to the labels of Fig. 1A and recall that we expect zero correlation between the fates of X and Y when the effect of Z is removed. For Spearman correlations, this assumption is valid if the correlations reflect monotonous relations between cells. Thus, ρ XY,Z = 0, and Eq. 3 implies: Likewise, we expect zero correlation between K and Y when the effect of X is removed, so ρ KY,X = 0. Eq. 5 is obtained directly from Eq. 4 by changing labels (X → K,Z → X), and by substitution of Eq. 4: Replacing the labels KY,KX,XZ and YZ by the cells relations (as in Fig. 1A) leads to: This derivation is valid for any cell fate measured on a lineage. In particular, Eq. 6 is the expected relation between the cell-cycle correlations ρ c−c , ρ s−s and ρ m−d (using Spearman coefficients), under the assumption of at most a monotonic dependence of a cell-cycle on its mother cell-cycle and on its sister cell-cycle (see also [35]). We note that ρ c−c is expected to be always smaller than ρ m−d . For the data sets of Corynebacteria and of E.coli, Eq. 6 holds quite well. However, for the Cyanobacteria, lymphocytes and EMT6 data sets we observe a large deviation from this expected behavior ( Fig.  2 and Table I), as quantified by the parameter ∆: Moreover, in several data sets, contrary to expectations we observe empirically an even stronger departure from Eq. 6, that we termed the cousin-mother inequality:

III.3. Fitting of the mathematical model to the data
The kicked cell cycle model offers a potential explanation for these results. As the intrinsic cell cycle period, τ 0 , the forcing amplitude k, and the mother-daughter coupling α vary, different dynamical behaviors and bifurcations occur. In Fig. 3A, the complex landscape of periodicity regions coming from the kicked cell-cycle model is shown for α = −0.5. Despite the classical nature of this problem, there is still comparatively little known about the bifurcations when α = 0 and α = 1. We  observed a strong similarity to the locking zones in the sine circle map [20,21,38], and further study is needed to investigate the differences for various α . From our context, we are most interested in the cousin-mother inequality. In Fig. 3B-C we plot the simulated values of the cousin-mother inequality for noise level ξ = 0.01 (Fig.  3B) and ξ = 0.1 (Fig. 3C). In a large part of the plot, the cousin-mother inequality holds, showing that the model accounts for the departure of the experimental data from Eq. 6, because of the non-monotonous sinusoidal term in the kicked cell-cycle model. This plot demonstrates that the Arnold tongue structure plays a strong role in defining the value of the correlations, even in the presence of a substantial noise level.
In order to further test the validity of the kicked-cell cycle model in describing the experimental data sets with high ∆ (L1210, EMT6 and Cyanobacteria), we extracted for each of these data sets the 6 measurements listed in Table I and fitted them to Eqs. 1-2 (see section V for more details). The noise level of each data set was directly evaluated from the variation between sister pairs (see ξ in Table I, and section V for the derivation). The fit resulted in values for τ 0 , k and α (Table II). The values of τ 0 are close to the measured mean values, but not identical, as expected from the model, where τ 0 has a strong effect on the mean value of the cell cycle and may depend on growth conditions. One interesting feature that comes out of the fitting procedure, is that the cell-cycle distributions over the population are close to the measured distributions. One striking example is the fitting to the Cyanobacteria data (Cyanobacteria 1) that displays a bimodal distribution, and this bimodality is apparent also in the simulated distribution without having been fed directly into the fitting procedure (Fig. 4A).

III.4. Correlations of cell-cycle duration in a mutant deleted for the clock genes
The model suggests that the coupling of the cell-cycle to the circadian clock can, depending on parameters values, either synchronizes the cells and reduce cell-to-cell variability [39], or increase the variability of the cellcycle duration in the population. This suggests, that driving the cell-cycle by the circadian clock outside the fixed points regions should lead to enhanced variability in cell cycle duration, and therefore reducing the coupling to the circadian clock should reduce this variability. In order to test this prediction, we compared cyanobacteria cell-cycle variability for the wild-type (WT) strain and for a mutant deleted for the circadian clock. Here we present new data sets of WT strain of Cyanobacteria and its derived clock mutant (Cyanobacteria 2 and Mutant in Tables I-II), as these two data sets were measured in the same conditions and have similar mean cell-cycle (Table I). Interestingly, we observe that the variability in the cell cycle is significantly reduced (F-test, p-value < 0.05) in the mutant strain (CV = 0.12) compared to the WT strain (CV = 0.22) (Fig. 4B-C). Accordingly, whereas the cousin-mother inequality was fulfilled in the WT strain, it was not significant in the mutant strain ( Table I). The parameters of growth of WT Cyanobacteria 2 are different from WT Cyanobacteria 1. Accord-  ingly, the distribution of cell cycle duration is different and predicted by the model to be unimodal (Fig. 4B).

IV. DISCUSSION
One main feature of the dynamics of the kicked cellcycle model is the observation that the Cousin-mother inequality is obtained for quite a wide range of parameters, even when the noise level is relatively high (Fig.  3C). However, there are still regions of parameters where the Cousin-mother inequality is not fulfilled, for example at fixed points. Also, higher noise level may eventually mask the deterministic periodic forcing as well as the inequality. Therefore, whereas the cousin-mother inequality strongly suggests the existence of a non-linear driving mechanism, its cannot disprove its existence. Although the absence of the inequality in E.coli and Corynebacterium data is consistent with the absence of a known internal clock, a putative clock may have been masked by the factors mentioned above. Data of E.coli grown on poorer medium [3,40] shows the Cousin-mother inequality and future work is needed to determine the underlying biological mechanism.
Comparison between the periodicity plot and the Cousin-mother inequality plot reveals that both display  Tables I-II), grown in different conditions leading to faster growth, with unimodal distribution. (C) Mutant of the strain shown in (B) deleted for the clock gene (red), displays a significantly smaller coefficient of variation (CV) compared to the WT (blue).
the "Arnold tongues" features. The periodicity plot presented in Fig. 3A shows the geometry of the locking zones as the period τ 0 and the amplitude of the periodic forcing k are varied. One prediction is that changing parameters to move from a period 1 (fixed point) to a period 2 zone, frequency doubling of the cell-cycle may be observed, which may be related to recent observations [41]. Comparing the periodicity plot to the plot of the Cousin-mother inequality (Fig. 3C) reveals that, even in the presence of high noise level (ξ = 0.1) the period 2 and fixed points features are retained. At lower noise levels, it is interesting to follow which additional features are robust to noise. For example, the regions of intersections are very prominent on Fig. 3B. Therefore, the strength of the Cousin-mother inequality may serve as a straightforward experimental measurement that can reveal non-trivial features characteristic of discrete maps in empirical data.
In contrast with the periodicity, Lyapunov exponent, Correlation Dimension or other indicators used to analyze chaotic behavior that require long data sets, the Cousin-mother inequality does not require extremely large quantities of data to detect deterministic contributions of the non-linear process to the dynamics of the system. Therefore, the Cousin-mother inequality could be useful as a mean to detect non-linear coupling on lineages of cells. However, the analysis should be done only after experimental artifacts that may lead to similar observations are ruled out. In particular, spurious spatial correlation between cells due to microenvironment should be ruled out by analysis of spatial dependencies. Furthermore, departure from Eq. 6 may be due to noise and insignificant. Therefore, analysis of the confidence interval for the departure from Eq. 6 should be done (see section V), as represented by the grey areas in Fig. 2. Once these external influences are ruled out, the Cousinmother inequality can be a powerful tool for revealing the effect of non-linear coupling on the cell cycle variability. However, it should be noted that for determining a significant inequality (i.e. outside the grey area in Fig. 2), a large enough number of independent lineages should be analyzed. For example, another data set of Cyanobacteria was analyzed [42], showing features consistent with vicinity to the fixed point in the model (mean cell-cycle close to 24h), however the small number of lineages (< 20) resulted in non-significant correlations and the fitting procedure could not be performed. It should be noted that our analysis is not restricted to measurements of the cell-cycle duration. In effect, any observable that is measured on a self-replicating entity may display the Cousin-mother inequality, provided that its variability is governed by a non-linear process. Therefore, we expect that the Cousin-mother inequality could be used as a general indicator of deterministic non-linear processes along lineages.
In this work, the cousin-mother inequality has been used to unveil a strong deterministic component in the variability of cell-cycle duration of Cyanobacteria. In agreement with our understanding, the deletion of the clock genes results in a significantly reduced cell cycle variability in the population. In most theoretical analyses, such variability is understood as the inevitable consequence of cellular noise [43,44]. Here we show that deterministic factors, such as periodic forcing, can increase the variability of the cell cycle in a way that can be controlled. In bacteria, the cell cycle variability has been shown to have important consequences for the survival of populations under stress [4], as well as for their ability to evolve resistance factors [45]. Therefore, understanding the source of this variability is important for predicting the behavior of bacteria under antibiotic treatment, as well as their ability to evolve. The cells were imaged in a Leica automated fluorescence microscope system, as previously described [10]. Briefly, a polydimethylsiloxane (PDMS) square mold was filled with medium (L-15) and sealed with a coverslip. Illumination was kept low enough to show no influence on total cell-cycle duration.
Analysis of division times was done using phasecontrast images and custom image analysis software (1). We used an automatic cell-tracking platform written in MATLAB (MathWorks). Cell cycle duration was determined from phase-contrast images acquired at 5 min intervals. Together with the sharp division process of L1210 cells, this resulted in less than 1 − 2 % experimental noise in T n .

V.1.2. E.coli
In the E.coli experiments we used the same system described for the lymphocytes. Here the PDMS square mold was filled with melted LBagarose, that was prepared from LB Broth, Lennox (Difco) (LBL). Images were acquired using a 100X 1.3NA oil objective and a cooled CCD camera (Orca; Hamamatsu). Microscopy was carried out at 37 • c. The phase-contrast images were acquired at 1 min intervals. The images were tracked semi-automatically for supervised analysis using a plugin developed for ImageJ [10].
Prior to microscopy single colonies were diluted into fresh LBL. Cells were grown overnight at 37 • c, with shaking. Cultures were supplied with 15% glycerol and stored at -80 • c. The frozen cultures were diluted into fresh LBL and grown at 37 • c, with shaking, for 3 hours.

V.1.3. Cyanobacteria 2
Parameters used in this study were extracted from time-lapse movies of Synechococcus elongatus wild type cells (ATCC strain 33912 T M ), and clock deletion mutant (∆kaiBC ) cells grown under constant light conditions [37]. Clock deletion cells carry an insertion of a gentamycin resistance cassette into the ORF of the kaiBC operon.
For time-lapse microscopy movies of Cyanobacteria, cells were first grown in liquid BG-11 media at 30 • c with constant rotation. The ∆kaiBC strain was supplemented with gentamycin at 2 µgml −1 . Constant light levels were maintained at approximately 20 − 25µEm −2 s −1 by cool fluorescent light sources. Cell cultures were kept at exponential phase and entrained by subjecting cells to a 12 h light : 12 h dark cycle. Nikon Ti-E inverted microscopes, equipped with the Nikon Perfect Focus System module, were then used to acquire time-lapse movies of growing micro-colonies over several days under constant light at 15µEm −2 s −1 , using a protocol adapted from [41]. Illumination for photoautotrophic growth was provided by a circular cool white light LED array (Cairn Research, UK), attached to the condenser lens. Images were acquired every 45 minutes using a CoolSNAP HQ2 camera (Photometrics, Arizona), and a 100X objective.
Movies were segmented and tracked using a modified version of Schnitzcells [46]. Finally, cell lineages were reconstructed by tracking individual cells across frames, and identifying mother-daughter relationships in division events. For full description of methods and data see [37].

V.1.4. Corynebacteria
Experiments were carried out (as in [47]) using an inverted time-lapse live cell microscope (Nikon TI-Eclipse, Nikon Instruments, Germany) equipped with a 100x oil immersion objective (CFI Plan Apochromat Lambda DM 100X, NA 1.45; Nikon Instruments, Germany) and a temperature incubator (PeCon GmbH, Germany). Cultivations were performed at 30 • C. Phase contrast images of growing microcolonies were captured every 2 minutes using an ANDOR Clara-E CCD camera (Andor Technnology, UK).
Time-lapse movies were analyzed using a custom, specialized workflow implemented as an ImageJ/Fiji plugin. Cell identification was performed using a segmentation procedure tailored to detect individual rod-shaped cells in crowded populations. Detected cells were subsequently tracked throughout all image sequences using an adapted single particle tracking approach as implemented in TrackMate. Derived quantities, i.e., growth rates, were computed using the Vizardous framework [48].
Cells were pre-cultured as 20mL cultures in 100mL baffled Erlenmeyer flasks on a rotary shaker at 120rpm orbital shaking at 30 • C. A first pre-culture in BHI (Brainheart infusion, Becton Dickinson, USA) complex medium was inoculated into a second pre-culture in CGXII mineral medium, which was finally inoculated at OD600 = 0.05 into CGXII mineral medium, the main culture.

V.2. Microfluidic devices-Corynebacteria
Polydimethylsiloxane (PDMS) (Dow Corning; Farnell GmbH, Germany) microfabrication was carried out [47] to manufacture single-use microfluidic devices with integrated 10µm high supply channels and cultivation chambers with a height of 1µm. A 100mm silicon wafer (Si-Mat, Silicon Materials, Germany) was spin-coated separately with two layers of SU-8 photoresist (Micro Resist Technology GmbH, Germany), processed by photolithography. This silicon wafer served as a reusable mold during subsequent PDMS casting. Thermally cured and separated PDMS chips were treated with oxygen plasma and permanently bonded to 170µm thick glass slides (Schott, Malaysia) just before the experiments. Manually punched inlets and outlets were connected with tubing (Tygon S-54-HL, ID=0.25mm, OD=0.76mm; VWR International) via dispensing needles (dispensing tips, ID =0.2mm, OD=0.42mm; Nordson EFD Germany).
Fluid flow into the microfluidic chip was controlled with a 4-fold NeMESYS syringe pump (Cetoni GmbH, Germany). A detailed setup protocol can be found in [50]. Prior to microfluidic cultivation, the microfluidic chip was purged with fresh and sterile-filtered CGXII medium at 200nl/min for 10min. Afterwards, the chip was infused with bacterial suspension for single-cell inoculation as described in full detail recently. Bacterial suspensions were withdrawn from the main culture at the exponential growth phase (OD600 between 0.5 and 1). As soon as sufficient single cells were inoculated into the microfluidic cultivation chambers, solely CGXII medium was infused at 200nl/min throughout the entire cultivation.  [47] Corynebacterium glutamicum ATCC 13032 V.4. Computation of correlation coefficients All correlation coefficients are Spearman coefficients (similar results were obtained using Pearson coefficients). ρ s−s and ρ c−c were calculated on the 3rd generation, while ρ m−d was calculated on the 2nd and 3rd generations (the mother from generation 2 and the daughter from generation 3). To avoid spurious dependencies, we were careful to include only one pair of cells chosen randomly from each cell lineage. The correlation coefficients were computed on the chosen cell cycles. This computation was repeated 1000 times, and the final correlation coefficients are the averages of those repetitions, while the errors of the correlations were taken to be the standard deviation, σ.

V.5. Significance computation
We simulated cell cycles of random data. The simulated cells were taken from a normal distribution with experimentally measured CV (coefficient of variation). Each simulation contained 60 lineages (typical number of lineages we track in each experimental set), and 7 cells per lineage (three generations: mother, 2 sisters and 4 cousins). By determining the covariance matrix of mothers and their two daughters, we matched ρ s−s and ρ m−d to the experimental results (separate simulation to each experiment). We computed all correlations from the simulated data, and made 100 simulations for each experiment. We determine departure from Eq. 6 as significant (with p-value< 0.05) if the experimental ∆ (or the experimental value of ρ c−c − |ρ m−d |) is larger than 2 standard deviations of the ∆ (or ρ c−c − |ρ m−d |) that is obtained from the random simulations described above. The grey area in Fig. 2 was computed using 500 simulations, similar to the described above. Here we used constant typical values for CV and ρ s−s (CV=0.15, ρ s−s = 0.6), and a range of values for ρ m−d . The grey and light grey areas display 1 and 2 standard deviations from ∆ = 0, respectively ( Fig. 2A), or from ρ c−c = |ρ m−d | (Fig. 2B). These areas indicate the region where the correlations might result from a random process. V.6. Evaluation of ξ from the experimental data In Eq. 1 ξ is the noise between sister cells that is related to the difference between sisters cell cycle, ∆T ss : Thus, the variance of ∆T ss is: V ar(∆T ss ) = 2ξ 2 .
We evaluate ξ for each experiment, by computing the left hand side of Eq. 10 from the data.

V.7. Fitting method
We simulated Eqs. 1-2 for a range of the parameters: τ 0 , k and , for 1000 lineages and 50 generations, and computed the correlations on the last generation. We fitted the model to 6 features of each experiment with ∆ significantly above zero (see Table I, Fig. 2, and Significance Computation): ρ s−s , ρ m−d , ρ c−c , ξ, the mean cell cycle, and the cell cycle CV . We found all parameters that provide close features to the experimental features. We chose the parameters that provide the closest simulated ρ s−s as the best fit, as ρ s−s is the most significant measured correlation that we have, and we reckon that ρ s−s determines the noise level. For Cyanobacteria 1 T osc was measured directly for every colony. We used the average on all colonies as the measured T osc . For all other experiments we fitted the model for T osc = 24 hours.

V.8. Periodicity computation
In order to find the periodicity shown in Fig. 3A we simulated Eqs. 1-2, without noise, for 10000 generations, and for a range of τ 0 and k, with α = −0.5. For the last 100 generations we checked whether the circadian phase (which in our case, is the absolute time modulo 1, as we normalized all the parameters by T osc ) is equal to the previous generations. If the phase is equal to the phase of 1 generation before, the period is 1. If not, but it is equal to the phase of 2 generations before the period is 2, and so forth. Periods above 10 were not computed.