Temporal scaling of aging as an adaptive strategy of Escherichia coli

See allHide authors and affiliations

Science Advances  29 May 2019:
Vol. 5, no. 5, eaaw2069
DOI: 10.1126/sciadv.aaw2069


Natural selection is thought to shape the evolution of aging patterns, although how life-history trajectories orchestrate the inherently stochastic processes associated with aging is unclear. Tracking clonal growth-arrested Escherichia coli cohorts in an homogeneous environment at single-cell resolution, we demonstrate that the Gompertz law of exponential mortality characterizes bacterial lifespan distributions. By disentangling the rate of aging from age-independent components of longevity, we find that increasing cellular maintenance through the general stress pathway reduces the aging rate and rescales the lifespan distribution at the expense of growth. This trade-off between aging and growth underpins the evolutionary tuning of the general stress response pathway in adaptation to the organism’s feast-or-famine lifestyle. It is thus necessary to involve both natural selection and stochastic physiology to explain aging patterns.


The biology of aging and senescence revolves around the duality of individual frailty and population resilience. Throughout the course of normal metabolism, components of living systems such as cells, lipids, proteins, and DNA inevitably suffer from wear and tear such as free-radical damage. This constant and collective decay eventually leads to the loss of vital functions and collapse of individuals, measured as increased mortality and decreased fertility rates for nearly any living organism or even engineered systems. As an example, for humans, the Gompertz empirical law of mortality, one of the earliest quantitative laws in life sciences published in 1825 (1), demonstrates that human death probability doubles every 8 years. While this functional senescence approach is broadly applicable to nearly all living organisms and even engineered systems, identifying the causal molecular processes for senescence of any particular organism remains a major challenge for the aging field. Beyond that, understanding how system-level failures emerge out of microscopic damage could even reveal general organization principles that give rise to self-maintaining individuals in the first place (2, 3).

Yet, for some organisms in the tree of life, aging does not lead to increased mortality and declined fertility. For organisms such as hydras, netleaf oaks, desert tortoises, and naked mole rats, mortality rates remain constant or may even decrease with age (4). More generally, the “immortal germline” theory (5) exemplified that organisms are able to repair or replace most of the damage to their components. Evolutionary biologists attribute the apparent senescence of metazoan somas to the inadequate investment in cellular maintenance as an adaptive strategy to maximize lifetime reproductive success, due to the trade-offs between the survival and reproduction of the young on one hand and maintenance for the benefit of the old on the other (6, 7). Presumably, natural selection has to operate through the molecular “levers” of damage accumulation and/or repair to achieve such life-history optimization.

Escherichia coli, a single-cell prokaryote, has historically served as a model organism that resolved many fundamental questions in biology. E. coli cells, as their metazoan counterparts, suffer from damage to their components, which lead to cellular senescence (810). In exponential growth, these damages are quickly diluted by de novo biosynthesis, and the effects of cellular senescence mitigated by rapid and robust reproduction (11). Yet, the natural life cycle of E. coli entails a much wider range of physiological conditions than exponential growth. Most bacterial cells spend much of their lives in resource-limited growth arrested conditions, where de novo biosynthesis rates are slower (12). While E. coli undergoes substantial changes in adaptation to these conditions, once the depletion of nutrients is complete, without growth-associated dilution, the accumulation of molecular damage such as protein misfolding and oxidation eventually leads to functional senescence (fig. S1). Despite the lack of fixed separation of germ and soma cells, the ability to survive the wear and tear of cellular components during growth arrests contributes to bacterial overall fitness as much as the ability for exponential growth. Therefore, we adopted a biodemographic approach (2, 3) to understand how modulation of molecular damage repair could shape aging dynamics in growth-arrested E. coli.


Longitudinal tracking of vitality at single-cell level

We first set out to acquire high-quality longitudinal physiological and life-history data at single-cell resolution of clonal populations of E. coli. For decades, limitations of traditional cell culture methods have frustrated quantitative physiologists’ attempts to understand bacterial maintenance (13, 14), due to both media-cell interactions and cellular interaction such as cross-feeding and cannibalism. Here, we designed a novel microfluidic device with cell-dimension chambers to measure individual E. coli cell lifespan in constant and homogeneous environmental settings (Fig. 1). Cells are trapped in an array of single-cell–wide dead-ended wells, with openings to a main flow channel. Constant flow of fresh media in this main channel provides necessary nutrients and eliminates metabolic waste, cell debris, and intercellular cross talk so that the environmental conditions are maintained constant over time. Cells with appropriately expressed fluorescent markers are imaged bottom-up and appear as fluorescent spots (Fig. 1B). At the end of the experiments, 70% ethanol was injected into the device to account for every cell in the cohort (see Material and Methods and fig. S2 for time-of-death determination). This experimental system allows easy tracking of cohorts of a large number of individual bacterial cells for prolonged periods (up to 7 days).

Fig. 1 Measuring E. coli lifespan distributions at single-cell resolution.

