Quantitation and Comparison of Phenotypic Heterogeneity Among Single Cells of Monoclonal Microbial Populations

Phenotypic heterogeneity within microbial populations arises even when the cells are exposed to putatively constant and homogeneous conditions. The outcome of this phenomenon can affect the whole function of the population, resulting in, for example, new “adapted” metabolic strategies and impacting its fitness at given environmental conditions. Accounting for phenotypic heterogeneity becomes thus necessary, due to its relevance in medical and applied microbiology as well as in environmental processes. Still, a comprehensive evaluation of this phenomenon requires a common and unique method of quantitation, which allows for the comparison between different studies carried out with different approaches. Consequently, in this study, two widely applicable indices for quantitation of heterogeneity were developed. The heterogeneity coefficient (HC) is valid when the population follows unimodal activity, while the differentiation tendency index (DTI) accounts for heterogeneity implying outbreak of subpopulations and multimodal activity. We demonstrated the applicability of HC and DTI for heterogeneity quantitation on stable isotope probing with nanoscale secondary ion mass spectrometry (SIP–nanoSIMS), flow cytometry, and optical microscopy datasets. The HC was found to provide a more accurate and precise measure of heterogeneity, being at the same time consistent with the coefficient of variation (CV) applied so far. The DTI is able to describe the differentiation in single-cell activity within monoclonal populations resolving subpopulations with low cell abundance, individual cells with similar phenotypic features (e.g., isotopic content close to natural abundance, as detected with nanoSIMS). The developed quantitation approach allows for a better understanding on the impact and the implications of phenotypic heterogeneity in environmental, medical and applied microbiology, microbial ecology, cell biology, and biotechnology.

Phenotypic heterogeneity within microbial populations arises even when the cells are exposed to putatively constant and homogeneous conditions. The outcome of this phenomenon can affect the whole function of the population, resulting in, for example, new "adapted" metabolic strategies and impacting its fitness at given environmental conditions. Accounting for phenotypic heterogeneity becomes thus necessary, due to its relevance in medical and applied microbiology as well as in environmental processes. Still, a comprehensive evaluation of this phenomenon requires a common and unique method of quantitation, which allows for the comparison between different studies carried out with different approaches. Consequently, in this study, two widely applicable indices for quantitation of heterogeneity were developed. The heterogeneity coefficient (HC) is valid when the population follows unimodal activity, while the differentiation tendency index (DTI) accounts for heterogeneity implying outbreak of subpopulations and multimodal activity. We demonstrated the applicability of HC and DTI for heterogeneity quantitation on stable isotope probing with nanoscale secondary ion mass spectrometry (SIP-nanoSIMS), flow cytometry, and optical microscopy datasets. The HC was found to provide a more accurate and precise measure of heterogeneity, being at the same time consistent with the coefficient of variation (CV) applied so far. The DTI is able to describe the differentiation in single-cell activity within monoclonal populations resolving subpopulations with low cell abundance, individual cells with similar phenotypic features (e.g., isotopic content close to natural abundance, as detected with nanoSIMS). The developed quantitation approach allows for a better understanding on the impact and the implications of phenotypic heterogeneity in environmental, medical and applied microbiology, microbial ecology, cell biology, and biotechnology.

