To generate large quantities of aged yeast cells, necessary to perform bulk measurements, we used a method, in which cells were immobilized on iron beads and trapped inside a column (Janssens et al., 2015). Briefly, cells were labelled with biotin and linked to streptavidin-coated iron beads. This iron bead bound cell culture was then grown in a column, equipped with an iron grid, in which the beads (and the cells attached to them) were trapped by a magnet. A continuous medium flow through the column washed out most emerging daughter cells and kept the mother cells in a constant, nutrient-rich environment. With the used flow rate of 170 mL h⁻¹, the glucose concentration stayed almost constant (only dropped from 21.7 to 20.1 g L⁻¹) and the concentration of major byproducts (pyruvate, succinate, glycerol, acetate and ethanol) never exceeded 1 g L⁻¹. Furthermore, the dissolved oxygen saturation never dropped below 75%. The precise instrumental as well as experimental setup for the column-based cultivation and harvest can be found in Janssens et al. (2015).

As samples harvested from the column still resembled a mixture of mother, daughter and dead cells and any subsequent sorting step, aiming at an absolutely pure mother cell fraction would have inherently led to a distortion of the metabolic phenotype, we opted for an approach also followed in our previous study (Janssens et al., 2015), to computationally infer the phenotype of each subpopulation. Specifically, we generated at each aging time point three samples with different proportions of mothers, daughter and dead cells (i.e. (1) from the column effluent, (2) from the column after an additional washing step, (3) from the washing solution (in the following referred to as mix 1, 2 and 3)). After harvesting and before the respective analysis (and for the physiological characterization additionally at the end of the growth experiment), the cell count specific fractional abundance of each subpopulation in each sample was determined by flow cytometry and a combined dye-staining with propidium iodide and avidin – FITC. Later the metabolite concentrations and the cellular physiologies of each individual cell population (i.e. mother, daughter and dead cells) were mathematically inferred from data originating from the mixed samples and the determined fractional abundance.

To allow the cells to recover from any possible stress during the sampling procedure, all samples were transferred in an Erlenmeyer flask containing 10 mL medium, adjusted to a cell density of 2 × 10⁷ cells mL⁻¹ and incubated for 20 min at 30°C and 300 rpm prior analysis.

A sample of 3 × 10⁷ cells was taken from the Erlenmeyer flask and immediately quenched in 10 mL −40°C methanol. The cells were separated from the organic solvent by centrifugation (5 min, 21’000 g, 4°C), washed with 2 mL −40°C methanol, separated again by centrifugation and stored at −80°C. For the following analysis, the cell pellet was re-suspended in 900 µL −40°C extraction buffer (methanol, acetonitrile and water, 4:4:2 v/v/v supplemented with 0.1 M formic acid) and an internal standard of ¹³C-labeled metabolites was added to the extraction. This standard was obtained and quantified from exponentially growing cell cultures prior to the experiment (Wahl et al., 2014). The extraction solution was agitated for 10 min at room temperature and thereafter centrifuged at maximum speed. The supernatant was transferred to a new vial and the cell pellet re-suspended in 900 µL −40°C extraction buffer and the extraction procedure was repeated a second time. The supernatants from both steps were combined and centrifuged for 45 min at 4°C and 21’000 g to remove any remaining non soluble parts. Thereafter, the supernatant was vacuum-dried at 45°C for approximately 1.5 hr and prior to the further analysis dissolved in 200 µL water.