(A) A three-dimensional (3D) model of microfluidic devices used to trap and isolate large number of single cells. Upper half of the picture with the arrows represents the main flow channel, where fresh, carbon-source–free media are supplied. The array of light blue rods represents E. coli single cells trapped in a 2D array of cell-sized dead-end chambers. (B) Fluorescence microscopy image of the microfluidic device loaded with E. coli cells. Each fluorescent dot corresponds to a single cell. Cyan, constitutively expressed fluorescent protein signal; red, PI viability staining. Horizontal scale bar, 10 µm; inset, 2× magnification. (C) Sample time-lapse PI images of mortality events, from early (top) to late (bottom), deaths. (D) A heatmap of all single-cell PI signal time series in an experimental cohort. The color values of each row correspond to the time series of each cell. Cells are sorted according to lifespan (top to bottom), from the shortest living to the longest living. See Materials and Methods for data processing.

To observe cell mortality, we included in the carbon source–free medium a red-fluorescent, DNA-binding, bacterial viability dye propidium iodide (PI), which penetrates the cells only when cellular membrane potentials are disrupted. PI staining has been established as an effective proxy of cellular death and correlates well with cell viability assessed by proliferating potential (15). We used PI at a concentration fourfold lower than the concentration that previously had been shown to have no effect on E. coli viability and growth (15). Automatic time-lapse fluorescence microscopy and fixed geometry of our devices allowed longitudinal quantification of the PI signal for every single cell in the population.

Despite being genetically clonal and environmentally controlled, cells did not share the same time or manner of death, as measured by the PI time series. The transition from life to death of a representative cohort (N = 4744) could be visualized by the boundary between dark and light in Fig. 1D. Because cells were sorted vertically according to estimated lifespan, the shape of this boundary represents the survival function of the population. Considering aging as a stochastic process of system reliability reduction, the observed time of death could be modeled as the first time to hit death as an absorbing barrier, in accordance with the classic first passage time problem for stochastic models (2, 3). In addition, we observed individual differences in death trajectories. Short-lived cells tended to have very sharp PI increases that are associated with abrupt losses of membrane integrity, while those that die late tended to suffer a type of “slow death” characterized by a gradual PI increase over the course of 10 to 15 hours. Dark strips after cell deaths indicate decaying DNA or empty wells after cell debris have been washed away.

E. coli lifespan distributions follow the Gompertz law

For both demographers and reliability engineers, the age-associated increase in death probability, also known as the hazard rate h(t), is considered to be a hallmark of aging (24). The detailed population statistics derived from our experiments are particularly useful in estimating the shape of the hazard functions over the whole lifespan. We find that E. coli lifespan distributions have, as their main feature, regimes with exponentially increasing hazard rates, h(t) ~ h(t0) eb(t−t0), where h(t0) is the initial hazard rate and b is the Gompertz aging rate, i.e., the Gompertz law of mortality (1, 16). For wild-type cells, this exponential regime spans from about 13 to 93 hours (Fig. 2A), corresponding to approximately 90% of all cell deaths (see below for the 10% cell deaths outside of the exponential regime), and ranges within at least 100-fold change in hazard rates (2 × 10–3 to 2 × 10–1 hour–1). Within this exponential regime, the doubling time of hazard rate is 9.4 ± 0.5 hours (Ntotal = 4744). In comparison, the doubling time of human mortality hazards is about 8 years (1).

Fig. 2 E. coli lifespan distribution can be characterized by exponential increases in mortality rates, i.e., the Gompertz law of mortality.

(A) Binomial estimators of the hazard rates, shown in log scale. Error bars on the y axis are the 95% CI. Cell deaths are binned into discrete time intervals, which are marked by the x-axis error bars. The time bins are chosen so that the Nelson-Aalen (N-A) CIs in (B) do not overlap with each other. Shading covers regions that deviate from the exponential regime. (B) Hazard dynamics in the phase plane, where instantaneous hazards are plotted against the cumulative hazards. Data, binning intervals, and y axis are the same as those in (A). Error bars on the x axis are 95% CIs of N-A estimator. Survivorships equivalent to the exponentials of negative cumulative hazards are shown in the top axis. The brown dashed lines in both (A) and (B) represent the same naive linear fit to the Gompertz regime, while the red straight line represents maximum likelihood parametric estimations using the Gamma-Gompertz-Makeham (GGM) model. The inset provides a zoom-out view of the whole data range, while the main figure zooms in on the first 63% of cell deaths. The vertical dashed line inside the inset approximates the end of the exponential regime.

Hazard rates not only define lifespan distributions but also can be thought of as surrogates for system vulnerability, whose dynamics through time reflect the physiological consequence of aging. In this light, the Gompertz law can be interpreted as a dynamic equation governing the aging process: dh(t)/dt = bh(t), where b is the Gompertz aging rate, which can be directly observed in the phase plane of the cumulative hazards H(t) without any free parameters (Fig. 2B). The fact that the theoretical trajectory (dashed line in Fig. 2B), including the aforementioned linear regime, lies within the confidence intervals (CIs) of almost every state coordinate lends strong statistical credence to the Gompertz law for describing our data. In addition, this dynamic interpretation of Gompertz law should also constrain stochastic models of biochemical decay (fig. S3).