INTRODUCTION
Ecosystems and microbial communities have traditionally been studied as an assembly of different monoclonal populations. Each population has been considered to be composed of genetically identical cells performing the same metabolic function. However, techniques with single-cell resolution were developed over the past decades, allowing us to narrow down the vision at the level of individual cells (Brehm-Stecher and Johnson, 2004;Pumphrey et al., 2009;Vasdekis and Stephanopoulos, 2015;Gao et al., 2016). This facilitates the understanding of monoclonal population's physiology, strengthening the concept that genetically identical cells can show cell-to-cell variability even when sharing the same environmental and nutritional conditions and phenomenon, generically known as phenotypic heterogeneity (Avery, 2006;Ackermann, 2015;Davis Kimberly and Isberg Ralph, 2016). By definition, the phenotype is every observable feature of an organism originating from the complex interaction between genotype, cellular biochemical mechanisms, and environment. Although a univocal definition of phenotypic heterogeneity is currently missing in literature, this usually embraces all intrinsic and extrinsic cellular noise (Simpson et al., 2009) as well as environment-induced differences arising between single cells belonging to a monoclonal isogenic population under homogeneous environmental conditions. Phenotypic heterogeneity arises from stochastic effects in random molecular and biochemical processes (Elowitz et al., 2002;Kaern et al., 2005;Kussell and Leibler, 2005;Kiviet et al., 2014), inequalities in gene expression, and transcriptional regulatory networks (Newman et al., 2006;Maheshri and O'shea, 2007;Fraser and Kaern, 2009) as well as constitutive single-cell features like cell cycle or cell aging (Sumner and Avery, 2002;Levy et al., 2012). All these effects contribute to the so-called cellular noise (Avery, 2006;Samoilov et al., 2006;Tsimring, 2014), which is extensively used in literature referring to cell-to-cell variability due to molecular mechanisms, regulatory pathways, and genetic cues (Balázsi et al., 2011;Levchenko and Nemenman, 2014). The unequal segregation of DNA hyper-structures transmitted to the daughter cells during division and the variation in populations' growth rate lead to the diversification in substrate assimilation by single cells in batch and chemostat cultures (Kopf Sebastian et al., 2015;Gangwe Nana et al., 2018). Metabolic/functional diversifications help microbial populations to cope better with stresses and Abbreviations: SIP-nanoSIMS, stable isotope probing combined with nanoscale secondary ion mass spectroscopy; K A , fraction of an element incorporated by a single cell during the incubation with isotope-labeled growth substrates; CV, coefficient of variation (equal to the ratio of standard deviation to the mean); CQD, coefficient of quartile deviation; COD, coefficient of dispersion; Q n , quantile with percentile n; Sk, skewness, indication of asymmetry in the distribution of a population; DW, distribution width on the histogram representing the distribution of a population; HC, heterogeneity coefficient developed in the present study; CSE, counting statistics error; HC corr , heterogeneity coefficient corrected for counting statistics error coming from nanoSIMS measurement; CSV R , CSV KA , counting statistics variations of isotope ratios (R) and relative assimilation (K A ); CSH KA , counting statistics heterogeneity of anabolic activity; DTI, differentiation tendency index defined as a slope (s) in a rank-activity distribution; CDTI, cumulative differentiation tendency index; SDT, strong differentiation tendency; WDT, weak differentiation tendency. fluctuations in their surrounding environment (Avery, 2006;Ackermann, 2015;Bódi et al., 2017). Indeed, differences at the metabolic level occur when cells respond to niche perturbations with different metabolic strategies gaining new ecological functions from which the whole population benefits (West and Cooper, 2016;Bódi et al., 2017). The environment affects considerably the metabolic fluxes inside bacterial cells via activation of different regulatory pathways of the central carbon metabolism, thus leading to different phenotypes within an isogenic population (Kotte et al., 2014). Yeast cells reveal multiple metabolic states under certain environmental conditions, regulating dynamically the glycolysis pathway and thus preventing metabolic imbalances (Van Heerden et al., 2014). The occurrence of cell-to-cell differences in metabolic traits results from an interplay of ecological, intracellular and extracellular factors, usually referred to as metabolic heterogeneity (Takhaveev and Heinemann, 2018). Thus, the term phenotypic heterogeneity encompasses all aspects of cell-to-cell variability including epigenetic, regulatory, metabolic, physical, or physiological aspects observed at the single-cell level.
A challenge in the resolution of single-cell phenotypic features as well as in the tracking of intracellular metabolic fluxes remains the detection limits of the applied experimental techniques and their throughput (Takhaveev and Heinemann, 2018). Nanoscale secondary ion mass spectrometry (nanoSIMS) provides high lateral resolution for single microbial cell imaging together with a mass resolving power (MRP) sufficient for tracking the assimilation of isotope-labeled substrates (Lechene et al., 2006). A combination of stable isotope probing (SIP) with nanoSIMS (SIP-nanoSIMS) has been applied to shed light on single-cell metabolic activity within microbial populations in environmental setups as well as under laboratory conditions (Musat et al., 2008;Pumphrey et al., 2009;Pett-Ridge and Weber, 2012;Sheik et al., 2015;Zimmermann et al., 2015;Jiang et al., 2016;Schreiber et al., 2016;Nikolic et al., 2017;Nuñez et al., 2017). Studies applying SIP-nanoSIMS have been focused on assimilation activity in the frame of specific metabolic functions, considering exclusively the heterogeneity in terms of isotopic enrichment at the singlecell level. Metabolic heterogeneity in monoclonal population has been shown to increase under limitation of nutrients and/or electron donors (Zimmermann et al., 2015(Zimmermann et al., , 2018Schreiber et al., 2016); nutritional and temperature upshifts as well as carbon-source competition enhance metabolic specialization among monoclonal cells (Sheik et al., 2015;Nikolic et al., 2017;Simşek and Kim, 2018). The above-mentioned studies have quantified the uptake of isotope-labeled substrates normalizing it with output from bulk measurements and applying empirical assimilation expressions. It is difficult to compare the results on metabolic activity from different studies, even when derived from SIP-nanoSIMS experiments, if a unified approach for quantitation of single-cell assimilation is not followed. A SIP-nanoSIMS-based approach has recently introduced the "relative assimilation" (K A ) as a measure to quantify the single-cell assimilation activity (Stryhanyuk et al., 2018). However, being exceptionally powerful in quantitation of isotopic composition at the single-cell level, nanoSIMS-based approaches require a long time to acquire chemical maps for a limited number of single cells (∼1-10 2 ), implying thus a considerable limitation in the method throughput.
Flow cytometry and fluorescence time-lapse microscopy have also been applied to study phenotypic heterogeneity (Davey and Kell, 1996;Müller and Babel, 2003;Balaban et al., 2004;Nikolic et al., 2013Nikolic et al., , 2017Lieder et al., 2014;Pratt et al., 2019). The output of these methods' application comprises multidimensional datasets providing information on, for example, gene expression, physiology, metabolism, and morphology of individual cells (Davey and Kell, 1996;Müller and Babel, 2003;Balaban et al., 2004;Nikolic et al., 2013Nikolic et al., , 2017Lieder et al., 2014;Pratt et al., 2019). Notably, the data are acquired within a short time while analyzing thousands of cells (∼10 3 -10 5 ). With single-cell resolution and high sensitivity, flow cytometry is however not an imaging technique, and it cannot be applied to resolve intracellular compartments and single cells in spatial arrangements.
Physical properties of single cells (e.g., their morphology, intracellular and extracellular features, and viability) and the 3-D structure of their arrangement can be characterized with optical, atomic force, or electron/ion probe microscopy techniques (Hawkes and Spence, 2007;Hlawacek and Gölzhäuser, 2016). The microscopy imaging delivers information on single-cell phenotypes (geometry, size, volume, etc.), showing that the size of individual cells differs from the population average, which in turn is modulated by the environment and the growth conditions (Grover and Woldringh, 2001;Taheri-Araghi et al., 2015;Nikolic et al., 2017;Westfall and Levin, 2018).
To our knowledge, a general approach to quantify phenotypic heterogeneity has not been developed to date. Nowadays, the interest on phenotypic heterogeneity is increasing due to its implications in medical microbiology, microbial ecology, and biotechnology. Phenotypic heterogeneity has been shown to have an impact on medical care, concerning antibiotic resistance, treatment persistence, and biofilm formation (Turner et al., 2000;Sumner and Avery, 2002;Balaban et al., 2004;Grote et al., 2015;Dhar et al., 2016;Sadiq et al., 2017;Van Den Bergh et al., 2017) but also in biomedical research for drug discovery, cancer therapy, and diagnostics due to the variable effectiveness of care treatments on different phenotypes (Almendro et al., 2013;Gough et al., 2017). Understanding cell-to-cell heterogeneity in bioprocessing is currently a big challenge due to its implications in processes' performance and therefore product yield in largescale production (Delvigne et al., 2014(Delvigne et al., , 2017. Hence, the quantitation of heterogeneity becomes relevant for an efficient process optimization in every field of biotechnology. Phenotypic heterogeneity has been so far expressed with the coefficient of variation (CV) for its comparison under different experimental conditions (Grover and Woldringh, 2001;Kopf Sebastian et al., 2015;Schreiber et al., 2016;Nikolic et al., 2017). The CV is calculated as the ratio of standard deviation over the arithmetic mean value, and its application as a measure of heterogeneity implies therefore a normal distribution of single cells in their activity or function. However, this is not always the case for a heterogeneous population, in which a large dispersion of data from the centroid is observed and the distribution shape is often asymmetrical (skewed). Measures of dispersion, such as coefficient of quartile deviation (CQD) (Bonett, 2006) or coefficient of dispersion (COD) (Bonett and Seier, 2006), are likewise used to calculate the scattering of the data around the centroids, but they are affected differently by fluctuations (outliers) of observations in a dataset. Additionally, there are currently no indices accounting for one of the distinctive features of phenotypic heterogeneity, that is, the uprising of clonal subpopulations (multimodality phenomenon) with different metabolic activities (Arnoldini et al., 2014;Li et al., 2018;Simşek and Kim, 2018). An attempt to quantify the growth-related heterogeneity of microbial populations has been recently done in a flow cytometry experiment (Heins et al., 2019), where authors considered the slope of cumulative distribution function in combination with skewness, peak width, and CV. The suggested combination of distribution parameters may represent an average or describe a most abundant subpopulation but does not provide a robust measure of heterogeneity for a multimodal distribution.
By extending the concept of CV and considering the distribution of single-cell anabolic activity measured in K A (Stryhanyuk et al., 2018), we developed a novel approach for heterogeneity expression which returns the "heterogeneity coefficient" (HC) involving the correction for the counting statistics error (CSE) coming from nanoSIMS analysis. Since neither the CV nor the HC accounts for the outcome of subpopulations, the challenge of quantitating heterogeneity in multimodal single-cell activity had to be tackled. With this purpose, a modification of power Zipf 's (1935) law describing the rank-frequency distribution of words in literary texts (Voloshynovska, 2011) was applied. The slope in rank-activity distribution of single cells, interpreted as differentiation tendency index (DTI), was suggested for heterogeneity quantitation. The DTI is independent of the population size, normalization procedures, and measure units, making its use general and universal.
In this study, the applicability of both developed indices, HC and DTI, was first tested on SIP-nanoSIMS datasets acquired for two different bacterial strains. The HC facilitates routine quantitation of heterogeneity for monoclonal populations following unimodality and can be easily calculated for any singlecell-resolved dataset with Supplementary Table S1 provided in the Supplementary Material. The DTI allows for recognition of subpopulations and expresses the heterogeneity showing multimodality. To show the wide applicability of the HC index and DTI, heterogeneity was also quantitated for datasets of flow cytometry and optical microscopy experiments. To our knowledge, this is the first time that such an approach was applied broadly for quantitation and comparison of phenotypic heterogeneity. Indeed, the two indices can be calculated for any numerically expressed feature at the single-cell level, since they are not bound to either any unit of measure or any particular technique.

Sample Preparation for NanoSIMS Analysis
Each sample was fixed overnight with 3% glutaraldehyde in sodium cacodylate buffer (0.2 M, pH 7.4, EM Science) at 4 • C. Filters (0.22 µm pore size, 25 mm, GTTP type, Merck) were coated with a 30 nm gold-palladium layer and mounted inside a stainless steel syringe 25 mm filter holder (Sartorius). Bacterial cells were filtered; washed twice with cacodylate buffer; incubated with 1% H 2 O 2 in cacodylate buffer, for 30 min; dehydrated via increasing concentrations of ethanol in water (30 50 70 80 90 96 100); and dried upon 20 cycles in a critical point drying device (Leica EM CPD 300a, Germany).

NanoSIMS Measurement
NanoSIMS analysis parameters were optimized in order to achieve the appropriate precision at the single-cell level [Supplementary Information (SI)]. From each filter of P. putida and P. stutzeri, a piece of 5 mm in diameter was cut and mounted on the 24-holes-holder for nanoSIMS measurement. Singlecell analysis was performed using a nanoSIMS-50L instrument (Cameca) acquiring the following seven molecular ion species (masses) 16 O − , 12 C − 2 , 13 C 12 C − , 12 C 14 N − , 13 C 14 N − , 32 S − , and 31 P 16 O − 2 . The MRP was set above 7,000. Samples areas of 100 × 100 µm 2 were pre-implanted with 100 pA of a 16 keV cesium (Cs + ) primary ion beam for 5 min before measurements. Smaller areas within 20 × 20 µm 2 field of view (FoV) were analyzed with a 4 pA primary Cs + ion beam, a nominal size of 100 µm for the entrance slit, 40 µm exit slits, and an energy slit cutting 20% of secondary ion intensity at the high-energy distribution tail. Samples were scanned with a 512 × 512-pixel raster size and a dwell time of 2 ms/pixel. In total, 30 planes were acquired, ensuring complete sample consumption in each FoV analyzed; 11 of those were accumulated and corrected for lateral drifting using the Look@NanoSIMS (LANS) software (Polerecky et al., 2012). Regions of interest (RoI) were drawn manually around each single cell using the 12 C 14 N − map as a biomass distribution template supported by the cell image acquired in secondary electron signal. Isotope ratio data were exported and further processed with OriginPro 2019 software for statistical analysis and graphing. The relative assimilation (K A ) values were calculated with the template provided in Stryhanyuk et al. (2018) and used for further numerical analysis and graphical representations.

Cell Cultivation and Preparation for Flow Cytometric Analysis
After pre-cultivation on peptone medium (Medium 1, DSMZ), P. putida KT2440 was inoculated in Erlenmeyer flasks with 50 ml of M9 leucine medium and 1 g/L of acetate as carbon and energy source and grown at 30 • C and 125 rpm. The strain grew with a µ max of 0.62 h −1 . After 24 h of incubation, acetate was added to the cultures again to reach a concentration of 1 g/l. Cell suspensions (1-3 ml) were sampled hourly and centrifuged for 10 min (3,200 × g at 4 • C). The resulting pellets were resuspended in 4 ml of formaldehyde solution [1 ml of 8% (wt/vol) formaldehyde/phosphate-buffered saline (PBS) in 3 ml of PBS] and incubated at 4 • C for 30 min. After a final centrifugation, the pellets were resuspended in 4 ml of 70% (vol/vol) ethanol/bi-distilled H 2 O and stored at −20 • C. For flow cytometric analysis, the cells were stained according to Koch et al. (2013). In short, 1 ml of the fixed sample was centrifuged, and the pellet was resuspended in 2 ml of PBS (6 mM Na 2 HPO 4 , 1.8 mM NaH 2 PO 4 , and 145 mM NaCl with bidistilled H 2 O, pH 7). The optical density (OD) of the well-mixed samples was measured (d = 0.5 cm; λ = 700 nm) and adjusted to 0.04 with PBS. After centrifugation of 1 ml of this solution for 10 min (3,200 × g at 4 • C), the supernatant was discarded. The cell pellet was resuspended in 0.5 ml of permeabilization buffer (0.1 M citric acid, 4.1 mM Tween 20, and bi-distilled H 2 O) and incubated at room temperature for 20 min. After a further centrifugation step, the supernatant was discarded, and the cells were resuspended in 1 ml of DNA dye solution for overnight staining at room temperature until cytometric measurement [0.68 µM 4 ′ ,6-diamidino-2-phenylindole (DAPI), Sigma-Aldrich, in 417 mM Na 2 HPO 4 /NaH 2 PO 4 buffer, 289 mM Na 2 HPO 4 , and 128 mM NaH 2 PO 4 with bi-distilled H 2 O, pH 7].

Flow Cytometric Analysis
The DAPI-stained P. putida KT2440 single cells were analyzed with a MoFlo Legacy cell sorter (Beckman-Coulter, Brea, California, USA) equipped with two lasers. The forward scatter (FSC) and side scatter (SSC) signals were measured using a blue laser (488 nm, 400 mW, Genesis MX488-500 STMOPS, Coherent, Santa Clara, California, USA). The FSC signal (bandpass filter 488 ± 5 nm, neutral density filter 1.9) is an optical characteristic related to cell size, whereas the SSC signal (trigger signal, band-pass filter 488 ± 5 nm, neutral density filter 1.9) is related to cell density. The DAPI fluorescence (band-pass filter 450 ± 32.5 nm) was excited by a UV laser (355 nm, 150 mW, Xcyte CY-355-150, Lumentum, Milpitas, California, USA). The resulting scatter and fluorescence signals were detected by photomultiplier tubes (Hamamatsu Photonics, models R928 and R3896; Hamamatsu City, Japan). The fluidic system was run at 56 psi (3.86 bar) with sample overpressure at a maximum of 0.3 psi (0.02 bar) and a 70 µm nozzle. The sheath fluid consisted of 10× sheath buffer ( Figure S11. The analog current signal delivered by the photomultiplier tubes was amplified and converted into a voltage signal within a range of 0-10 V in a preamplifier. This voltage signal was then amplified onto a 4decade logarithmic range. The analog-to-digital converter assigns each event in this signal to one of 1,024 intensity channels according to the event peak height.

RESULTS AND DISCUSSION
In the present study, the heterogeneity quantitation approach was developed with considerations of single-cell anabolic activity derived from nanoSIMS data. The applicability of the suggested heterogeneity indices was proved for the analysis of (i) cellular DNA content measured with flow cytometry of DAPI-stained single cells and (ii) length of single cells acquired with fluorescence microscopy.

Quantitation of Single-Cell Anabolic Activity and Evaluation of Activity Distribution
The two Pseudomonas strains were grown in batch cultures with defined starting concentrations of 13 C-labeled and unlabeled acetate. Afterwards, SIP-nanoSIMS experiments were performed to follow the biosynthetic (anabolic) activity and to compare the metabolic differences within the two monoclonal populations. Single-cell activity was expressed as relative assimilation (K A , Equation 1) (Stryhanyuk et al., 2018): where R f and R i are the final and initial cellular isotope ratios, respectively, and D gs is the fraction of heavy isotope in the growth substrate during incubation. The distribution of single cells in K A was considered for the analysis of anabolic heterogeneity within each population. The applicability of existing dispersion expressions for quantitation of heterogeneity is limited by the shape of the distribution. When in a certain population single cells fit a normal distribution in their anabolic activity, the dispersion in K Ai (i ∈ [1, n]) values around the mean (K A ), is expressed for a population of n cells as standard deviation (σ ): The 2σ range represents the distribution width (DW) comprising 68.2% of the population distributed equally (±34.1%) around K A . The probability density function of a normal distribution reveals the symmetrical bell-shaped profile with the probability maximum (mode) and the median centered at the mean value (K A ). The variability inside the population can be thus quantified by the CV calculated as In previous studies, the CV has been applied for the evaluation of population heterogeneity and biological noise (Grover and Woldringh, 2001;Bar-Even et al., 2006;Newman et al., 2006;Simpson et al., 2009;Schreiber et al., 2016;Nikolic et al., 2017). However, the biosynthetic activity of individual cells may deviate from a normal distribution even in monoclonal populations, which calls for more comprehensive approaches to evaluate heterogeneity at the single-cell level.

Quantitation of the Distribution Width for Heterogeneity Measurement
Even in actively growing populations, the single-cell anabolic activity distribution cannot always be described with a symmetric probability density function. When a distribution of single cells in anabolic activity is asymmetric, the K A is displaced from the distribution maximum (mode, Figure 1), rendering the CV (Equation 4) unsuitable. The median valuẽ is more appropriate since it represents exactly the centroid of an asymmetric distribution rather than the average of its terms (like K A does). Thus, to express a robust measure of heterogeneity, considering the distribution asymmetry, the centroid of the activity distribution derived asK A has to substitute the K A in the denominator of Equation (4). Additionally, to keep the quantitation of heterogeneity consistent with the CV, half of the distribution width (i.e., DW/2) has to be considered in the nominator, exactly as half of the 2σ range is considered in the nominator of the CV expression (Equation 4). With the fulfillment of these two objectives, that is, robustness and consistency, the HC has been developed. By means of this index, the anabolic heterogeneity of a population can be expressed with the DW in K A normalized to theK A as follows: The DW can be measured in different ways, for example, as standard deviation (σ ); median absolute deviation (MAD); and min-max, interquantile, or interquartile range. With the MAD, the deviation of a single-cell K Ai value from the median centroid (K A ) of activity distribution is expressed as In the case of normal distribution, the 2 × MAD range (MAD ≈ 0.67449 × σ ) comprises 46% of the population equally distributed (±23%) aroundK A . With the DW measured as 2 × MAD, the HC expression (Equation 6) returns the COD (Bonett and Seier, 2006): Thus, as a measure of heterogeneity, COD substantially reduces sensitivity to the distribution asymmetry. Expressing the DW in terms of quantiles Q (P) provides the possibility of tuning the sensitivity of a heterogeneity index to the distribution asymmetry. The Q (P) quantile splits the distribution into two parts comprising P% and 100 -P% of a population. A Q (P) quantile equals the K A value, which is higher than those values revealed by P% of cells in the population. The percentile P (%) varies within a range of 0-100%, where 100% corresponds to the whole population. Thus, quantile Q (50) of K A distribution equals the median valueK A . With a set of three quantiles (Q 1 , Q 2 , and Q 3 ; Figure 1), the DW can be expressed as an interquantile range: When the percentile P equals 25%, the quantiles Q 1 , Q 2 , and Q 3 become Q (25) , Q (50) , and Q (75) , respectively, and are called quartiles dividing a cell population into four quarters with an equal cell number. The following measures of DW are expressed in terms of quartiles: 25) ). The expression of DW in quartiles provides a reduced sensitivity to the distribution asymmetry, since it considers just 50% of the population within the interquartile range. To tune the asymmetry contribution into a heterogeneity index, one can expand or reduce the interquantile range by varying the percentile Pvalue. Hence, the maximum sensitivity is achieved with DW expressed with Q (0) and Q (100) as the [min-max] range. The calculation of heterogeneity indices with different definitions of DW is implemented in Supplementary Table S1. However, to render the HC (Equation 6) consistent with the CV (Equation 4), the interquantile range 1−3 = Q (84) -Q (16) was considered as DW in Equation (6) for further steps of the heterogeneity quantitation. In the case of normal distribution, the HC value equals approximately the CV (Equation 4). Indeed, with a limited number of single cells, quantiles are defined by real single-cell values, whereas the calculated mean value does not necessarily match any experimentally observed one (because it represents the average of the observed values). For a discrete distribution, the quantiles Q (16) , Q (50) , and Q (84) return single-cell activity values that are not matching exactly the mean and the borders of the 2σ range (i.e., the Q (16) and Q (84) values are approximately matching the ±34.1% points around the mean); therefore, CV and HC values cannot be exactly the same.