The extracted metabolite samples were analyzed using a UHPLC-MS/MS system. The chromatographic separation was performed on a Dionex Ultimate 3000 RS UHPLC (Dionex, Germering, Germany) equipped with a Waters Acquity UPLC HSS T3 ion pair column with precolumn (dimensions: 150 × 2.1 mm, particle size: 3 μm; Waters, Milford, MA, USA). The injection volume was 10 μL and the samples were permanently cooled at 4°C. A binary solvent gradient was employed (0 min: 100% A; 5 min: 100% A 10 min: 98% A; 11 min: 91% A; 16 min: 91% A; 18 min: 75% A, 22 min: 75% A; 22 min: 0% A; 26 min: 0% A; 26 min: 100% A; 30 min: 100% A) at a flow rate of 0.35 mL min⁻¹ where solvent A was composed of 5% methanol in water v/v supplemented with 10 mM tributylamine, 15 mM acetic acid and 1 mM 3,5-heptanedione and isopropanole as solvent B. The detection was done using multiple reaction monitoring (MRM) on a MDS Sciex API365 tandem mass spectrometer, upgraded to EP10+ (Ionics, Bolton, Ontario, Canada) and equipped with a Turbo-Ionspray source (MDS Sciex, Nieuwerkerk aan den Ijssel, Netherlands) with the following source parameter: NEB (nebulizing gas, N2): 12 a.u., CUR (curtain gas, N2): 12 a.u., CAD (collision activated dissociation gas): 4 a.u., IS (ion spray voltage): −4,500 V, TEM (temperature): 500°C.

The concentrations of intracellular metabolites were determined from samples harvested after 10, 20, 44, and 68 hr. The samples were measured in six replicates and the average of this replicates was used for the mathematical inference. To validate the interference approach we independently determined the intracellular metabolite concentrations of biotin labeled cells before loading them onto the column.

The general idea of the in the following described mathematical inference rests on the concept that a system of linear equations can be solved if the number of independent equations is greater or equal than the number of unknowns. This was implemented by generating at each time point three samples (i.e. mix 1, 2 and 3, cf. Methods 1). The measured concentration in each of these three samples is constituted as the sum of the two unknown concentrations in mother and daughter cells, weighted by their respective known fractional abundance.

Specifically, the in each sample (with n^cell cells) measured amount of metabolite, n^meas, contains metabolites originating from mother (mo) and daughter (da) cells. As dead cells were considered to be lysed and their metabolite content accordingly leaked into the medium, we assumed that their contribution to the total metabolite pool can be neglected. With taking the respective volumes of mother and daughter cells (Method 5 and Figure 2—figure supplement 1), and the fractional abundance of each population into account, the amount of substance of each metabolite in each cell is given by,

 where n^meas_i,j,k is the measured amount of substance (unit mol) of the metabolite i in the sample j (i.e. mix 1, 2 or 3) at the aging time point k (i.e. 10, 20, 44 or 68 hr), n^cell_j,k the total amount of cells in the respective sample, α_j,k and β_j,k the cell count specific fractional abundance of mother and daughter cells, V^mo_k and V^da the cell volume (unit L cell⁻¹) of mother and daughter cells and c^mo_i,k and c^da_i the unknown metabolite concentration (unit M) in mother and daughter cells. Note that c^da_i and V^da are not indexed over the aging time points k, as we assumed that the daughter cell phenotype does not change over time (i.e. daughter cells produced by young mothers are identical with daughter cells produced by old mothers). To infer the intracellular metabolite concentrations c^mo and c^da from the measurements, n^meas, we formulated a non-negative least square regression problem of the form,

 where the matrix A contains all fractional volumes α_j,k V^mo_k and β_j,k V^da in every sample j at every aging time point k, the vector c the unknown concentrations c^mo_i,k and c^da_i of the metabolite i in mother and daughter cells at every aging time point k and the vector n all metabolite measurements, n^meas_i,j,k, normalized by the total amount of cells in the sample, n^cell_j,k, in every sample j at every aging time point k.

The regression problem in Equation 2 was implemented in MATLAB (Release R2013, MathWorks, Inc, Massachusetts, USA) and the unknown metabolite concentrations, c, in mother and daughter cells were identified using the function ‘lsqnonneg’. The uncertainty of the estimation was then determined by leave-one-out cross-validation, where we one-by-one removed data points from the set and repeated the estimation procedure (Figure 2—figure supplement 4).

The three samples obtained from the cultivation column (i.e. mix 1, 2 and 3) as well as the two reference samples (i.e. 0 hr and unlabeled) were transferred each in a 250 mL Erlenmeyer flask (or RAMOS flasks) containing 25 mL medium, adjusted to a cell density of 2 × 10⁷ cell mL⁻¹, and incubated at 300 rpm and 30°C.