The observed hazard rates deviate from the exponential regime at the beginning (t < 15 hours, 3% cell deaths) and the end (t > 93 hours, 7% cell deaths; see Fig. 2B, inset, top x axis) of the total lifespan. These observations are reminiscent of similar deviations from Gompertz law in human mortality data (17, 18). In our case, the additional mortality at the early age might result from harvesting and transferring exponentially growing cells inside our microfluidic chip. Late-life hazard decelerations are also common phenomena for both model organisms and human populations and are thought not to reflect decelerations of physiological aging on the individual level, but result from changes in composition of heterogeneous populations (18). Despite being genetically and environmentally nondistinguishable, due to the stochastic nature of cellular biochemistry, we think it is reasonable to assume that physiological heterogeneity exists in our single-cell populations for such selection effect to manifest itself. We found that an extension to the Gompertz law with four parameters, named the Gamma-Gompertz-Makeham (GGM) modelhggm(t)=λ+bs1+(β1)ebtcould fully model our data satisfactorily across the whole lifespan (see Materials and Methods for parameter explanations and Fig. 3, fig. S4, and tables S1, S2, and S3 for statistics). We find that the GGM model maintains good fit to all of our data and is flexible enough to account for various types of hazard dynamics. For this reason, the GGM model enabled us to evaluate genetic impacts on the aging process, in the next section.

Fig. 3 The GSR of E. coli modulates aging rate.

(A) A scheme representing relevant regulatory features of the GSR and, in particular, the functions of the genes rpoS and rssB. (B) Experimental and GGM model survivorship. The lifespan distributions for the wild-type (wt), ΔrpoS (lacking GSR), and ΔrssB (overexpressing the GSR) strains are measured multiple times by independent microfluidic experiments. Representing the experimental survivorship, color bands are the 95% CI of the Kaplan-Meier estimators. Colored dashed lines are GGM models whose parameters are estimated from maximum-likelihood (ML) methods. (C) Hazard rates estimated using only cell deaths within discrete time intervals (error bar markers), and GGM hazard models estimated from the whole dataset using ML methods. Error bars are similar to those in Fig. 2A. Parametric comparisons and GLM-based models (see “Experimental design and experimental variation analysis” section) indicate that ΔrssB prolongs longevity entirely by reducing aging rate and that ΔrpoS both accelerates aging and increases basal mortality rate. (D) Aging rates for each strain, estimated by three independent experiments and GLMs. Error bars represent 95% CI.

Modulation of aging rate by the general stress response

Several nutrient-sensing pathways in both bacteria and metazoans have been shown to control stress resistance and reduce mortality (1921). Yet, their roles in delaying aging-related damage accumulations are often controversial. After all, it is relatively easy to prolong lifespan and reduce unnecessary mortality in an age-independent way in the laboratory by changing environmental or genetic factors unrelated to aging. Disentangling effects on aging from age-independent components of longevity (22) requires large cohort size and detailed mortality statistical resolution as described in the previous section. Both the Gompertz regime in hazard curves and the fully parametric GGM model can be used to quantify the age-dependent effects as the aging rate b, which is the exponential rate of increase of mortality risk. It can be estimated by examining the slope of the hazard curves on a semi-log scale (Fig. 2A). For example, the original Gompertz law postulation that the risk of death for humans doubles every 8 years (1) is equivalent to bhuman=ln 28 years=0.087 year1.

Many bacteria species regulate the level of cellular maintenance through a genetic pathway called general stress response (GSR), controlled by the master transcriptional regulator rpoS that is activated by nutrient deprivation among other signals (Fig. 3A) (23). Evolutionary theories of aging predicted the existence of aging-modulating mechanisms and their activation by nutrient limitation (24). To empirically test this idea, we measured and compared the full lifespan distributions of two GSR mutants, ΔrpoS and ΔrssB, with that of the wild-type strain. ΔrpoS is devoid of GSR, while ΔrssB displays an elevated GSR due to increased RpoS stability (25).

Consistent with the well-documented role of GSR in stress resistance (23), we observed that a higher GSR level promotes longevity in the microfluidic experiments, whereas the absence of GSR results in shortened longevity (Fig. 3B and fig. S5A). In contrast to similar studies, where only the average lifespan or the survival curve was measured, our microfluidic experiments can measure the hazard dynamics (Fig. 3C) of each strain and extract the rates of aging and not just longevity.

Remarkably, we found that increased GSR reduces the rate of aging (Fig. 3D and fig. S5C). Given how well the GGM model fitted for all three strains, the impact of genotypes on aging parameters such as aging rates could be extracted using generalized linear models (GLMs). With enhanced GSR, ΔrssB cells double their mortality risk every 14.1 hours, with a 95% CI ranging from 13.1 to 15.3 hours, compared to 9.4 hours for the wild-type (CI, 8.9 to 9.9 hours) and 7.3 hours for the null strain ΔrpoS (CI, 6.8 to 7.7 hours). We assessed the variability of aging rate measurements using three independent experimental cohorts for each strain (fig. S4). We visualized nonparametrically the overall experimental variations (fig. S5, A and B), and also tested parametric differences in aging rates among experimental repeats using GLMs and the Akaike information criterion (AIC). We confirmed that experimental repeats shared similar aging rates (H0: Same aging rate for experimental repeats. Difference in degrees of freedom: H0H1 = −2; ΔAICrssB = −2.56, ∆ΔAICwt = −2.26, ΔAICrpoS = −2.30; Nwt = 6867, NrssB = 6969, NrpoS = 4793; for details, see tables S1 to S4).

The systematic increase in vulnerability of E. coli in our experiments is likely driven by the catabolism of preexisting macromolecules and dissimilation of biomass (26), which is necessary to provide energy to express housekeeping genes and maintain physiological homeostasis (12, 27). Our finding that GSR modulates aging rate suggests that optimizing this maintenance energy requirement is likely one of the physiological functions of the RpoS regulon (see fig. S3).