HC for Skewed Distributions
The heterogeneity in the anabolic activity of monoclonal population causes the asymmetry of the K A distribution profile (Figure 1). For this reason, an asymmetry measure was considered in the expression of anabolic heterogeneity. The asymmetry is revealed as a distribution skew toward higher values (positive skew) or lower values (negative skew). The skew causes the displacement of the centroid (Q 2 ) from the maximum (mode), gaining the difference = 2−3 -1−2 between the interquantile distances (Figure 1). The degree of distribution asymmetry is measured with skewness (Sk) expressed in interquantile distances as According to the interquantile range considered in the nominator of Equation (9), the difference is equal to Taking Equation (10) into account, can be expressed in terms of Sk and 1−3 as follows: To operate with positive values and thus make Equation (12) applicable for positively and negatively skewed distributions, the absolute value of skewness |Sk| was considered.
With the difference (Equation 12), the skewness was taken into account for the quantitation of heterogeneity. Expression of DW as the sum of the interquantile range 1−3 and the value allows for the enhancement of the HC sensitivity to the data points scattered over the distribution tails (Figure 1).
The expression in Equation (12) allows for rewriting of the DW (Equation 13) in terms of Sk and 1−3 as follows: With the DW expressed in terms of skewness and interquantile distance (Equation 14), the HC expression (Equation 6) can be now rewritten as The introduction of the Sk weighting factor ε into Equations (15 and 15 ′ ) as provides the possibility of adjusting the HC sensitivity to the skew by varying the DW within the 1−3 ± range with ε ∈ [−1; 1]. With ε = 1, the DW = 1−3 + (Figure 1) and Equation (16) turns to Equation (15). By setting the weighting factor as ε = −1 (Figure 1), we can narrow down the DW to 1−3 -, and the sensitivity is reduced. With ε = 0, the term ε × |Sk| is nullified, and Equation (16) turns to Equation (9); even in this case though, it has to be considered that the skewness affects the HC value since 1−3 andK A values are intrinsically influenced by the distribution skew. If the K A distribution is normal (ε × |Sk| = 0, Q 2 = K A , 1−3 = Q 3 -Q 2 = 2σ ), the HC expressed with Equation (16) becomes approximately equal to CV (Equation 4).