The cell count was measured every 20 min between 1 and 3 hr after inoculation using a BD Accuri C6 flow cytometer (Becton, Dickinson and Company, Franklin Lakes, NJ). The samples were diluted with PBS at pH seven to <10⁶ cell mL⁻¹ and 20 µL sample were counted at 'medium' flow. The FSC-H thresholds was set to 80’000 in order to cut off most of the electronic noise. To correct the measured dry weight for the mass of iron beads in the sample, the iron beads were gated separately and counted as well. The data were analyzed using the Accuri CFlow Plus software.

As the cell volume and thus the cell specific dry weight (i.e. the weight of one cell) of mother cells changes with age, towards converting the measured cell counts to dry weight (biomass), we first determined the cell specific dry weight of mother/dead, m^mo/de, and daughter cells, m^da. After 3 hr, at the end of each batch cultivation, 20 mL of culture were filtered through a pre-weighed nitrocellulose filter with a pore size of 0.2 µm. The filter was washed once with distilled water, dried at 80°C for two days and afterwards weighed again. The total weight of iron beads attached to mother cells (here we assumed that one mother cell is attached to one iron bead; Janssens et al., 2015) and free beads, which was determined from the counted number of iron beads in the sample and the weight of one individual bead, was subtracted from the total dry weight of each sample. The bead weight had been determined to be 8.49 × 10⁻¹³ g per bead by filtration and weighting of a known amount of beads. Next, the cell specific dry weight of mother/dead and daughter cells was inferred from the measured population-average dry weight in the samples, m^meas, by following an in principle similar approach as done for the intracellular metabolite concentrations. Specifically, we assumed that dead cells (i.e. died mother cells) and mother cells have the same dry mass and that the dry mass of newly formed daughter cells does not change over the aging time points. Taking the fractional abundances of each cell population into account, the measured cell specific dry mass in each sample is given as,

 where m^meas_j,k is the measured population-average dry mass (unit g) after 3 hr cultivation in the sample j at the aging time point k, n^cell_j,k the total amount of cells in the respective sample, α_j,k the cell count specific fraction of mother cells, γ_j,k the cell count specific fraction of dead cells, m^mo/de_k the unknown cell specific dry mass (unit g) of mother or dead cells, β_j,k the cell count specific fraction of daughter cells and m^da the unknown cell specific dry mass (unit g) of daughter cells. Next, we formulated a least square regression problem of the form,

 where the matrix A contains all fractional abundances α_j,k + γ_j,k and β_j,k in every sample j at every aging time point k, the vector m the unknown cell specific dry weights m^mo/de_k and m^da at every aging time point k and the vector n all measured cell dry weights, m^meas_j,k, normalized by the total amount of cells in the sample, n^cell_j,k, in every sample j at every aging time point k. The regression problem in Equation 4 was implemented in R (Release 3.2.0) and the unknown cell specific dry weights, m, of mother/dead and daughter cells were identified using the function ‘lm’.

The inferred cell specific dry weights of mother/dead and daughter cells were then used to convert the measured cell counts to dry weight. At the beginning of each cultivation (t = 0) the total dry weight, Xt_t=0, is constituted of mother/dead and daughter cells, taking their fractional abundance into account, while in the following all new emerging cells are daughter cells. The total dry weight at every time t, X_t, is then given as,

 where X_t,j,k is the dry weight of the mixed population sample j of the aging time point k at time t, α_t=0,j,k + γ_t=0,_j,k and β_t=0,_j,k the cell count specific fractional abundances of mother/dead and daughter cells at the beginning of the cultivation, n_t=0,j,k the cell count at the beginning of the cultivation and n_t,j,k the cell count at the time t. Note that k refers to the cell age (i.e. aging time point) and t refers to the cultivation time at each aging time point (between 0 and 3 hr).

Additionally, the inferred cell specific dry weights of mother/dead and daughter cells were used to convert the cell count specific fractional abundances, α_j,k, β_j,k, and γ_j,k, in the dry mass specific fractional abundances of mother, daughter and dead cells, α^dw_j,k, β^dw_j,k, and γ^dw_j,k, in every sample j at every aging time point k:

0.3 mL samples were taken every 20 min from 1 to 3 hr after inoculation. To separate the cells from the medium, the samples were centrifuged at maximum speed for 3 min, the supernatant transferred onto a filter column (SpinX, pore size 0.22 µm), again centrifuged at maximum spend and the flow through was further analyzed. The glucose, pyruvate, glycerol, acetate and ethanol concentration was detected using an Agilent 1290 LC HPLC system equipped with a Hi-Plex H column and 5 mM H₂SO₄ as eluent at a constant flow rate of 0.6 mL min⁻¹. The injection volume was 10 µL and the column temperature was kept constant at 60°C. Glucose, glycerol, ethanol and acetate were detected by refractive index and pyruvate by UV (constant wave length of 210 nm) and the respective concentrations were determined using an external standard with known concentrations. The data were analyzed using the Agilent Open Lab CDS software.

The oxygen transfer rate (OTR) and carbon dioxide transfer rate (CTR) were determined from exhaust gas analysis using a respiration activity monitoring system (RAMOS) (Hansen et al., 2012). The RAMOS measurement flask, containing 25 mL medium, was inoculated with 2 × 10⁷ cell mL⁻¹ and the cultivation conditions were identical to the batch cultures used to determine the other physiological parameters. One RAMOS measurement cycle encompassed a 10 min measuring phase and a 20 min rinsing phase. The total consumption of oxygen and the production of carbon dioxide in a time interval were calculated from the mean of two consecutive OTR and CTR measurement cycles multiplied by the time.

To infer the physiological parameter of mother (mo), daughter (da) and dead (de) cells from the mixed population measurements, we formulated an ordinary differential equation model describing the dynamic change of biomass and extracellular metabolites during the 3 hr cultivation in each sample. To this end, we assumed that the physiology of daughter cells stays constant over all aging time points and that within the 3 hr cultivation the physiology of the mother cells stays constant. Finally, due to the short experiment time the evaporation of water and metabolites was neglected.

The total biomass in the sample is constituted of mother, dead and daughter cells and thus the differential mass balance can be formulated as,

Due to the short experiment time (3 hr) compared to their life span (>50 hr), we assumed that the amount of initial mother and dead cells stays constant (i.e. no new mother cells emerge and no mother cells die during the experiment). Thus,

and

 where X_j,k is the total biomass and X^mo_j,k and X^de_j,k the biomass of mother and dead cells in sample j at the aging time point k.

From Equation 9, 10 and 11, and follows that the change in total biomass is only due to the change in daughter cell biomass, X^da_j,k, which in turn can be either due to the emergence of new daughter cells originating from mother cells (i.e. budding of mother cells) or originating from daughter cells (i.e. budding of daughter cells). Thus, the change of the total biomass is given as,

 where µ^mo_k is the growth rate (unit h⁻¹) of mother cells and µ^da is the growth rate (unit h⁻¹) of daughter cells.

Reformulating the partial derivatives in Equations 10 and 11 and adding Equation 12 yields the change in dry mass specific fractional abundance of mother and dead cells as,

and

and plugging Equations 13 and 14 and in the differential biomass balance (Equation 9) yields the change in fractional abundance of daughter cells due to budding of mother and daughter cells as,

Next, the change in glucose concentration in the medium can be due to the uptake by mother and daughter cells as in,

 where c_glc,j,k is the measured glucose concentration (unit g L⁻¹) in sample j at the aging time point k, q_S^mo_k and q_S^da the specific uptake rates of mother and daughter cells and Y_XS^mo_k and Y_XS^da the biomass yields (unit g g_GLU⁻¹) of mother and daughter cells.

In a similar way, the mass balance for oxygen, carbon dioxide and other fermentation products can be formulated:

 where q_O2^mo_k, q_O2^da, q_P^mo_k and q_P^da are the biomass specific oxygen uptake and product (including carbon dioxide) excretion rates (unit g g_DW⁻¹ h⁻¹) of mother and daughter cells at the aging time point k and Y_O2S^mo_k, Y_O2S^da, Y_PS^mo_k and Y_PS^da the respective oxygen and product yields (unit g g_GLU⁻¹) of mother and daughter cells.