Evolutionary trade-offs mediated by the GSR

Is slower aging an adaptive life-history trait? This was one of original questions raised by Sir Peter Medawar (28) that motivated much of later aging research. Optimal life-history theory suggested that trade-offs and constraints among fitness components shape metazoan aging rates (7). The possibility of modulating aging rate through GSR offers us the opportunity to test these ideas in a fast-evolving organism as E. coli.

The relatively well-understood GSR pathway provides a clear molecular mechanism for a trade-off between growth and maintenance. The master regulator rpoS encodes the RNA polymerase (RNAP) sigma subunit σS, which competes with the other sigma factors including the vegetative σD to recruit the core RNAP and direct the transcription and translation machinery toward the RpoS regulon (Fig. 4A). By titrating protein synthesis activity away from metabolic and ribosomal genes controlled by σD, RpoS activity inhibits growth and nutrient assimilation (29). We measured quantitatively the growth impact of modulating GSR levels and modeled its effect using a simple coarse-grained model of proteome sectors (30). We found that the proportion of protein synthesis devoted to the RpoS regulon increases the time scale of growth (top axis in Fig. 4B).

Fig. 4 Trade-offs between growth and maintenance mediated by the rpoS pathway and its fitness consequences.

(A) Scheme of the ecological processes (dashed line) and regulatory relationships (solid line) involved in the trade-offs mediated by rpoS. (B) Experimentally measured fitness of ΔrpoS (red), wild-type (blue), and ΔrssB (green) strains as functions of time spent in feast (top) and famine (bottom) conditions. Fitness is defined as the logarithmic change of population sizes (for measurements, see Materials and Methods). Ranges of environmental conditions could be modeled by two transition rates between feast and famine, k1 and k2. (C) Fitness simulation results identifying environmental regimes favoring faster and slower aging strategies, exemplified by the three strains measured in (B). The color-coded regions identify environmental conditions under which one strain dominates over the other two. Environmental conditions are parameterized by lifestyle ratio k1/k2, the ratio between time spent in famine and feast conditions, and the average duration of famine episodes 1/k2 (see Materials and Methods for simulation details).

To assess the life-history optimality of different aging rates, we integrated numerically the experimentally derived growth and mortality rates into fitness, defined as the long-term population growth rates. In contrast to metazoans, for whom fertility and mortality schedules are connected through fixed age structures, our model of E. coli life history consisted of alternating environmental episodes of varying durations, in which E. coli populations either grow (“feast”) or decline (“famine”) (Fig. 4B). These feast-or-famine cycles were parameterized by two transition rates, controlling respectively the average lengths of feast and famine episodes. To identify selective pressure for aging rates, we directly compared the fitness of the three strains with different GSR phenotypes, each representing a different strategic position in the growth-maintenance trade-off. We identified the environmental regimes selecting for faster or slower aging. We characterized the boundaries between these regimes along two axes (Fig. 4C): lifestyle ratio, defined as the ratio of time spent in famine over feast; versus average duration of famine episodes.

There are two necessary conditions for selecting slower aging strategies represented by ΔrssB. First, populations have to spend much more time in famine rather than feast so that, over the long term, populations decline rather than grow. Second, but no less important, given the same lifestyle ratio, famines should consist of longer episodes rather than short but more frequent ones. This condition is necessary due to the exponential mortality dynamics described by Gompertz law: Investing in cellular maintenance only becomes beneficial at old age when the exponentially increasing benefits of slower aging eventually overcome the more immediate cost on growth. It is the typical time scale of famine that provides the selective pressure for aging rates.

By representing the complex regulations of GSR with two mutants, we can now understand the ecological role of GSR activation and its adaptive consequences. Activated by declining nutrient availability, GSR directs resources toward internal maintenance to counter the adverse conditions, whose lengths determine the optimal activation level. Previous observations from experimental evolution of E. coli support our predictions. In continuous cultures of E. coli where the populations do not pass through prolonged growth-arresting bottlenecks, mutations that attenuate or knock out RpoS activity are among the first to arise (31). In contrast to the isolated populations under constant environmental conditions in our chip experiments, E. coli populations, in nature, influence their environments and also interact with each other. These interactions may give rise to frequency-dependent selection and evolutionary game dynamics between slower and faster aging strategies, as is observed in experimental evolution (32).


We demonstrated that aging of growth-arrested E. coli follows the Gompertz law of mortality. Moreover, by regulating the GSR, E. coli could modulate its aging rate and rescale the lifespan distribution temporally. We further articulated, in a demographic model, the trade-offs and selective pressure driving the evolution of aging rates.

To rescale the lifespan distribution without changing the shape of hazard functions, GSR has to orchestrate a coordinated response to manage various damages that accumulate on different time scales (2). RpoS regulates hundreds of genes conferring resistance to various internal and external stresses (23). For this reason, the aging rate, in our case, reflects the general level of macromolecular catabolism. We hypothesize that the maintenance energy (14), i.e., the energetic cost of homeostasis, is paid by the loss of biomass, leading to the gradual increase in the probability of death. It is thus the rate of living (33) that correlates with the aging rate in our system (fig. S3).