Correction of the HC for the CSE
In the calculation of anabolic activity (K A ) from SIP-nanoSIMS data, the changes in cellular isotopic composition upon incubation with stable isotope labels are considered at the singlecell level (Stryhanyuk et al., 2018). During nanoSIMS analysis, the quantification of cellular isotope composition is based on the isotope ion counts accumulated in imaging mode for each pixel within a single cell. Therefore, it is important to evaluate the effect of CSE of the acquired data on the derived K A values considered for quantitation of heterogeneity.
The CSE is derived as the square-root of secondary ion counts (N) (Sprawls, 1995).
Taking as example Carbon (C) isotopes, the calculation of their ratio (R) considers the counts of 13 C − and 12 C − isotope ions as and implies the propagation of CSE (Fitzsimons et al., 2000) into the Counting Statistics Variations of isotope ratios (CSV R ).
Thus, the presence of the isotope ratios (R i and R f ) in the expression of K A (Equation 1) implicates the CSE propagation (Fitzsimons et al., 2000) into the variation of calculated K A values (CSV K A ).
Consequently, the propagation of CSE introduces the term of Counting Statistics Heterogeneity (CSH K A ) as an error into the anabolic activity derived from the SIP-nanoSIMS experiment that can be expressed as: The correction for the CSE propagation is especially necessary when the HC value (Equation 16), derived for the cell anabolism (K A ), approaches the CSH K A term (Equation 19). . Two values of the total carbon ion counts are represented: N = 15 × 10 3 (thin lines) and N = 5 × 10 4 (thick lines). The initial isotope ratio R i = 0.011, corresponding to 1% of 13 C cellular abundance, and the CSV Ri = 0 are considered.
The CSV K A and the corresponding CSH K A depend on K A (Figure 2) and are defined with i) accumulated ion counts (N) for both isotopes, ii) isotope-labeled substrate content (D gs ) in the growth substrate.
Considering different CSV K A values for Q i quantiles of a certain activity distribution, the CSE propagation into the HC can be expressed more precisely with HC in the following way: (Equation 18) corresponding to Q i quantiles (i = 1, 2, 3) in an experimentally derived distribution of cellular anabolic activity. With the denotation of constant term in the HC expression (Equation 16), the partial derivatives are calculated as follows: The correction of HC for the CSE propagation was implemented via the subtraction of the CSV K A from the dispersion of experimental data in the following steps: 1) calculation of the CSV K A for K A values corresponding to Q 1 and Q 3 quantiles (i.e., CSV {Q 1 } and CSV {Q 3 }) with Equation (18) accounting for the isotope ion counts (N12 C and N13 C ) and the fraction of labeling isotope (D gs ) in the growth substrate; 2) calculation of the corrected interquantile range 1−3corr by subtracting the CSV defined by CSV {Q 1 } and CSV {Q 3 } (derived in the previous step) from the interquantile range 1−3 = Q 3 − Q 1 in the following way: 3) calculation of the HC corr with Equation (16) substituting the 1−3 with the 1−3corr in the numerator and considering the experimental median Q 2 =K A in the denominator: The HC corr error due to the CSE propagation into the HC corr calculation can be expressed in the following way: Figure 2 shows the comparison of CSH K A (solid lines) and CSV K A (dashed lines) for different experimental conditions (acquired ion counts N and concentrations of isotope-labeled growth substrate D gs ). The CSH K A and the CSV K A terms are enhanced with lower ion counts (N = 1.5 × 10 4 , thin lines) and lower D gs (10 at%, Figure 2A) as compared with the case when higher N and D gs values are considered, that is, N = 5 × 10 4 (thick lines) and D gs = 20 at% ( Figure 2B). For K A > 0.5, the steepness of the CSV K A and CSH K A profiles increases with the reduction of N and D gs . High CSH K A values at K A < 0.5 are due to the division of CSV K A by a smallK A value (Equation 19). With the K A approaching zero, the CSV K A becomes comparable with the DW and the CSH K A reaches its asymptote. These trends cause a huge error in HC calculation rendering the correction for CSE impossible and the HC value unreliable, when small changes in cellular heavy-isotope enrichment (small R f − R i difference in Equation 1; K A → 0) are considered. Consequently, the anabolic heterogeneity cannot be expressed in terms of HC when cellular isotope content is close to its natural abundance (e.g., short-time incubation or time prior incubation with isotopelabeled substrates). Therefore, in this study HC values were not provided for samples before incubation with 13 C-labeled acetate (hereafter T0).
The increase in CSV K A and CSH K A is steeper when the cellular isotope-label content (D f ) approaches the D gs value, In biological systems, the cellular isotope enrichment D f may approach but not exceed the D gs because the cells retain the original unlabeled material along cell division. The only exception is possible when metabolic reactions reveal an isotopic fractionation factor smaller than the unit ( R gs R f = α < 1) (Stryhanyuk et al., 2018). The values of K A , CSV K A , and CSH K A show a strong increase at higher D f especially when approaching their asymptote at D f = D gs . It is therefore not recommended to quantify the anabolic activity and the anabolic heterogeneity when the cellular isotope-label content D f exceeds 0.6 × D gs (Stryhanyuk et al., 2018). An increase of D f requires high amount of labeled substrate to be assimilated when D f approaches the asymptote at D gs . The D f proximity to the corresponding K A asymptote causes a strong increase in ∂K A ∂R f Equation (18) providing a higher error in K A calculation as well as an enhancement of the CSE propagation into the CSV K A . In turn, the increase in CSV K A intensifies the CSH K A term (Equation 19) of heterogeneity inaccuracy originated from the CSE.