To increase robustness in the estimation, we stated that the mother and daughter cell physiology needs to fulfill the carbon balance within a certain range.

 where q^C_S and q^C_P are the specific carbon uptake and excretion rates (unit C-mol g_DW⁻¹ h⁻¹) of mother and daughter cells.

All three datasets were combined into one parameter estimation problem subject to the Equations 12–19. All parameters (including initial conditions) as well as the associated uncertainties were estimated using Maximum Likelihood estimation implemented in the software gPROMS ModelBuilder (Release 4.0, PSE software systems) with the MINLP solver SRQPD where a constant variance (error model) was assumed for all measurements.

To determine the intracellular fluxes at different cell ages from the inferred metabolite concentrations and physiologies, we made use of a recently published computational inference method (Niebel et al., 2019). This method rests on a combined thermodynamic and stoichiometric network model of cellular operation, M(v,lnc)≤0 (Equation 20), consisting of a mass balanced metabolic reaction network, in which the reaction directionalities are constraint by the associated changes in Gibbs energy – as a function of the metabolite concentrations c – through the 2^nd law of thermodynamics. Additionally, the Gibbs energy, which is dissipated through metabolic operation (i.e. the sum of all metabolic processes, MET) is balanced with the Gibbs energy exchanged with the environment through exchange processes (i.e. the production and consumption of metabolites, EXG),

 where S_ij is the stoichiometric coefficient of the i^th reactant (i.e. metabolite) in reaction j, v_j the rate of the reaction j (i.e. the flux through this reaction), Δ_rG’(ln c_j) the Gibbs free energy of reaction of the metabolic process j and Δ_fG’(ln c_i) the Gibbs free energy of formation of the reactant i.

The published, and here used, model for Saccharomyces cerevisiae encompasses the metabolic processes of glycolysis, gluconeogenesis, tricarboxylic acid cycle, amino acid-, nucleotide-, sterol-synthesis and considers the processes’ location in the cytosol, mitochondria and extracellular space. To account for cofactor turnover due to the combatting of reactive oxygen species, which is known to occur at high replicative ages (Ayer et al., 2014), the model was extended by reactions describing the oxidation of NADH and NADPH through glutathione in the cytoplasm as well as the glutathione exchange (i.e. a sink and a source). This exchange does not represent any direct metabolic process but needed to be included since the glutathione metabolism is not part of this model.

             nadh[c] + gthox[c] => nad[c] + (2) gthrd[c]
             nadh[c] + gthox[c] => nad[c] + (2) gthrd[c]
                               gthox[c] <=>
                               gthrd[c] <=>

A more detailed description of this model and its implementation can be found in Niebel et al. (2019).

Using this model and the inferred age-dependent metabolite concentrations and physiologies, we formulated a regression problem minimizing the weighted residual sum of squares, RSS(y) (Equation 21). As data we used (i) the inferred yields, ${\displaystyle {\tilde{Y}}_i^{(k)}}$ (i∈PY… physiological yield), (ii) the inferred metabolite concentrations ${\displaystyle {\tilde{c}}_i^{(k)}}$ (i∈MC1∪i∈MC2… metabolite concentration set 1 or 2 (see below)), both of daughter and aged mother cells at a replicate age of 0, 10, 20, 44 and 68 hr and (iii) standard Gibbs energies of reaction, ${\displaystyle {\mathrm{\Delta }}_r{\tilde{G}}_j^{\prime o}}$. The later were determined (including uncertainty) using the component contribution method (Noor et al., 2013) and as this was not possible for all standard Gibbs energies, to prevent overfitting, the regression was regularized by the Lasso method (Hastie et al., 2011).

To ensure the same thermodynamic reference state (i.e. the same standard Gibbs energies of reactions) in all experimental conditions, we bundled all datasets in on regression problem and indexed the model (Equation 20) over the experimental conditions k.