We provide rare experimental support for antagonistic pleiotropy (6), a key concept in evolutionary theories of aging, despite bacteria’s vast differences from animals. Sigma factor competition, a well-understood mechanism in the physiology of bacterial stress response, provides a molecular basis for a trade-off between growth and maintenance. Activated by nutritional deprivation and now shown to decrease the aging rate and prolong lifespan, GSR has immediate analogies in calorie restriction–induced longevity. Whereas many metazoan studies have only shown that calorie restriction promotes longevity, our demographic data indicate that GSR changes the rate of aging in addition to longevity promotion. It provides not only resistance against stress but also insurance against prolonged growth arrests, or more generally, protection against withering through time.

In our work, the proximate and ultimate causes of aging are integrated and applied to the iconic model organism, E. coli. The proximate cause explains aging as the stochastic and inevitable erosion of organismal order and biochemical redundancy. The ultimate cause views aging as a component of the optimized life-history strategy. In bacteria, aging and mortality processes are indeed stochastic but still could be characterized quantitatively by the Gompertz law of mortality in terms of lifespan distributions. Regulatory and selective forces do not modulate age-specific mortalities independently, as traditionally modeled (34), they rescale the whole distributions. Our empirical observations and integrated perspective in bacteria call for similar approaches in understanding aging as a general phenomenon in living systems.


Experimental methods

Microfluidic chip fabrication. Our polydimethylsiloxane (PDMS) lab-on-a-chip system consists of two layers, containing the flow channel and the array of cell-sized chambers. The two layers were fabricated separately using soft lithography technology and then bonded together to produce the microfluidic device (fig. S6).

To fabricate the negative master for the piece containing the flow channel, SU8 3050 photoresist (MicroChem, MA, USA) was patterned on a silicon wafer using photolithography. SU8 3050 was spin-coated on a silicon wafer at 4000 rpm for 30 s, baked at 95°C for 15 min, and then subjected to ultraviolet exposure (25 s, 10 mW/cm2). After post-exposure baking (95°C for 5 min), the master was developed using SU8 developer (MicroChem, USA), rinsed with isopropanol, and dried with filtered nitrogen.

The negative master for the layer with cell-sized chambers was fabricated using reactive ion etching (R.I.E.) technology on a silicon wafer. The mask was patterned on a silicon wafer using photoresist AZ5214 (MicroChem, USA), and a 100-nm layer of nickel was sputtered onto the substrate. A lift-off procedure was applied to remove the photoresist layer yielding the metal mask for the R.I.E. process. By adjusting the R.I.E. parameters [SF6, 4 sccm (standard cubic centimeter per minute); CHF3, 16 sccm; pressure, 10 mTorr; power, 30 W], we managed to achieve a large array of micropillars with high aspect ratio (diameter, 1.2 μm; height, 6 μm).

To form the device, PDMS mixtures (RTV615; Momentive Performance Materials Inc., Waterford, NY) were poured (flow channel) and spin-coated (array) onto the masters to a thickness of 5 mm (main channel) and 80 μm (array), respectively. Heat curing initially formed solid PDMS layers with patterned surfaces. After drilling inlets and outlets through the flow channel layer, and mounting the array layers onto cover glasses, the two layers were then bonded together using oxygen plasma (90 s, 1000 mTorr). Last, the assemblies were cured at 80°C overnight to produce the integrated microfluidic chips. On the day of use, the wetted surfaces of the PDMS chip were first activated by exposure to oxygen plasma (90 s, 250 mTorr), immediately followed by infusion of 20% (v/v) polyethylene glycol 400 (PEG400) solution to prevent bacterial adhesion and biofilm formation.

Medium preparation. All equipment used for medium preparation, sterilization, and infusion were made of nonleaching materials (glass, polytetrafluoroethylene i.e., PTFE, or similar perfluoropolymer material) to avoid contamination with trace level carbon sources from leachable plastic additives (table S5 for details). Medium was filter-sterilized (0.2 μm) to avoid volatile organic contamination during autoclaving, and glassware was sterilized by dry heat. Carbon-free minimum medium mentioned below refer to those prepared in this fashion.

Strain information. All lifespan distributions described above were measured for strains of the Keio E. coli BW25113 strain (“wild-type”) single-gene knockout collection (35). For the knockout strains, the presence and location of genomic inserts were verified by kanamycin resistance and polymerase chain reaction amplification. The GSR phenotypes of ΔrpoS and ΔrssB were verified using the catalase test, where ΔrpoS correctly showed minimal catalase activity and ΔrssB greatly increased catalase activity compared to wild-type (36). In addition, in developing the microfluidic device and validation of our method, we used an MG1655-derived E. coli strain with a chromosomally integrated cyan fluorescent protein under the P2rrnB constitutive promoter (Fig. 1B).

Cell culture and loading. Single isolated colonies of the bacterial strains E. coli wild type, ΔrpoS, and ΔrssB were grown overnight in minimal medium (1× M9 salts, 2 mM MgSO4, and 0.1 mM CaCl2) supplemented with 20% (w/v) glucose (final concentration of 0.4%). The following day, the overnight cultures were diluted 200-fold in 50 ml of fresh medium in 250-ml Erlenmeyer flasks and grown to early exponential phase [optical density at 600 nm (OD600) = 0.2]. This growth phase was chosen to guarantee that the variation of birth time among cells was less than one cell cycle and thus to minimize the uncertainty of lifespan measurement. Cells were concentrated by centrifugation (4000 rpm for 15 min at 37°C) and washed by three cycles of gentle resuspension with carbon-free minimal medium and centrifugation before injection into the microfluidic channels. Cells were then trapped into the dead-ended wells by centrifugation at 2000 rpm for 15 min at 37°C with a surface density up to 6.25 × 104 cells/mm2. The main channels were then thoroughly washed with carbon-free minimal medium.