HC Applicability to NanoSIMS Data
To test the applicability of the derived HC (Equation 16) on quantitation of anabolic heterogeneity, two Pseudomonas strains were incubated with 13 C-labeled acetate, sampled at different time points during their exponential growth and analyzed by nanoSIMS. The acquired isotope distribution maps were used to calculate the single-cell isotope content D f .
The dependence of K A on D f (Equations 1 and 22) is not linear (dotted line in Figure 3A, right Y-axis); for example, to increase the 13 C content (D f ) from 7 to 10 at%, the cells have to assimilate twice (×2.08) more carbon (higher steepness causing larger K A interval), as compared with the same (3 at%) increase in D f from 2 to 5 at% (lower steepness causing smaller K A interval). To illustrate this, the distribution of P. stutzeri cells in their 13 Cfraction and K A (Figures 3A,B) were shown. Due to this nonlinearity, the dispersion of cells (DW) in K A scale ( Figure 3B) does not reproduce the one in D f scale ( Figure 3A) although the same P. stutzeri cells are represented in both D f and K A plots (Figures 3A,B). In comparison with the D f distribution, the K A distribution appears more compressed (Figures 3A,B) in the range of small K A values ( Figure 3A, 15 min); instead it is more stretched and better reveals the intra-population metabolic diversifications between cells with higher K A values ( Figure 3A, 60 min).
To compare the heterogeneity in anabolic activity of the two strains, the K A distribution of P. putida cells ( Figure 3C) were plotted with the same integration interval (0.02) as the one used for P. stutzeri histogram plots ( Figure 3B). The histograms display different distribution in cellular anabolic activity between the two bacterial populations; however, this representation does not provide any quantitative measure of the heterogeneity in single-cell activity. For this reason, the developed HC-computation approach was applied for quantitation of heterogeneity and used afterwards for comparison between incubation time points of both strains.
As one would expect from an actively growing culture, the single-cell K A values increase over the time of incubation with the isotope-labeled substrate. Meanwhile, the diversity in cellular anabolic activity causes the broadening of K A distribution; that is, the DW varies over time. An increase in DW does not necessary indicate any rise of heterogeneity; namely, the HC value (Equation 15) is preserved when DW increases or decreases proportionally to the change inK A centroid. The scattering of single-cell anabolic activity (K A ) revealed as DW increase and skew of K A distribution, is indeed induced by the diversification in single-cell anabolism; that is, each cell incorporates a different amount of substrate at that specific time point.
The DW values were derived from the experimental data according to Equation (14). To account for the CSE, the corrected interquantile range ( 1−3corr , Equation 20) was considered in Equation (14) (instead of 1−3 ). The CSE does depend neither on the skew nor on the width of the distribution, but is defined by the CSV K A (Equation 18) depending on K A term, collected isotopeion counts (N12 C and N13 C , Equation 17) and content of heavy isotope (D gs ) in the growth substrate (Figure 2, dashed lines). Therefore, to account for CSE, the CSV K A values (CSV {Q 1 } and CSV {Q 3 }) were subtracted from the 1−3 interquantile range (Equation 20) without taking the skewness into account. The CSV K A dependence on K A corresponding to the nanoSIMS applied experimental conditions for both Pseudomonas strains (N = N12 C + N13 C = 15 × 10 3 counts, D gs = 20 at%) is shown in Figure 2B with the thin dashed line; thin solid line ( Figure 2B) (20) from a relatively small 1−3 interquantile range, as observed for P. stutzeri ( Figure 3B; Table 1). Instead P. putida revealed strong dispersion (large 1−3 ) of single-cell anabolic activity ( Figure 3C; Table 1) values exceeding considerably the CSV K A , thus making the CSEcorrection effect negligible. The reduction of CSE is possible via optimization of the SIP-nanoSIMS experimental conditions (D gs , D f , N12 C and N13 C ) in order to achieve minimal CSV K A and CSH K A . Increase of raster density and current of primary ions (PI) improve the counting statistics. However, the increase of PI current causes also the reduction in focus quality and may render the derived isotope-ratio values unreliable when the counting of major isotopes reaches saturation. Thus, a compromise in SIP-nanoSIMS conditions has to be found to deliver heterogeneity values with acceptable errors; that is, HC corr ± HC corr (SI, Supplementary Figures S1, S2).
The HC calculation developed in the present study for nanoSIMS-derived data, including CSE correction, were incorporated in the up-to-date version of the supplementary Excel template (Supplementary Table S1). The developed HC expression can be applied as a measure of heterogeneity not only of cellular anabolic activity, but also of any parameter measured at the single-cell level within a population: for example, length, volume, fluorescence yield, and gene expression (Nikolic et al., 2013;Gangwe Nana et al., 2018;Simşek and Kim, 2018;Heyse et al., 2019). The supplementary Excel template was therefore extended for HC calculation on the data acquired with other single-cell-resolved techniques.