where #n_Y and #n_c are the number of inferred yields and metabolite concentrations, #n_CCM the number of standard Gibbs energies of reaction, which could be estimated by the component contribution method and #n_unk the number of reactions where no standard Gibbs energy of reaction could be calculated. The residuals were weighted by the respective prediction uncertainty, indicated by the superscript SE. Metabolites can occur in the cytoplasm and/or in the mitochondrial space (MC1… metabolites occurring in one compartment and MC2… metabolites occurring in two compartments). Thus, we stated that the sum of the metabolite concentrations in the respective compartments, weighted by the fractional compartmental volume (0.9 for the cytoplasm and 0.1 for the mitochondrial space), had to be equal to the inferred (cell-averaging) concentration. Last, to facilitate the convergence of the optimization and for an easy conversion of reaction rates to yields, the glucose uptake rate, v_{glc-D_EX}, was constraint to a value of 1 mmol gDW^-1 h^-1.

The regression analysis was implemented in the mathematical programming system GAMS (GAMS Development Corporation; General Algebraic Modeling System (GAMS) Release 24.2.2. Washington, DC, USA).

To obtain a picture of the intracellular flux distribution, we formulated the solution space, Ω^reg (Equation 22), of the optimal regression solution, indicated by an *,

Within this solution space we then minimized the ‘absolute sum of fluxes’,

The optimization problem in Equation 23 was implemented in the mathematical programming system GAMS (GAMS Development Corporation; General Algebraic Modeling System (GAMS) Release 24.2.2. Washington, DC, USA).

For microscopy experiments, cells from exponentially growing batch cultures were used to load a microfluidic device (Huberts et al., 2013; Lee et al., 2012). Individual cells were monitored using an inverted fluorescence microscope (Eclipse Ti-E; Nikon) housed in an custom-made microscope incubator (Life Imaging Services GmbH) that retained at a constant temperature of 30°C. During the experiment, cells were continuously fed with fresh medium. An LED-based excitation system (pE2; CoolLED) was used for illumination, and images were recorded using an Andor 897 Ultra EX2 EM-CCD camera. NAD(P)H autofluorescence (excitation at 365 nm using a 357/44 nm filter and a 409 nm beam-splitter, 200 ms exposure time, 15 % light intensity, 435/40 nm emission, EM gain 1) was recorded every 60 min to minimize phototoxic effects, and brightfield images every 10 min to reliably track individual cells and determine their division times. A CSI S Fluor 40x Oil (NA = 1.3; Nikon) objective was used for NAD(P)H. Automated hardware (PFS, Nikon) was used for correction of axial focus fluctuations during imaging.

Cell segmentation for estimation of cell volume and fluorescence intensity took place in a semi-automated manner using the ImageJ plugin BudJ (Ferrezuelo et al., 2012). For cell volume estimation, brightfield images captured with the 60x objective were used. Fluorescent intensity measurements were corrected for background fluorescence using the Rolling Ball Radius algorithm of ImageJ. For budding rate estimations on the basis of single-cells, the doubling time, td, (time from bud emergence to bud emergence) was measured for each cell in 60x brightfield images, and the budding rate for each doubling event (ln(2) t_d^-1) was calculated.

Cells from an exponentially growing culture (minimal medium; Verduyn et al., 1992) supplemented with 1 % (w/v) glucose were loaded in two identical microfluidic devices located on one cover glass. Minimal media supplemented with 0.5 % (w/v) glucose with and without 0.1 % (v/v) ethanol were constantly supplied into the two microfluidic devices, respectively. The cells in the microfluidic devices were monitored simultaneously by taking bright-field images every 10 minutes for more than 5 days (halogen lamp with a UV-blocking filter, 60x objective). The time points of budding, death and washout loss were recorded for individual cells using a custom macro in ImageJ. The number of budding events and fate (death or washed) of the individual cells in both microfluidic devices were used to assess the replicative age-associated survival via the Kaplan-Meier estimator. The analysis was implemented using the Lifelines (0.9.4) module in Python (2.7.13). The mean survival and its standard error were calculated using the Survival (2.43-3) package in R (3.4.1) integrating the survival curves until 44 buds (the maximal number of buds per cell in two conditions).

• Matthias Heinemann

• Vakil Takhaveev
• Matthias Heinemann