Experimental setup and microscopy. A constant flow of carbon-free M9 minimal medium at 20 μl per hour was provided to the microchannels using a high-precision syringe pump (Harvard Apparatus PHD 2000 Programmable) and Hamilton GC-grade glass/PTFE syringes (Gastight 1000 Series). PTFE tubing was used to connect the syringes to the microfluidic chip. The medium was supplemented with 1.5% (v/v) PEG400 to prevent unspecific adherence of cells to the channels, and PI (5 μg/ml), as a fluorescent indicator of cell viability, was added. Cell viability was monitored using temperature-controlled (37°C) automatic time-lapse microscopy (Zeiss AX10, 63× oil-immersion objective, controlled with MetaMorph software). Focus was maintained by a Z-scanning maximum-contrast procedure using phase-contrast illumination so that dead ends of the microwells are imaged. Focus was automatically readjusted before each imaging cycle for every position and maintained within a Z range of 0.2 μm around the maximum-contrast Z-position. Phase-contrast and fluorescence images (PI signal excitation, 546/12 nm; emission, 605/75 nm) were acquired for every stage position once every hour for up to 150 hours.

Growth phenotypes. Growth phenotypes of aforementioned E. coli strains on selected carbon sources in minimal media were measured in the 96-well format using a TECAN Spark microplate reader. To induce the appropriate metabolic enzymes before growth curves could be measured, strains were first grown for 24 hours in minimal medium (M9, as in those used for microfluidic experiments) supplemented with the assayed carbon sources. The optimal densities of these cultures were determined and diluted into fresh medium identical to those used for the overnight cultures. The dilution ratios were chosen so that all experimental cell cultures have an initial optimal density of OD600 of 0.002. For each experimental well on the microplates, 50 μl of mineral oil was added to 100 μl of cell culture to prevent evaporation during the experiments. The microplates were then maintained at 37°C and shaken constantly in double-orbital motion at 150 rpm by the plate reader. Microplates with flat-bottomed wells were used to maximize agitation. OD600 readings were taken every 10 min. The growth phenotypes used in Fig. 4 are based on minimum medium culture supplemented with 60 mM acetate as the carbon source.

Statistical and computational methods

Image analysis. The cells in our experiments were trapped in an evenly spaced 2D grid. The fluorescence signals of every cell throughout their whole lifespan could be extracted at fixed positions, once images in the time-lapse stack were properly registered. A simple registration procedure might misidentify one cell for another because the cells were vertically imaged and looked very similar to each other if only local features were considered. To register the images based on global features, such as the presence/absence of cells at individual grid positions, we devised a two-pass, coarse-to-fine registration strategy. Cells were first identified and segmented within the images using the point spread function (37). These segmented binary image stacks containing global information were registered in a coarse pass using least-square minimization. The obtained 2D translations were applied to the original images. A second fine registration pass was executed on these preregistered original images using the Pyramid Approach (38). After registration, the salient positions of cells were detected on the Z-projected images of the whole time-lapse stacks, and fluorescence time series was extracted from these positions.

Regularized estimation of PI signal. To determine the true PI signal of each individual cell and to remove the noisy effects of focusing fluctuations, we designed and implemented a correction algorithm to the raw fluorescence intensity time series (fig. S7). Our three-step algorithm estimated the focusing noises of each imaging position and deduced them from the raw time series to arrive at the true PI signals of each cell. The first step of our algorithm took advantage of the fact that focusing fluctuations should change synchronously the intensities of all fluorescent objects within an imaging position, while the true PI signal from the cells should move independently of each other. By averaging the fold changes in intensities over all cells within a given time-lapse image stack, the focusing noise was enhanced while cell-specific signals were spread over hundreds of independent time series. In the second step, the averaged fold changes were decomposed into focusing noises and population-wide PI trends. This was possible because noise from focusing should have very quick fluctuations (small autocorrelation time on the order of imaging cycles) yet no long-term trends (focusing was maintained with a narrow Z range of 0.2 μm). We applied a total variation regularization algorithm (39) to effectively denoise the average fold changes to produce the population-wide long-term trends. In the last step, the population-wide PI trends were combined with the cell-specific signals to recover the true PI signals used to determine the times of death. See fig. S7 for the details and the effects of the algorithm.

Time-of-death determination. To establish the time of death or time of censorship, we used a thresholding method (fig. S2). A threshold specific to each cell was defined as the half-way point between peak PI signal of each cell and background fluorescence. The first time that PI signal exceeded the threshold was considered to be the time of death. Our method was designed to capture the mortality events of cells trapped at the ends of cell-sized chambers. Throughout our experiments, the focusing plane of the microscope was maintained at the bottom ends of these 6-μm-long chambers. Because the cells trapped there were in focus, their deaths would generate the biggest PI upticks in their respective PI time series. Thus, it would be most likely the deaths of these cells that take the PI signal across the half-way threshold between peak and background.