Dispersion of Metabolic Activity in Monoclonal Subpopulations
Different studies dealing with phenotypic heterogeneity have shown that a population splits in subpopulations with different functions and/or activities (Simpson et al., 2009;Arnoldini et al., 2014;Kotte et al., 2014;Lieder et al., 2014;Li et al., 2018;Simşek and Kim, 2018). The same outcome was found in our results, in particular for P. putida strain. The histograms of K A distribution of both strains investigated here (Figures 3B,C) revealed several peaks (subpopulations) with different centroids of anabolic activity. For example, for P. putida cells after 60 min of incubation (Figure 3C), several subpopulations could be clearly resolved in the histogram (withK A at around 0. 03, 0.27, 0.48, 0.76, and 0.95). Even without considering the outer two subpopulations, the cell number in the first three subpopulations is comparable, rendering the overall DW of the entire distribution not representative for the K A dispersion of the single subpopulations. Also, the consideration of a single centroid value (K A ) renders the HC calculation inappropriate when the cell activity is distributed over several subpopulations. Thus, a reliable approach for quantitation of heterogeneity upon multimodal anabolic activity is necessary, ensuring the elucidation of (i) the activity dispersion in monoclonal subpopulations and (ii) its propagation into the entire microbial population.
Empirical approximation can be applied to study the relation between cellular properties, surrounding environments and experimentally observed development of heterogeneity in singlecell functions. Empirical analysis has been employed in allometry, for example, to describe the relation between animal body mass and its metabolic rate with the power function according to Kleiber's law (Kleiber, 1947). Another example of empirical approximation is the relation between cell dry mass and its volume described with a power function (Loferer-Krössbacher et al., 1998). Furthermore, the distribution of single-characteristic values ranked in their descending order has been also described with the power function. George K. Zipf had shown (Equation 23) (Zipf, 1935) that the frequency (f ) of word appearance in a text of natural language is inversely proportional to the rank (r) of a word in the word frequency table sorted in frequency descending order (Powers, 1998), where N is the number of words in a text.
The power function f (r; s, N) Equation (23) represents the Zipf 's law and thanks to the variable power index (s) is applicable in different fields whenever the power law is obeyed. George U. Yule had applied the analysis of rankfrequency distribution to study new genera development from monotypic genus during evolution (Yule, 1925). Power law approximation has been employed also in ecological applications to quantify the relative abundance of species in different ecosystems (Mouillot and Lepretre, 2000). The statistical properties of selfish spreading DNA repeats were also studied with a power-low approach, explaining how the high abundance of long DNA elements does not depend on the coded functions but rather on the ability of selfishly spreading in the host genome (Sheinman et al., 2016). The rank-frequency distribution of guanine/cytosine nucleotides in mitochondrial DNA has been approximated with Zipf 's law to distinguish families and genera for taxonomic classification (Rovenchak, 2018).
In the present study, we applied the power law approach (i) to measure the ability of single cells to differentiate from each other in their anabolic activity and (ii) to account for the subpopulations inside a monoclonal population. The derivation of the s index in single-cell rank-activity distribution allowed for quantitation of anabolic heterogeneity as explained below.

Differentiation Tendency in Cellular Anabolic Activity
The development of metabolic heterogeneity was resolved with the rank-activity distribution of single cells plotted in doublelogarithmic scale. The analyzed cells were ranked according to their anabolic activity (K A ) sorted in descending order; the cells with a low r show higher anabolism as compared to those at high r ranges.
The rank-activity plot of P. putida cell at 60-min time point (Figure 4, open circles) shows the distribution of single cells along multiple steps. Each step shows a specific slope describing the tendency of the cells to differentiate in their anabolic activity. This differentiation tendency of single cells causes the population heterogeneity and can be measured for each subpopulation by the steepness of the corresponding slope in the rank-activity distribution. In the present study, common trend in the activity differentiation, that is, fitting a single slope, was considered as criterion for assignment of single cells to a certain subpopulation. The K A distribution slope revealed in the double-logarithmic coordinates plot FIGURE 4 | Rank-activity distribution (bottom X-axis) of Pseudomonas putida single cells after 60 min of incubation with an isotope-labeled growth substrate. The relative assimilation of each single cell is represented in the rank-activity plot with hollow circles. For comparison, the corresponding histogram from Figure 3 was overlaid (top X-axis).
(log (K A ) ; log (r)) can be described with s exponent of a power function as: predicting the anabolic activity (K A ) of single cells with their rank r in the cells series sorted in K A descending order. The K A stays constant (K A = C ) with the power index s = 0 (Equation 24) describing the case in which cells grow without differentiation in their anabolism and their heterogeneity is approaching zero. Higher power index values s indicate a steeper slope of rank-activity distribution, implying an increase in the differentiation of anabolic activity and a heterogeneity gain. Thus, the power index s was considered as a measure of anabolic heterogeneity and is referred hereafter as DTI (s).

Quantitation of Differentiation Tendency in Anabolic Activity
Rank distributions of real experimental data, following a power law, show usually a single slope in the range of intermediate r values, but its profile in the range of low and high r values shows a considerable deviation from a single slope. Numerous modifications of the Zipfian function (Yule, 1944;Mandelbrot, 1953;Simon, 1955;Lavalette, 1996) have been suggested for the improvement of experimental data fit at either low or high r values. Poor fit accuracy and lack of meaningful parameters' interpretation limit the application of power law functions for approximation of experimental data. In the present study, the expression of the Zipfian function Equation (25), previously suggested by Voloshynovska for the analysis of word frequencies in texts (Voloshynovska, 2011), was adopted for the approximation of the rank-activity distribution of our experimental data (Figure 4): The r + q term, suggested by Mandelbrot, provides the possibility to improve the Zipfian fit in the low-rank range (small r and high K A values). Lavalette's term ( N N−r+1 ) facilitates a more accurate fit of rank-activity distribution in the high r range (high r approaching N and low K A values) and takes explicitly the size (N) of a population (or a subpopulation) into account.
In the case of cell rank-activity distribution, the s in the power exponent (Equations 24 and 25) stays invariant to linear normalization and to population size (Voloshynovska, 2011), allowing therefore for the following: i) comparison of heterogeneity in subpopulations with differentK A centroids; ii) quantitation of differentiation tendency in anabolic activity of subpopulations with low cell abundance; iii) resolution of subpopulations possessing close median values; and iv) elucidation of differentiation in anabolic activity even after a short incubation with stable-isotope-labeled growth substrates.
Importantly, the CSE propagation does not influence the slope of rank-activity distribution as it does on the DW considered in the HC expression (Equation 6). Therefore, DTI is less sensitive to the CSV K A variation (Equation 18) as long as D f is kept well below D gs and a sufficient number of isotope ion counts is accumulated during the data acquisition with nanoSIMS; otherwise, high CSV K A value smears the rank-activity distribution and may enhance the uncertainty in the slope determination, lowering the accuracy of the experimental data fit (Figure 4) considered below. For a unimodal anabolic activity, the rank-activity distribution was expected to show a single slope, namely, without a subpopulation outcome. To underpin this assumption, simulated normal distributions (Supplementary Figures S3A, S4) were approximated with the Zipfian function (Equation 25). Normal distributions were simulated, keeping the mean value fixed and changing the σ values; then DTI and CV values were calculated for each of the distribution, revealing the same trend in σ (Supplementary Figure S3B). This behavior of DTI provided proof of its suitability as a heterogeneity index under unimodal activity. Moreover, the Zipfian approximation delivers the slope values s = 0.1624 ± 0.0006 for 1,000 cells, s = 0.1631 ± 0.0011 for 200 cells, and s = 0.1619 ± 0.0082 for 50 cells, proving the invariance of s (DTI) to the population size, as shown in Supplementary Figure S5.
For the instance in which the population is split in multiple subpopulations (multimodal distribution, Figure 4), the rankactivity distribution was implemented here with the K A (r; q, s, N) expression (Equation 25) represented as a combination of subpopulation activities K Ai (i ∈ [ 1, m]):  26) provides a set of slope values s i , each of them characterizing the DTI in the anabolic activity of a subpopulation. In order to obtain a unique value, able to describe how the differentiation in anabolic activity of single subpopulations is propagated into the entire monoclonal population, the CDTI (S) was introduced. The CDTI was calculated as the distance in the m-dimensional Euclidean space considering s i values of m subpopulations as differentiation propagation vectors using the following expression: The coefficients γ i were introduced to weight the contribution of each subpopulation into the cumulative differentiation. The weighting term γ i is expressed as a logarithm of a cell number (n i ) in a subpopulation relative to the logarithm of the total number of cells in the whole population (N).
Taking the weighting term into account, the CDTI was expressed in the following way: The propagation of the fit errors s i into the calculated CDTI (Equation 27) results in S expressed as Frontiers in Microbiology | www.frontiersin.org The CDTI calculations (S and S) developed in the present study were implemented in the supplementary Excel template (Supplementary Table S2).