In cases where no PI signal declines were observed before the end of the experiments when alcohol was injected, we could not establish for the cells in question reliable peak intensities for the PI time series nor times of death. Only lower bounds of both quantities could be established. In the language of survival analysis, these cells were considered to be censored at the time that their PI signals cross the half-way point between the observed PI maximum and the background. Censoring here referred to a technical procedure in survival analysis, indicating that their survival up to the censoring time still contributed to the overall likelihood of the model, while their precise time of death were not considered. The average proportions of censorship for each of three strains were 0.0% (ΔrpoS), 2.7% (wild-type), and 28.4% (ΔrssB).

Nonparametric estimation of hazard functions. The conventional method of estimating cohort hazard rate is to estimate the cumulative hazard H(t)=0th(t)dt using the Nelson-Aalen (N-A) estimator H˜NA(t) (40). This is because H(t) can be estimated independent of binning, avoiding unnecessary reductions in statistical power. However, in most of our experiments, the sizes of cohorts allowed us to estimate h(t) directly, by binning mortality events within discrete time intervals (τi, τi + 1). h˜i, the hazard rate within (τi, τi + 1), could be estimated with binomial error h˜i ~ di/Yi/Δτi, where Yi is the number of individuals at risk at time τi, and di is the number of cell deaths within (τi, τi + 1).

Because of the exponential nature of Gompertz mortality, we could observe the Gompertz law nonparametrically without any free parameters. The instantaneous hazard rate h(t) and cumulative hazard rate H(t) could be related to each other by a simple linear relationship h(t) = bH(t) + h(t0). This relationship linked the instantaneous rate at any given time t, with the state of the whole population before time t [because H(t) = −ln[S(t)], where S(t) is the survival function]. Thus, when the two independent hazard estimates were plotted against each other in Fig. 2B, the x-axis coordinates h˜i were only determined by events occurred within each time interval (τi, τi + 1), while the cumulative hazard estimates H˜NA(τi) were only determined by events occurred before τi (40). Given the independence of these two estimates [H˜NA(τi), h˜i], the linear regime in Fig. 2B lent strong statistical credence to the Gompertz law.

Parametric models of hazard functions and lifespan distributions. We used the GGM model, hggm(t)=λ+bs1+(β1)ebt, with four parameters to fully model our data across the whole lifespan. In addition to the Gompertz law parameters, the aging rate b and β, we added two parameters to account for deviations from the Gompertz regime at the beginning and the end of the lifespan.

The additional mortality at the early age might result from harvesting and transferring exponentially growing cells from batch culture directly to growth-arresting conditions inside our microfluidic chip, in a way analogous to infant mortality in human mortality data. The standard way of modeling this is the Gompertz-Makeham model, by adding an age-independent component λ to the hazard function: hgm(t)=λ+bβebt, where β controls the age at which the Gompertz aging regime overtakes the age-independent Makeham term.

When hazard rates are high late in life, mortality removed notable portions of the populations from the cohort. Decelerations in hazard increases may reflect in changes in composition of heterogeneous populations, where the more fragile individuals are more likely to die and be removed from the cohort. An extension to the Gompertz law named Gamma-Gompertz, originally used to model old-age human model data (41), could satisfactorily model late life decelerations in our data, for both wild-type and mutant strains (see Fig. 3 and fig. S4 for the results). It could be understood as accounting for the compositional changes by an additional parameter, s, controlling the level of frailty heterogeneity: hgg(t)=bβebts[Sgg(t)]1/s, where Sgg(t) is the survival function of Gamma-Gompertz, and other parameters are as before.

Survival analysis. The lifespan data of each cohort were fitted parametrically using the family of lifespan distributions described above. Maximum-likelihood (ML) estimators of the model parameters and their CIs were obtained using the R package “flexsurv.” We tested the two-, three-, or four-parameter hazard models mentioned in the previous section, with the four-parameter GGM model being the most general, and chose the best among these candidate models according to the AIC at their ML parameters.

Goodness of fit of the ML models were tested using the one-sample version of sup{|Y(t)|}, a Kolmogorov-Smirnov (K-S) type statistic adapted to account for right censorship (42), which occurs in our experiments for cells whose PI peaks were observed after alcohol injection. In the one-sample case of goodness-of-fit testing, Y(t) is the normalized K-S distance between the empirical survivorship and the survival function of the fitted model, i.e., the fitting residues. The best models and their fitting residues Y(t) were plotted in fig. S4. For comparison and visual inspection, we also plotted alongside the Kaplan-Meier estimators for survivorship, N-A estimators for cumulative hazards, and their respective 95% CI.

Experimental design and experimental variation analysis. Populations within the same microfluidic channels were followed at multiple imaging positions, distributed evenly throughout the flow channel. We tested the statistical consistency and homogeneity of lifespan distributions among subpopulations at different imaging positions within the same channel (see figs. S5A and S8 for examples of these subpopulations). The PI time series data and lifespan distributions from each imaging position were visualized in the style of Fig. 1D (fig. S8A). Statistically, we tested for nonparametric differences between subpopulation lifespan distributions using a two-sample K-S statistic sup{|Y(t)|}, i.e., the supremum of empirical distribution distances |Y(t)| normalized to account for censorship (fig. S8B) (42). After passing the consistency test, the subpopulations from the same channel were merged to form the experimental cohorts.