Power Law Fit of Rank-Activity Distribution Derived From NanoSIMS Data
The rank-activity distributions of P. putida cells after different incubation times were approximated with the multicomponent Zipfian function (Equation 26; Supplementary Figures S6A, S7). First, DTI (s i ± s i ) values were derived as fitting parameters for all the revealed subpopulations; then the CDTI (S ± S) values (Equations 27 and 28; Supplementary Table S2) were calculated and compared with the corresponding HC corr ones (Figure 5).
Even when calculated with ε ∈ [−1; 1] (Equation 21), the HC corr (Figure 5A) did reveal neither a clear correlation nor a common trend with the corresponding CDTI ( Figure 5B) for P. putida. To clarify this discrepancy, the rank-activity plot was overlaid with the corresponding histogram (Figure 4). The representation of K A distribution with a rank-activity plot allows for distinction between cells with close K A values, providing a single-cell resolution, which is not achievable with a histogram plot. The profiles of rank-activity distribution (Figure 4;  Supplementary Figures S7, S9) show the cells dispersed into two main groups. In general, the cells in the lower-rank range (high K A values) are fitting moderate slopes, whereas those in the high-rank range show considerably steeper slopes. Thus, the population is clearly split into two cell-activity trends: (i) weak differentiation tendency (WDT) and (ii) strong differentiation tendency (SDT). The SDT cells (steeper slopes) can also be observed in the lower-rank range (high K A values) as revealed by P. putida populations at time points of 90 and 120 min (Supplementary Figure S7). The SDT cells showing high anabolic activity (high K A values) are metabolically different from those showing low anabolic activity (low K A values). Therefore, high-SDT cells (high K A values) were differentiated from low-SDT cells (low K A values), and the CDTI values were calculated separately for the corresponding SDT subpopulations of P. putida at 90-and 120-min time points (Figure 5C). At a time point of 60 min, P. putida populations did not show high-SDT cells; consequently, the CDTI values were calculated for WDT and low-SDT cells. Considering the CDTI calculated separately (Figure 5C), the SDT subpopulations reproduced the descending trend of the HC corr derived with ε ≤ 0 ( Figure 5A), providing the major contribution to the heterogeneity of the entire population (CDTI, Figure 5B).
The activity of P. stutzeri single cells (Figure 3B), contributing to the main peak in the histogram, was considered to be unimodal. Those cells were ascribed to the WDT subpopulation, and their rank-activity distribution was initially approximated with the single-component Zipfian function (Equation 25), delivering DTI as s ± s values. However, this approximation provided considerably high relative errors ( s/s × 100 [%]; Supplementary Figure S8). Consequently, the rank-activity distributions of entire P. stutzeri populations were approximated (Supplementary Figures S6B, S9) with the multicomponent Zipfian function (Equation 26), and the CDTI values were derived as S ± S (Supplementary Table S2) according to Equations 27 and 28 ( Figure 6B). The low-SDT and WDT subpopulations were recognized in the rank-activity distribution, and the corresponding CDTI values were calculated separately ( Figure 6C). The increase in CDTI from 15-to 30-min time points is driven by low-SDT subpopulations (Figure 6C, triangles), which correspond to the single cells showing low K A values in the histogram (Figure 3A). Those cells are included into the HC corr calculation when ε ≥ 0 (Figure 6A), resulting in the HC trend following the one revealed by CDTI (Figure 6B) for the entire P. stutzeri population. Hence, the SDT subpopulations were shown to contribute mostly to the overall heterogeneity in anabolic activity of P. putida and P. stutzeri. The CDTI trends over time of both Pseudomonas strains are compared in Supplementary Figure S10.
As mentioned in Correction of the HC for the CSE, the calculation of HC on SIP-nanoSIMS data is unreliable for T0 when small changes in cellular heavy-isotope enrichment (small R f -R i difference in Equation (1); K A → 0) are considered and cellular heavy-isotope content is close to the natural abundance. Instead, rank-activity plot in double-logarithmic coordinates resolves also the low-abundance subpopulations (i.e., a slope can be defined with three cells already), allowing for calculation of CDTI and therefore the quantitation of heterogeneity prior to labeling or at a very short incubation time. The results of multicomponent Zipfian fit of rank-activity for both strains at T0 are presented in Supplementary Figure S6. Table S2) were introduced into the CDTI plots for both strains (Figures 5B,  6B; Supplementary Figure S10) to elucidate the development of heterogeneity, expressed as CDTI, starting from T0. Thus, at T0, both Pseudomonas populations show already anabolic heterogeneity, revealing multimodality; afterwards, each species develops distinct tendency in heterogeneity according to a strainconstitutive strategy.

Method Applicability
To demonstrate the broad applicability of the suggested heterogeneity indices, HC and DTI/CDTI were derived for two different datasets: (i) the distribution of DAPI-stained single cells of P. putida acquired with flow cytometry and (ii) the distribution of Escherichia coli single-cell length upon different growth conditions reported by Nikolic et al. (2017).

Heterogeneity Quantitation on Flow Cytometry Data
The P. putida cells were analyzed with flow cytometry at eight time points upon cultivation in batch for 26 h. The cellular DNA content was measured as DAPI fluorescence intensity and visualized in dot plots ( Figure 7A; Supplementary Figure S11) together with the cell-size-related FSC intensity. The cells were distributed within a subpopulations' pattern that changed dynamically during cultivation. The boundaries between the five subpopulations G1, G2, G3, G4, and Gx were defined with this pattern according to the local minima in the DAPI fluorescence intensity histogram ( Figure 7B) at inoculation (0 h). The distribution of cells along this histogram represents the different chromosome numbers per cell that can occur during different states of cell growth. During the cell cycle, cells increase in size (frequently, but not always) while duplicating their DNA, and consequently, their FSC intensity increases (Müller, 2007).
The observed changes in cellular DNA contents result in heterogeneous distributions that can be described by the HC index and the CDTI. For this purpose, the flow cytometric raw data in FCS 3.0 format (Seamer et al., 1997) were treated as described in Supplementary Figure S11. In short, the values in logarithmic scale were translated into equally distributed 1,024 channels. In a next step, only the range of channels that represented cells was plotted. Finally, the cells were ranked according to their DAPI fluorescence intensity, and the multicomponent Zipfian function (Equation 26) was fitted to the resulting rank-distribution curves (Figure 7). The changes in the Zipfian slope ( Figure 7C) matched the established subpopulation segregation, deduced from the dot plots ( Figure 7A) and shown by the corresponding histograms ( Figure 7B). The heterogeneity in cellular DNA content was expressed with the HC index (despite the multimodality) and the CDTI calculated (Supplementary Table S2) with the results of multicomponent Zipfian approximation. The development over the 26-h growth curve of both indices is represented in Figure 8.
For the initial and final growth time points, the HC was highly dependent on the ε value, which can be used to adjust the index sensitivity to the distribution skew. The HC decreased from the time of inoculation, when the population is constituted by cells with varying chromosome numbers, to the 1-h time point with most cells containing only one chromosome. During cell growth, the HC showed the heterogeneity boost due to an increasing number of cells harboring two, four, and more chromosomes.
Interestingly, the HC values do not change much with the variation of ε values during the exponential phase. The ε value weights the asymmetry (skew) contribution in the HC. The reduced effect of ε can be considered  as an evidence for a negligible skewness. In the case of multimodality, a small skewness value may be achieved with the centroid (median, Q 2 ), located equidistant between two quantiles (Q 1 and Q 3 ), even though a single Q 2 value is not representative in the case of a multimodal distribution. Finally, at the 24-h time point, with the population in the stationary phase, the HC values differentiate according to the chosen ε values again. In the present study, the error in graded fluorescence intensity was estimated to be 2 of 1,024 channels (<0.5%) and therefore neglected upon the HC calculations.
The HC describes the general dynamic of heterogeneity during the population growth, but its explanatory power is restricted in case of subpopulation segregation displayed by the P. putida culture. In this case, the CDTI is more suitable for the heterogeneity quantitation since it accounts for multimodality. The CDTI G1−Gx was initially calculated for the entire population with all slope values derived from the Zipfian fit (Figure 8B, rectangles). In order to avoid overweighting the influence of low-probability events in Gx that likely comprise cell aggregates (0 h, 6.60%; 1 h, 0.45%; 2 h, 0.18 %; 4 h, 4.44%; 6 h, 6.80%; 8 h, 9.49%; 24 h, 6.00%; and 26 h, 3.28%), the CDTI G1−G4 (Figure 8B, circles) was additionally calculated for the cells in the clearly segregated subpopulations G1-G4. The CDTI G1−G4 value increased from 0 to 1 h when cells started to duplicate their chromosomes after a short lag phase. At 2 h, most cells have doubled their DNA content and subsequently start the duplication resulting in (i) the majority of cells having a single chromosome and (ii) heterogeneity decrease. In the exponential growth phase, the population showed increasing CDTI G1−G4 at 4 and 6 h with higher chromosome numbers and uncoupled DNA synthesis. The confinement of cell distribution in the G1-G4 subpopulations resulted in the reduction of CDTI G1−G4 at 8 h when the culture gets into the stationary phase. In the advanced stationary phase at 24 h, the cell distribution pattern is very similar to the one upon inoculation (0-h time point). This pattern similarity is supported by the close CDTI G1−G4 values derived for 0-and 24-h time points. The restarting cell activity after the acetate feeding at 24 h was captured by the slight increase of CDTI G1−G4 at 26 h. Thus, the CDTI G1−G4 (Figure 8B, circles) describes the dynamics of the population heterogeneity, whereas the HC ( Figure 8A) represents its general trend missing the resolution of subpopulations.

Heterogeneity Quantitation on Cell Length Distribution
The heterogeneity quantitation was also applied to cell length distributions from a published dataset (Nikolic et al., 2017). In their work, E. coli cells, grown in the presence of two different carbon sources, did not specialize in a bimodal fashion but rather followed unimodality, consuming both substrates simultaneously at different rates. The highest degree of heterogeneity of the monoclonal population (measured with CV) was found during co-feeding with glucose (Glc) and arabinose (Ara) under carbon limitation as compared with the growth with one single carbon source under nitrogen limitation (Nikolic et al., 2017). This behavior was a consequence of different gene expression and transcription levels as well as different rates of singlecell growth in chemostat. Nutrient-dependent changes in the cultivation environment influence the growth rate, which in turn influences the cell size (Chien et al., 2012;Westfall and Levin, 2018). Bacterial cells have to find a compromise between the maintenance of a certain DNA amount and their cytoplasmic size in order to keep the DNA-to-cell mass ratio constant upon different growth rates. However, when external conditions start to be unfavorable, while cells continue to duplicate their DNA via multiple duplication forks, the mechanism of division via the FtsZ ring and other intracellular molecules is inhibited, thus resulting in a larger size of cells (Chien et al., 2012). The data on cell length provided by Nikolic et al. (2017), confirmed these experimental observations. This dataset was analyzed (Figure 9) to elucidate the heterogeneity in cell size upon different conditions reported in their study (Nikolic et al., 2017). The derived CV, HC index, and DTI are shown in Figure 10. The highest median value of cell length (Figure 10A,  circles) was found with the "3 µM Ara + 3 µM Glc" substrate composition, which represents the strongest carbon limitation. The HC value calculated with −1 ≤ ε ≤ 0 ( Figure 10B) followed the CV trend ( Figure 10A, rectangles), as expected for unimodal distribution. The approximation of rank-length distributions with the single-component Zipfian function delivered the slope values (s, DTI; Figure 10C), bringing two interesting outcomes: (i) the heterogeneity of population in cell size upon "3 µM Ara + 3 µM Glc" is not much higher than that upon "10 µM Ara + 10 µM Glc" as revealed with HC (ε = 0; Figure 10B) although this difference was overestimated with the CV (Figure 10A) and (ii) the length of single cells showed a unimodal distribution confirmed by small s (and high goodness of fit), which is in agreement with the unimodality reported in Nikolic et al. (2017).

Guidelines on the Heterogeneity Quantitation
The summary of the approaches (HC and DTI/CDTI) developed in this study for heterogeneity quantitation is provided in Table 2. An accurate quantitation of cellular anabolic activity requires a reliable data acquisition with single-cell resolution. The aspects described for the SIP-nanoSIMS experiment (step 1 in Table 2) are, in general, relevant for experiments implying data acquisition in counting mode.
Each step in Table 2 contains a link to the corresponding section, equation, or figure in the main text or in SI. The calculation of HC including the CSE correction is implemented in the Excel template (Supplementary Table S1). The results of nanoSIMS data processing with the LANS software together with other input parameters (like D gs , D i , and N) have to be pasted into the appropriate green-marked fields, and the Excel template calculates the outputs including K A , HC, and HC corr values. If the analyzed population of single cells follows unimodal anabolic activity, heterogeneity can be measured with HC or HC corr .
Unimodal rank-activity distribution can be also approximated with the single-component Zipfian function expressing the heterogeneity with the Zipfian slope as s ± s. The approximation of multimodal rank-activity distribution with single-component Zipfian delivers poor fit accuracy, resulting in large s values. In such a case, the rank-activity distribution has to be approximated with the multicomponent Zipfian function, and heterogeneity is measured with CDTI (S ± S). The CDTI calculation was also implemented in the corresponding Excel template (Supplementary Table S2). In this study, the OriginPro 2019 software was used for Zipfian approximation of rank-activity distributions.
When heterogeneity has to be quantitated on a set of epifluorescence microscopy data, the error in the HC and DTI/CDTI indices may arise from the following factors: (i) non-linearity of fluorescence signal detection; (ii) contribution of fluorescence background; (iii) presence of contaminants and extracellular substances. These factors may modulate the fluorescence intensity affecting the resulting HC and DTI values. Both indices are influenced by background subtraction while the DTI remains invariant to linear normalization and cell number.

CONCLUSIONS
In the present study, indices to measure metabolic heterogeneity were subjected to a comprehensive consideration. The expression of HC was developed for the quantitation of heterogeneity within a single-cell dataset. The adjustment of HC sensitivity to the distribution asymmetry was discussed and implemented in HC expression. When calculated with the same "sensitivity settings, " HC allows for heterogeneity comparison among different datasets. Importantly, the implemented HC expression returns the CV value if applied on datasets that reveal normal distribution. For the HC calculated from the nanoSIMS data, the correction for CSE was implemented. The influence of experimental conditions on CSE propagation

Heterogeneity quantitation steps
Tools applied 1. Acquisition of data with a single-cell resolution technique • Stable isotope probing (SIP) with nanoscale secondary ion mass spectrometry (nanoSIMS) (SIP-nanoSIMS) experiment (as a specific case) 1.1. SIP: the dynamic range for cellular isotope enrichment (D f ) is defined by isotope label content in the growth substrate (D gs ); relative assimilation with acceptable error (K A ± K A ) is achieved with D f ≤ 0.6 × D gs ( §3.1.3) 1.2. nanoSIMS: maximal precision of single-cell isotope content and pixel-average ion counts per cell can be achieved with the optimization of nanoSIMS conditions (SI §1, Supplementary Figures S1, S2). Generally valid for data acquisition in counting mode into the calculated HC was discussed, and appropriate recommendations for a SIP-nanoSIMS experiment setup were made. Anabolic heterogeneity of monoclonal populations can be also measured with the slope s of a rank-activity distribution, characterizing the tendency of cells to differentiate in their activity within several subpopulations. The DTI (s) is derived for each subpopulation as an exponent of the Zipfian power law approximation function. The CDTI (S) is derived from DTI values of a subpopulation set and measures the anabolic heterogeneity of the entire population. The power exponent s (DTI) is invariant to linear normalization and to the number of cells in a subpopulation/population. Both HC and DTI can be used as heterogeneity indices when a population shows unimodal anabolic activity. If the anabolic activity becomes multimodal, the DTI of each subpopulation (s i ) has to be first derived with multicomponent Zipfian approximation, and thereafter, the calculation of CDTI (S) delivers a measure of anabolic heterogeneity of the entire monoclonal population. Besides SIP-nanoSIMS, the applicability of both developed indices was shown on flow cytometry and epifluorescence microscopy datasets. The HC and DTI/CDTI were proven to be robust indices for quantitative and comparative analyses of heterogeneity in single-cell-resolved studies. Thus, these indices are not restricted to particular techniques, have a broad range of applications, and measure phenotypic heterogeneity in all its aspects: metabolic, physiological, and morphological.

DATA AVAILABILITY STATEMENT
The flow cytometry dataset generated and analyzed in this study is available in the FlowRepository (ID: FR-FCM-Z2AS).
within the Research Program Terrestrial Environment of the Helmholtz Association.