Three independent experimental cohorts were tracked on separate dates for each strain. Two statistical methodologies were used to assess the level of experimental variability and test for significant differences in lifespan distribution and its parameters among the strains. Overall, variations in lifespan distribution were analyzed using hierarchical clustering. The two-sample K-S statistic sup{|Y(t)|} (2, 42) mentioned in the previous section was used as a distance measure in complete linkage clustering. Preferential aggregation of experimental cohorts of the same genotype within the lower branches of the clustering tree was taken as evidence that lifespan distributions of the strains were significantly different.

Because the GGM model adequately described our data, we also analyzed variability of its parameters and the sources of the said variability. GLMs with GGM as the probability distribution were built to examine the explanatory power of categorical covariates such as experimental cohorts (exp) or strain genotypes (strain). Specifically, the null hypothesis H0 that cohorts with the same genotype have the same aging rate b was tested against the alternative hypothesis H1 that experiment cohorts all have significantly different aging rates. Since H0 is a special case of H1, both ΔAIC and likelihood ratio test were used, and H1 was rejected for all three strains. In addition, we also tested whether the three strains have significant differences in aging rates. See tables S1 to S4 for more details on model selection and hypothesis testing.

Fitness models. Fitness is defined over a time period from t1 to t2 as f(s, t1, t2) = Δln (population size)/(t2t1), where s denotes one of the three strains. We calculated fitness based on time spent in feast or famine episodes. These fitness functions were extrapolated from experimental data. For famine episodes, changes in logarithmic population size were simply the negatives of cumulative hazards H(t1) − H(t2), so that ffamine(s, t1, t2) = [Hs(t1) – Hs(t2)]/(t2t1). Extrapolations from experimental data were done using the fitted GGM models. For feast episodes, ffeast(s, t1, t2) = 1t2t1t1t2μs(t)dt. Specific growth rates μs(t) and maximum growth rates μsmax were determined from the experimental growth curves smoothed by cubic B-splines. For feast episodes longer than the time at which μsmax was reached, we assumed exponential growth at μsmax.

Overall fitness f(s, k1, k2) of each strain s at a given set of ecological rates k1 and k2 was calculated by averaging the episodic fitness ffamine and ffeast over 5000 episodes of both feast τe,i and famine τm,i, where i = 1, …, 5000. Episode lengths τe,i and τm,i were independently drawn from exponential distributions parameterized by the two ecological rates, k1 and k2, as shown in Fig. 4B. In summary, the formula for over fitness is f(s, k1, k2) = 0τe,iμs(t)dtΔ0τm,iHs(t)iτe,i+iτm,i.

To generate the color map in Fig. 4C for visual comparison of fitness, for each set of environmental parameters k1 and k2, the fitness of the three strains f(s, k1, k2) were converted into RGB color values using the softmax function: eTfsseTfs, where s = ΔrpoS, wild-type, or ΔrssB. The parameter T was immaterial to the model conclusions as it did not affect the boundaries between different environmental regimes. It was used for visualizing the sharpness of transitions around these boundaries (for Fig. 4C, in particular, T = 200 hours).

Programming codes and data availability. The image analysis procedure described above was implemented in Java as an ImageJ plugin. Fitness modeling, plotting, and time series analyses including noise correction and time-of-death determination were implemented in Python2. We relied on R package “flexsurv” (v1.1) for the core algorithm of survival analysis and GLM, and Python package “rpy2” (v2.8.6) was used for interoperability between the Python and R codes and custom Python codes for other statistical analysis. All source codes and data are made available through GitHub (


Supplementary material for this article is available at

Fig. S1. The feast-and-famine life cycle of E. coli.

Fig. S2. Time-of-death determination from PI fluorescence time series.

Fig. S3. Gompertz law of mortality informs and constraints biochemical models of aging and mortality in growth-arrested E. coli cells.

Fig. S4. Mortality statistics and goodness-of-fit of three independent experimental replicates for the three strains under study.

Fig. S5. Experimental repeatability and variances.

Fig. S6. Fabrication and implementation of the two-layer microfluidic chip.

Fig. S7. Noise removal from focusing fluctuations using regularized estimation of PI signal.

Fig. S8. Visualization and nonparametric statistics of variability among imaging positions within one experimental cohort.

Table S1. AIC-based model selection among GLMs for the three experimental replication cohorts (i = 1, 2, 3) of wild-type strain.

Table S2. AIC-based model selection among GLMs for the three experimental replication cohorts (i = 1, 2, 3) of ΔrssB strain.

Table S3. AIC-based model selection among GLMs for the three experimental replication cohorts (i = 1, 2, 3) of ΔrpoS strain.

Table S4. Testing for aging rate differences among the three strains using AIC and GLM.

Table S5. List of plasticwares that leach carbon-supplying contaminants and their nonleaching replacements.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank X. Song for contributions on image analysis, R. Milo for insightful comments, D. Misevic and other members of INSERM U1001 for continuous advice and discussions. Funding: This study was supported by the AXA Research Fund Chair on Longevity. Author contributions: Conceived and designed the project: Y.Y., L.X., F.T., and A.B.L.; conceived and designed the experiments: Y.Y., L.X., and A.B.L.; performed the experiments: Y.Y., A.L.S., and C.L.; data analysis: Y.Y.; wrote the paper: Y.Y. with contribution from A.B.L. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All source codes and data are made available through GitHub ( Additional data related to this paper may be requested from the corresponding authors.

Stay Connected to Science Advances

Navigate This Article