Stochastic and Age-Dependent Proteostasis Decline Underlies Heterogeneity in Heat-Shock Response Dynamics

Caenorhabditis elegans , the aging process shows a high degree of variability among individuals sharing the same genetic and environmental background. In this study, we investigated population heterogeneity of the aging process using the model organism C. elegans . The nematode C. elegans is a small ( ≈ 1 mm in length) and fast Significant non-genetic stochastic factors affect aging, causing lifespan differences among individuals, even those sharing the same genetic and environmental background. In Caenorhabditis elegans , differences in heat-shock response (HSR) are predictive of lifespan. However, factors contributing to the heterogeneity of HSR are still not fully elucidated. Here, the authors characterized HSR dynamics in isogenic C. elegans expressing GFP reporter for hsp-16.2 for identifying the key contributors of HSR heterogeneity. Specifically, microfluidic devices that enable cross-sectional and longitudinal measurements of HSR dynamics in C. elegans at different scales are developed: in populations, within individuals, and in embryos. The authors adapted a mathematical model of HSR to single C. elegans and identified model parameters associated with proteostasis—maintenance of protein homeostasis—more specifically, protein turnover, as the major drivers of heterogeneity in HSR dynamics. It is verified that individuals with enhanced proteostasis fidelity in early adulthood live longer. The model-based comparative analysis of protein turnover in day-1 and day-2 adult C. elegans revealed a stochastic-onset of age-related proteostasis decline that increases the heterogeneity of HSR capacity. Finally, the analysis of C. elegans embryos showed higher HSR and proteostasis capacity than young adults and established transgenerational contribution to HSR heterogeneity that depends on maternal age.


Introduction
Medical breakthroughs, such as antibiotics and vaccines, as well as improved health care have led to a doubling in human life expectancy in developed countries over the last 200 years. [1]he proportion of European population older than 65 years has reproducing hermaphrodite that is particularly amenable to large scale studies because of the ease in maintaining a genetically identical (isogenic) population under tightly controlled conditions.Moreover, a typical lifespan of between three to four weeks makes the organism highly suitable for aging research.Heterogeneity in the aging process related to non-genetic and other stochastic factors in a C. elegans population can be studied by keeping the genetic background and environment uniform.10][11][12][13] Interventions that alter the mean lifespan of C. elegans have also been shown to affect the heterogeneity in the lifespan by a linear temporal scaling [14] -an increase in the mean lifespan will translate to the same fold increase in the standard deviation of lifespan.
Heterogeneity in the heat-shock response (HSR)-an evolutionarily conserved molecular mechanism to maintain protein homeostasis (proteostasis)-has been tied to the variability in lifespan of the model organism C. elegans.A stronger HSR, specifically the expression of hsp-16.2[17] This suggests that the capability of nematodes to respond to heat stress predicts longevity.[20] But, whether the variability in these factors is set at the start of life and stays constant or whether it increases with age remains unclear.
In our study, we used TJ375 [15] C. elegans, carrying GFP reporter for hsp-16.2expression, to study the non-genetic and stochastic factor(s) contributing to heterogeneity of HSR in isogenic C. elegans populations.To characterize HSR capacity among individuals in large C. elegans cohorts, we developed two novel microfluidic systems; one for screening fluorescent protein reporter expression in large numbers of individuals at specific time points, and another for sorting the population according to the characteristics of such expression.We also employed a previously developed microfluidic device for highresolution time-lapse GFP measurements of HSR dynamics in single C. elegans. [21]At the population level, we found that the heterogeneity in HSR activity in isogenic TJ375 (Phsp-16.2::GFP)populations generally follows a Gaussian distribution.Through sorting and screening of TJ375 after heat-shock, we observed that isogenic C. elegans exposed to the same environment exhibit diverse HSR dynamics.We confirmed the heterogeneity in HSR dynamics by time-lapse fluorescence imaging of individual TJ375 after heat-shock.By applying mathematical modeling and sensitivity analysis to time-series HSR data of single worms, we identified the variability of protein turnover parameters (translation and degradation rate constants) as a major driver of population heterogeneity of HSR dynamics.Guided by model simulations, we sorted TJ375 populations based on their HSR capacity, and confirmed by lifespan assay that higher HSR capacity is associated with a longer lifespan.In addition, by characterizing the heterogeneity of HSR in TJ375 populations of two different ages (day-1 and day-2 adulthood), we observed a stochastic age-associated decline in protein translation and degradation rate constants that explain a decrease in HSR capacity and an increase in population heterogeneity of HSR dynamics.Compartmental analysis further revealed differences in HSR dynamics along the body of TJ375 C. elegans.While agedependent decrease of HSR capacity coupled with an increase of heterogeneity is also seen in the compartments, the differences of HSR dynamics across the compartments lower with age.Finally, our characterization of HSR dynamics in embryos of TJ375 from mothers of different ages (day-1, day-2, and day-4) indicated a maternal age-dependent effect of transgenerational contribution to the heterogeneity of HSR dynamics, where day-2 mothers produced embryos with the highest protein turnover and HSR capacity.

Cross-Sectional Population Measurements show Heterogeneity of HSR in Isogenic C. elegans
We used the transgenic nematode TJ375 (Phsp-16.2::GFP)carrying a plasmid with GFP expression under the control of heat-shock protein HSP-16.2promoter, [15] for in vivo assessments of HSR in C. elegans (Figure 1a,b).Consistent with previous reports, [22] we found an agreement between the time-course expression of the endogenous HSP-16.2 and that of the GFP via Western Blotting (Figure S1, Supporting Information).Here, we developed a novel microfluidic system for measuring the expression of fluorescent protein reporters in C. elegans in a high-throughput and biocompatible manner (Figure 1c, Figure S2, Supporting Information), as detailed in the Experimental Section.We employed this microfluidic system to study the heterogeneity of HSR in an isogenic nematode population by GFP expression quantification of transgenic nematodes at multiple time points after heat stress.To the best of our knowledge, a time-course (dynamic) characterization of HSR for a constant nematode population upon heat-shock has not been previously performed.
The details of the microfluidic system and the experimental procedure are described in the Experimental Section.Briefly, adult TJ375 C. elegans suspended in M9 buffer enters the microfluidic device through main inlet and flows in a serial and non-overlapping manner through a main straight channel, supported by secondary buffer inlets (Figure 1c).The optical detection system integrates a detection volume, orthogonal to the direction of the flow stream.In this manner, when an animal passes through the detection volume, GFP expression (fluorescence intensity) along the body of individual C. elegans is recorded (Figure S2, Supporting Information).To assess the heterogeneity of HSR, we heat-shocked isogenic day-2 adult TJ375 C. elegans for 1-h at 37 o C. We then passed the nematodes through the microfluidic system and measured the whole-body GFP expression of the population at 2, 4, 8, 12, 24, 32, and 48 h after the heat-shock (Figure 1d).We recorded the whole-body GFP fluorescence intensities of 2280 individual transgenic C. elegans.Noteworthy, non-heat-shocked C. elegans did not express any detectable GFP in our system.For all time points, the HSR heterogeneity in the nematode population follows a Gaussian distribution (Figure 1e and Figure S3, Supporting Information).The standard deviation of fluorescence intensities increases initially and then remains approximately constant beyond 4 h after heat-shock (Figure S3, Supporting Information), which is in agreement with a previous study using the same transgenic nematodes. [15]The HSR heterogeneity does not differ markedly between TJ375 and a mutant with epigenetically suppressed HSR, LSD2088 (Phsp-16.2::GFP,jmjd-3.1(gk384)) [23]b) Schematics of the GFP expression on the TJ375 C. elegans strain following heat-shock.GFP is only apparent on C. elegans hours after heat-shock.c) Schematic of the microfluidic device for C. elegans fluorescent protein screening.The "C. elegans inlet" i) incorporates wallembedded PDMS blades to untangle animals as they enter the device.Additional buffer inlets ii) enable efficient spacing of C. elegans in the main channel.iii) Fluorescence detection volume where the optical detection system is focused onto.Animals exit via the device outlet (iv) and calibration is performed using the reference channel (v).The above measurements of the population give crosssectional information on the heterogeneity of time-dependent GFP expression on a TJ375 C. elegans cohort.The overt heterogeneity may arise from differences in the HSR dynamics among individuals (i.e., each animal following a different HSR time-course) or from static variability (i.e., each animal following the same HSR time-course, but with some variation in the intensity of GFP expression around the time-course trajectory).

Sorting and Tracking of Isogenic C. elegans Subpopulations Indicate Heterogeneity of HSR Dynamics
To shed light on the source of the observed HSR heterogeneity among TJ375 C. elegans, we adapted the aforementioned microfluidic system to allow for the sorting of nematode populations into up to three groups according to their whole-body fluorescence signal in a high-throughput and biocompatible manner (Figure 2a, Figure S2, Supporting Information), as detailed in the Experimental Section.Using this microfluidic system, we sorted a day-2 adult TJ375 C. elegans population after a 1-h 37 o C heat-shock (see Figure 2b).First, the heat-shocked population was sorted into two groups based on the whole-body GFP fluorescence at 4 h after heat-shock: a High sub-population for nematodes with a whole-body fluorescence above the median level, and a Low sub-population for those with a whole-body fluorescence below the median level.Subsequently, at 8 h after heat-shock, we repeated the sorting process for each sub-populations.For the High sub-population, we separated the top 16% from the bottom 84% of fluorescence intensities, and respectively labeled them as High-High and High-Low .elegans inlet" i) incorporates wall-embedded PDMS blades to untangle animals as they enter the device.Additional buffer inlets ii) enable efficient spacing of animals in the main channel.iii) Fluorescence detection volume where the optical detection system is focused onto.Sorted animals exit via one of the device outlets (iv) and calibration is performed using the reference channel (v).Lateral on-chip valves are labeled (vi).b) Experimental layout of C. elegans population sorting experiments.c) Double sorting and subsequent GFP intensity measurement of a C. elegans population (n = 2600), following layout shown in (b), at different time points after heat-shock.Each measured population is shown with violin plots.Initial sorting is performed 4 h after heatshock, dividing the population into a top and bottom 50%.A secondary sorting step (8 h after heat-shock) further separated the top 50% into its top 16% and bottom 84%.Similarly, the initial bottom 50% of the population is further sorted into a bottom 16% and a top 84%.The four resultant groups are screened at posterior time points (12, 25, and 32 h).The change in the mean value of each group at each sorting or screening time is shown with a distinct color line per group, that is, green: "high-high" group, blue: "high-low", yellow: "low-high", and red: "low-low".
(see Figure 2c).Mirroring this sorting strategy, we sorted the Low sub-population at 8 h after heat-shock to separate the bottom 16% from the top 84%, and respectively labeled them as Low-Low and Low-High (see Figure 2c).Separation in this manner, captures the upper (lower) tail of the High (Low) subpopulation at one standard deviation above (below) the mean.We further screened GFP expression of the four resulting sub-populations at several time points after the sorting experiments-at 12, 25, and 33 h after heat-shock (colored in red, yellow, blue, and green; Figure 2c).By performing these doublesorting experiments, we were able to follow dynamically distinct sub-populations and better understand the heterogeneity of HSR dynamics in the TJ375 C. elegans.
We observed that within 4 h of the first sorting (i.e., 8 h after heat-shock), the whole-body GFP average fluorescence of each sub-population returned to an approximately Gaussian distribution (see Figure 2c).Comparing the High and Low subpopulations at 8 h after heat-shock, we also noted a large overlap between their whole-body fluorescence average distributions.That said, the High sub-population maintained a statistically significant higher average whole-body fluorescence signal than the Low sub-population (two sample, two-tailed, t-test p-value < 1 × 10 −5 ).Among the four sub-populations resulting from the second sorting experiment, the High-High sub-populationthe strongest responders to heat-shock-maintained the highest average fluorescence over the remaining duration of the experiment (green line; Figure 2c), having a maximum recorded fluorescence at 8 h after heat-shock.The Low-Low subpopulation-the weakest responders-had the lowest average fluorescence (red line; Figure 2c), with a maximum recorded GFP expression at 12 h after heat-shock.The High-Low and Low-High sub-populations displayed strikingly similar HSR dynamics post-sorting.But interestingly, the Low-High group had a slightly (but statistically significant) higher average fluorescence than the High-Low group at and beyond 12 h after heatshock (two sample, two-tailed, t-test p-value = 3.99 × 10 −4 ; yellow and blue lines in Figure 2c, Table S1, Supporting Information).The cross-over in the average fluorescence between the High-Low and Low-High, as well as between High-High and High-Low, suggest that a number of animals have a slow response to heatshock.Finally, differences among the four sub-populations were observed to diminish over time.The results from these double-sorting experiments support the notion that the heterogeneity of HSR in TJ375 C. elegans is a consequence of variability in the HSR dynamics among the animals, that is, each nematode follows a different time-course HSR trajectory.

Time-Lapse Imaging of Individual C. elegans Reveals the Characteristics of HSR Dynamics and their Variability in an Isogenic Population
Our initial sorting-screening experiments suggest a significant animal-to-animal heterogeneity in HSR dynamics, even within an isogenic population that is kept in the same environment.To better understand the key molecular mechanisms that give rise to this heterogeneity, we obtained a longitudinal observation of the HSR time-course by measuring GFP fluorescence in single TJ375 C. elegans using time-lapse microscopy.For this purpose, we adapted a microfluidic platform for long-term and high-resolution imaging of immobilized C. elegans that was previously developed in our laboratory [21] (Figure 3a,b; details on the device in the Experimental Section section).We confirmed that the device itself does not induce HSR response in TJ375 C. elegans for the same duration of the experiment (Figure S5a, Supporting Information).We characterized the HSR dynamics after 1-h 37 o C heat-shock for fifteen day-2 adult TJ375 C. elegans (Figure 3c, Figure S5b, Video S1, Supporting Information), measuring the GFP expression intensity hourly over 48 h.The results of this experiment are presented in Figure 3d, where the different colored lines illustrate the whole-body fluorescence of individual C. elegans as well as the average GFP emission across the 15 C. elegans (bold black line).In general, GFP expression in individual worms shows a sharp increase after heatshock, peaking between 6 and 24 h after heat-shock and then decreasing with a typical exponential decay dynamic.However, it is noticeable that the rate of the GFP expression increase following heat-shock differs between individuals.In summary, the longitudinal time-course measurements for single animals confirm a significant animal-to-animal heterogeneity in the HSR dynamics, as suggested by our screening and sorting experiments.

Mathematical Modeling and Sensitivity Analysis Identify Translation and Degradation Rate Constants as a Major Driver of Individual Variability of HSR Dynamics
In the following, we employed a mathematical model of HSR to shed light on the possible underlying biological mechanism(s) driving the observed heterogeneity among isogenic C. elegans.[26][27][28][29][30][31][32][33] Here, we adapted an Ordinary Differential Equation (ODE) model of HSR originally developed for human HeLa cells by Scheff et al. [34] The model was among the most comprehensive models of HSR, built on previous modeling works by Petre et al. [35] and Szymanska and Zylicz. [36]riefly, the HSR model describes the network of molecular reactions that responds to an increase in the level of protein denaturation caused by heat insults.The original model parameters were obtained by fitting to literature compiled time-course HSP data of HeLa cells.For the specific purpose of our study and also for a direct comparison to our data of time-course GFP fluorescence in TJ375 (Phsp-16.2::GFP)C. elegans, we added a balance equation for the GFP reporter to the original model (details in Experimental Section, Table 1, Table 2, and Figure 4a).Here, we assumed that the rate constants of the GFP transcription and translation mirror those of the HSP, based on the fact that they are controlled by the same promoter.We also set the degradation rate constant of the GFP to equal that of the HSP, since both proteins followed similar dynamics after their peaks in our western blotting experiments (Figure S1, Supporting Information).
We fit the model parameters to the GFP time-course data of individual C. elegans in two stages.In the first stage, we fit the GFP concentration predicted by the model to the average C. elegans whole-body fluorescence measurements of the 15 day-2 adult TJ375 C. elegans (Figure 3d, Table 3 and Experimental Section).We performed the parameter estimation by using a hybrid scatter search method-a stochastic optimization method-100 times and generated 100 model parameter combinations with various quality of data fit, the best of which was used for the second stage of single-animal modeling.Figure 4b illustrates the model simulation of the key components of HSR in C. elegans after a 1-h 37 °C heat-shock (y-axis not to scale).Misfolded protein (MFP) increases sharply during the 1-h heat-shock (not shown).HSP captures MFP for refolding and MFP levels drop immediately after heatshock.The binding of HSP to MFP frees HSF-1 that would have otherwise formed a complex with HSP.The free HSF-1 undergoes stepwise polymerizations to produce HSF-1 trimers, and after translocating to the nucleus, the trimer binds to the heat-shock element (HSE) to upregulate the expression of HSP (and GFP).The level of HSP mRNA rises after heat-shock until MFP returns to its pre-heat-shock level.The amount of HSP increases along with its mRNA, though with a slight delay.The GFP signal follows a time-course similar to that of the HSP, consistent with the Western blot data (Figure S1, Supporting Information).
In the second stage, we performed parameter fitting to the GFP time-courses of individual C. elegans to obtain 15 different parameter sets, one for each C. elegans in the experiment.Here, we assumed that the parameter values of individual animals differ amongst each other and from the average C. elegans only for a few parameters.More specifically, we set the majority of the single-worm model parameters to the values obtained from the fitting of the average whole-body GFP fluorescence above, and re-fit only a select set of parameters to the single-animal whole-body GFP data.To this end, we first performed the parametric sensitivity analysis to identify the key model parameters that affect the time to reach maximum GFP and the maximum GFP value.We ran the sensitivity analysis for each of the top 10 out of the 100 parameter combinations obtained from the average worm parameter fitting, and ranked the parameters based on the average sensitivity coefficients (details in the Experimental Section).By using the top 10 parameter combinations, rather than just the best parameter, the result of the sensitivity analysis would not be highly dependent on a specific model parameterization.The sensitivity analysis (Figure S6, Supporting Information and Experimental Section) points to the basal level of HSF (HSF0) and parameters related to the protein turnover as the key model parameters that govern the HSR dynamics-the speed and strength of HSRimmediately after heat-shock.The result of our sensitivity analysis agrees with a recent study showing variability in basal protein expression, notably HSF, determines the overt heterogeneity of HSR dynamics in HeLa cell population. [18]39] Based on the results of sensitivity analysis and other considerations related to the nature of various parameters (see the Experimental Section for more detail), we performed single-worm model parameter fitting by optimizing parameters associated with protein translation and degradation rate constants.Note that changes in protein translation and degradation affect protein expression capacity of the organism [20] which we accounted for by scaling the protein the total protein and HSF amount by the ratio of the translation and degradation changes (see Experimental Section).By only varying these two parameters, we were able to obtain accurate model fits to the GFP time-courses of individual C. elegans (Figure 4c).The distribution of the protein turnover parameters for the individual C. elegans is shown in Figure 4d, where the numbers next to the parameter points correspond to the identity of the single C. elegans in Figure 3d.Thus, the results suggest that differences in the protein turnover capacity are likely a major driver of the heterogeneity in the HSR dynamics among individual C. elegans in isogenic populations.

The Stochastic Collapse of Proteostasis during Early Adulthood Increases the Heterogeneity of HSR in Isogenic C. elegans Population
Protein translation, trafficking, and degradation are essential for the functionality of cellular processes.[46] In the above experiments, we used day-2 adult C. elegans that may have potentially already experienced such proteostasis decline.We thus hypothesized that younger C. elegans should demonstrate a better (faster and stronger) response to the same heat-shock stress.To confirm this hypothesis, we repeated our assessments of HSR dynamics for a cohort of day-1 adults exposed to the same 1-h 37 °C heatshock treatment.The HSR assessment using our high-throughput microfluidic system indicated that the average whole-body GFP expression of day-1 adult TJ375 population reaches its maximum earlier than day-2 adults (8 h versus 12 h after heat-shock; Figure 5a,b).In order to further confirm that the source of the differences between these two cohorts is associated with proteostasis decline, we characterized the single-animal HSR dynamics of day-1 adult C. elegans (n = 10) using timelapse fluorescence microscopy (Figure 5c, Figure S7a, Supporting Information).We found that the GFP time courses of day-1 adults have, on average, a higher maximum fluorescence value than day-2 adult nematodes (195.11± 63.78 a.u.for day-1 versus 147.74 ± 35.16 a.u.for day-2 adults, two-sample t-test p-value = 0.025), and also that they peak earlier (8.7 ± 1.5 h for day-1 versus 14 ± 3.48 h for day-2, two-sample t-test p-value = 1.56 × 10 −4 ).We then applied our mathematical model to the combined GFP expression time course curves from day-1 and day-2 TJ375 adults (n = 25) and repeated the two-stage parameter fitting process (Figure S7b, Supporting Information).Again, we were able to explain the variability of HSR dynamics among individual nematodes by varying protein translation and degradation rate constants among individual nematodes.
The distribution of the protein translation and degradation rate constants for day-1 and day-2 adult C. elegans are depicted in Figure 5d (Table S2, Supporting Information).The results of the model fitting show that younger C. elegans (day-1 adults) have on average higher rate constants for protein translation and degradation, both of which are indicators of a higher proteostasis capacity, confirming our hypothesis.Although a positive correlation is observed for day-2 adults between translation and degradation rates (ρ = 0.80) this correlation is not observed for adult day-1 (ρ = −0.02);suggesting that the processes might not be coupled until later stages in the life of C. elegans.Besides the lower protein turnover parameters, the day-2 cohort also displays a higher heterogeneity (variance) in the rate constants (Table S2, Supporting Information).The increased heterogeneity agrees with a previous study of change in protein turnover in C. elegans aging. [47]Moreover, when comparing the GFP expression curves of day-1 and day-2 adult C. elegans (Figure 5c), we observed an overlap, where a fraction of day-2 adult C. elegans (especially, day-2 C. elegans number 7, 11, and 12) exhibited GFP expression curves that are similar to the typical GFP expression curves of day-1 adult C. elegans.A further comparison of the protein turnover parameters shows that these specific day-2 adult C. elegans indeed have translation and degradation rate constants that are higher than those of the other day-2 animals and closer to those of day-1 adults (Figure 5d and Table S2, Supporting Information).These comparisons suggest that aging-associated decline of  protein turnover does not occur at equal rates for nematodes in a given population.Indeed, other studies have reported a significant stochasticity in the onset and the rate of deterioration in protein turnover within isogenic C. elegans populations. [48]e further investigated the impact of protein turnover heterogeneity on HSR capacity.For this purpose, we focused on the early time-course of HSR as an indicator of its capacity.More specifically, we associated the HSR capacity in individual C. elegans with the speed and strength of the GFP expression increase in response to a heat-shock.We used the HSR model to simulate the time to reach the maximum GFP expressiona marker of speed-and the maximum level of the GFP expression-a marker of strength of the HSR.Model simulations show that the speed of the HSR in individual C. elegans is positively correlated with their protein translation rate constants (Figure S7c, Supporting Information), that is, a higher translation rate constant gives a faster GFP increase in response to heat-shock.Meanwhile, the maximum GFP expression level is proportional to the ratio between protein translation and degradation rate constants among single C. elegans (Figure S7d, Supporting Information).
We subsequently defined a single combined marker for HSR capacity in individual C. elegans which would account for both the speed and maximum level of GFP expression and computed the slope of GFP expression by taking the ratio between the maximum GFP level and the time to reach the maximum and used this value as an indicator of the HSR capacity that each C. elegans has. Figure 5e depicts the isolines (contours) of the GFP slopes as a function of the rate constant of protein translation (x-axis) and protein degradation (y-axis) (details in Experimental Section).The isolines show that a faster translation corresponds to a higher slope, whilst a higher degradation corresponds to a lower slope.Thus, the effect of declining protein translation on the HSR capacity can be abrogated by lowering protein degradation.Furthermore, the protein translation rate constant exerts a stronger influence on the GFP slope than the protein degradation rate constant.The dependence of the slope on protein translation strengthens with lower level of protein degradation, as indicated by denser isolines at lower degradation rate constants.By comparing the locations of the parameters of day-1 adult C. elegans with those of day-2 adult C. elegans with respect to the isolines in Figure 5e, we noted that day-1 adult C. elegans have generally higher GFP slope than day-2 adults (see Table S3, Supporting Information, two-tailed, two-sample t-test p-value < 0.0001).Thus, this suggests that the decline in the protein turnover during early adulthood is detrimental to the HSR capacity of the animal.
In order to confirm the influence of the protein turnover rate on the HSR capacity, we crossed the ncl-1(e1942) knockout mutants that have increased amounts of ribosomes and almost double the amount of general protein translation, [49] into TJ375 (Phsp-16.2::GFP),producing a transgenic strain LSD2089 (Phsp-16.2::GFP,ncl-1(e1942)).The transgenic strain is analogous to almost doubling the protein translation rate constant in the HSR model.We measured the GFP expression time-course of heat-shocked day-1 and day-2 adults LSD2089 C. elegans and compared them with that of TJ375 (Figure 5f).We found that while the times to reach the maximum GFP expression in the populations are similar between same age LSD2089 and TJ375, the mean GFP expression values of LSD2089 are 1.8 times compared to TJ375.That is, the transgenic strain LSD2089 has a GFP slope that is approximately double that of TJ375 strain, an increase that is consistent with the doubling of protein translation rate constant and with the trend indicated by the GFP slope isolines in Figure 5e-a higher protein translation rate constant corresponds to a higher slope.The deviation in the GFP time-courses between day-1 and day-2 adults of LSD2089 C. elegans mirrors that of TJ375.These data thus suggest that despite an improved HSR activity conferred by a higher protein translational capacity, LSD2089 C. elegans still experience a similar proteostasis collapse during their early adulthood.

Mathematical Model Simulations of HSR Identify Subpopulations of C. elegans with Early Predictive Markers of Lifespan
Previous studies showed that the strength of HSR in TJ375 C. elegans as measured by their HSP expression levels is predictive of lifespan. [15,17]Based on our modeling results, which point out that assayed C. elegans of the same chronological age appear to have different protein-turnover capacities, we hypothesize that the HSR capacity, and thus the protein-turnover machinery of each animal, is the underlying element of this lifespan predictor.For hypothesis validation, we performed model simulations to determine the sorting criterion that demarcates the top and bottom HSR capacity among day-2 adult C. elegans.Briefly, we generated an ensemble of 10 000 in-silico C. elegans models, where the models differed from each other only in the values of protein translation and degradation rate constants (Figure 6a,b).The model parameters were generated by a random sampling of the parameters, assuming a bivariate Gaussian distribution  with the mean and variance equaling to those of single C. elegans parameters in Figure 5c.We simulated each singleanimal model to produce a GFP expression time-course.The collection of 10 000 such time courses is depicted in Figure 6b.The simulated GFP curves agree well with our screening measurements from the high-throughput experiments, as shown in Figure 1e (Figure S8, Supporting Information), as well as with the sorting-screening experiment shown in Figure 2c (Figure S9, Supporting Information).
Based on the simulations of the model ensemble in Figure 6a, we found that taking the top and bottom 5% GFP expression of the population at 5 h post heat-shock will result in significantly different HSR capacity with the average GFP slope of the top 5% being fourfold higher than that of the bottom 5% (25.54 a.u.h −1 for "high" group versus 6.41 a.u.h −1 for "low" group).The sorting criteria above would capture C. elegans corresponding to the parameter space depicted in Figure 6c (the top 5% in green and the bottom 5% in blue) and to the GFP expression time-courses shown in Figure 6d.The model simulations shown in Figure 6d suggest that the two sub-populations should comprise C. elegans with starkly different GFP slopes and thus HSR capacities.We implemented the sorting criterion as determined by the above model simulations in a sorting experiment using a day-2 adult TJ375 C. elegans population following heat-shock.We used our highthroughput sorting microfluidic device (Figure 2a) to sort and isolate the two sub-populations.Figure 6d gives the combined survival curves for the top and bottom 5% GFP expressing day-2 adult C. elegans at t = 5 h after heat-shock for 3 biological repeats; 2 of those repeats were obtained using a fullyautomated lifespan machine [50] and one by manual inspection, as detailed in Experimental Section and Supporting Information S1.In combination, the top 5% sub-population has a 4.4-day (25%) higher mean lifespan than the bottom 5% in these independent experiments (p-value < 0.001 Figure 6e, Supporting Information S1b-d).The result of our lifespan study using model-based sorting criteria thus validates our hypothesis that the HSR capacity-as measured by the GFP slope-is an early biomarker for longevity among a genetically identical C. elegans population-as early as adult day-1.

Heterogeneity of HSR Dynamics Exists within Individual C. elegans
Individual tissues within an organism have been shown to respond not only to their own proteostasis disruption but also to that of adjacent tissues in a coordinated fashion. [51]Tissuespecific HSR regulation has further been tied to distinct cellular proteome and specific profiles of heat-shock inducible genes, which give rise to specific tissue-dependent contributions to the whole-organismal proteostasis machinery. [52]To examine tissue-specific HSR dynamics in C. elegans, we re-analyzed the time-lapse images of individual day-1 and day-2 adult nematodes upon heat-shock.For this purpose, we sectioned the images into five equal-sized compartments and monitored the temporal expression of GFP in distinct tissue regions, namely head, anterior intestine, central intestine, posterior intestine, and tail (Figure 7a).In each of these tissue sections, we obtained the mean GFP expression over 30-h time course for day-1 adult (Figure 7b) and day-2 adult C. elegans (Figure 7c). Figure 7d (Table S4, Supporting Information) compares the GFP expression of day-1 and day-2 adult datasets in different tissue sections.
Comparing across compartments, we observed that the central intestine compartment exhibits the quickest response to heat-shock, judging from the initial rise (slope) of the GFP.Comparing between day-1 and day-2 cohorts, except for the head section, the GFP expression of day-1 worms peaks sooner than that of day-2 worms.In the central intestine section, day-1 nematodes on average take 8.4 h while day-2 animals on average take 12.6 h to reach the GFP peak (Figure 7b-d).In the anterior intestine, the average time to reach the maximum GFP is 10 h for day-1 versus 13.5 h for day-2 (Figure 7b-d).In the posterior intestine compartment, the average peak times are 9 h for day-1 and 16.7 h for day-2 (Figure 7b-d).In the tail section, the average time to achieve the maximum GFP expression is 10.1 h for day-1 and 16.7 h for day-2 (Figure 7b-d).Meanwhile, the head section of day-1 and day-2 cohorts have their maximum GFP expression on average at approximately the same time (9.7 h versus 10.3 h after heat-shock for day-1 and day-2, respectively).Besides a quicker HSR, day-1 nematodes also have a higher average GFP maximum expression across the compartments, except for the tail.Among the compartments, the early adulthood decline of HSR in the intestine is notably stronger than the other tissues.
We again used the mathematical model above to describe the compartment-specific HSR dynamics, following the same parameter estimation procedure aforementioned.As before, we found that the heterogeneity of GFP expression dynamics across the compartments can be explained by varying the model parameters associated with protein translation and degradation.Figure 7e depicts the rate constants of protein translation and degradation estimated using the average GFP dynamics of each compartment for day-1 and day-2 nematodes and the contours of the GFP slope (see Figure S10, Supporting Information, for the full parameter estimation of individual compartments in single worms).The central intestine section has turnover parameters that are consistently on the highest GFP slope in comparison to the other compartments, followed by the anterior intestine, the posterior intestine, the head, and lastly the tail.Comparing the parameters of day-1 and day-2 adults, we observed an age-related decline in the protein translation and degradation (protein turnover) and consequently in the GFP slope and thus the HSR capacity (Figure 7e).The decline in the GFP slope differs across compartments with those toward the middle of the animal having the larger drop, in agreement with the direct comparison of the average GFP fluorescence time profiles above.Notably, in comparison to day-1, the discrepancies of the GFP slopes among the compartments become less pronounced for day-2 cohort.In other words, the HSR capacity becomes more uniform across the compartments from day-1 to day-2.Note that while HSR capacity (dynamics) of various compartments in an individual worm become more alike with age, we observed again an increase in the heterogeneity of protein turnover parameters with age, consistent with the whole-body HSR modeling above (Table S4, Supporting Information).

Time-Lapse Imaging of C. elegans Embryos Reveals Maternal Contribution to HSR Heterogeneity
Given the above observation of proteostasis decline with age in adult nematodes and its effects on HSR capacity, we reasoned that proteostasis and HSR capacity in embryos should be reset to optimal efficiency and individual embryos' HSR should be similar to each other.To test this hypothesis, we quantified the GFP expression of embryos from TJ375 (Phsp-16.2::GFP).Briefly, we submitted gravid day-1, day-2, and day-4 C. elegans hermaphrodites to heat-shock for 1 h at 37 °C and immediately after collected their embryos.We placed the embryos in a microfluidic device for a longitudinal observation, as shown in  in Figure 8c, where the average responses of each dataset are indicated in bolder lines (magenta for embryos of day-1 adults, cyan for embryos of day-2 adults, and green for embryos of day-4 adults).
In comparison to adults, HSR dynamics in embryos are significantly faster, with an average GFP fluorescence peaking in 5 h after heat-shock.(Figure 8c).The faster HSR dynamics suggest a higher HSR capacity in embryos than in adults, confirming our hypothesis.That said, there is still heterogeneity in the HSR time-courses among individual embryos (Figure 8c).We applied our mathematical modeling to the embryo dataset, estimating the model parameters that reproduce the average response among the embryos (Table 4) and then finding the protein translation and degradation rate constants for individual embryos (Figure 8d, Figure S11, Table S5, Supporting Information).As for adults, variation of protein translation and degradation rate constants was able to capture the heterogeneity in the HSR dynamics among embryos.The embryos of day-2 mothers have faster protein turnover rate and thus a higher HSR capacity than those of day-1 mothers (k 6 : p-value = 4.2 × 10 −4 , k 11 : p-value = 2.4 × 10 −3 , two-sample t-test on mean difference).Also, in comparison to embryos of day-4 mothers, embryos of day-2 mothers have faster protein turnover and HSR dynamics (k 6 : p-value = 9.3 × 10 −3 , k 11 : p-value = 1.8 × 10 −2 , two-sample t-test on mean difference), suggesting that the HSR capacity in embryos is linked with the mother's age. [53]Finally, embryos from day-1 and day-2 mothers exhibited a similar level of heterogeneity in protein turnover rate constants k 6 and k 11 , with day-2 mothers producing the least heterogeneous embryos (two sample, two-tailed, t-test p-value < 0.003) (Table S5, Supporting Information).Meanwhile, day-4 mothers had the most heterogeneous embryos in protein turnover rate constants k 6 and k 11 (CV = 0.26 and 0.33 respectively).Thus, in the context of HSR capacity and its heterogeneity, day-2 mothers produced optimal embryos.

Discussion
Lifespan varies greatly across individuals, even among genetically identical animals reared in the same environment. [13]hile the heterogeneity of HSR and other stress response systems have been associated with individual longevity, [8][9][10][11][12][13] the underlying molecular determinants of such heterogeneity are still insufficiently understood.Using bespoke microfluidic systems that enables high-throughput and high-resolution measurements of HSR dynamics in TJ375 C. elegans and a molecular level mathematical modeling of the HSR, we identified important molecular aspects of HSR dynamics in C. elegans at single-organism level.First, individual TJ375 C. elegans display a significant non-genetically driven heterogeneity in their HSR dynamics, as indicated by GFP expression regulated by the HSE.Interestingly, cross-sectional screening of heat-shocked C. elegans populations at different time points showed that this heterogeneity is highly reproducible.Specifically, we found little change in the population heterogeneity of GFP expression in TJ375 as defined by the coefficient of variation [54,55] between 4 and 48 h after heat-shock (Figure 1e, Figure S3, Supporting Information).In eukaryotes, random fluctuations in chromatin states or in the number of ribosomes can be a major source of variation in gene expression. [54,55]However, altering the chromatin remodeler which is important for the HSR (i.e., jmjd-3.1 mutants; [23] Figure S4, Supporting Information) or using mutant backgrounds with higher numbers of ribosomes (i.e., ncl-1 mutants; Figure 5f) did not significantly affect the standard deviation nor the overall heterogeneity.
To identify the molecular underpinning of the observed heterogeneity in HSR dynamics, we applied mathematical modeling and sensitivity analysis to our high-resolution population and time-lapse HSR data of individual C. elegans.We carried out sensitivity analysis using an ensemble of model parameter sets (n = 10) that provided good fit to the average GFP fluorescence data, that is, the average worm (see Experimental Section).Aggregating results from multiple model parameter sets alleviated a potential issue that the outcome of the sensitivity analysis is an artifact of a specific model parameterization.The sensitivity analysis pointed to model parameters related to the maintenance of protein amount, that is, proteostasis, as the controlling factors of HSR dynamics.By changing model parameters associated with protein translation and degradation rates and their effects on total protein and HSF amount, while keeping all other parameters the same as the average worm, we were able to fit the fluorescence data of HSR dynamics in individual worms (Figure 4).Thus, our single worm model fitting suggested that the overt heterogeneity of HSR dynamics in the isogenic TJ375 C. elegans population can be explained by differences in proteostasis capacity among the animals.Considering the complexity of the HSR model in our study, there could be other parameter combinations that could also be used to fit the single-worm data.[20] Furthermore, we also performed validation experiments by crossing TJ375 with ncl-1 mutant strain that has an almost double general protein translation, to generate a transgenic C. elegans LSD2089 with altered protein turnover rates-analogous to doubling the protein translation rate constant.High-throughput time-series measurements of HSR in LSD2089 population after the same heat-shock stimulus showed that increasing protein translation capacity leads to a proportional rise in the HSR activity-that is, a higher slope of GFP-in agreement with the mathematical modeling.
Guided by model simulations of HSR dynamics in 10000 in silico worms, we sorted and isolated subpopulations of TJ375 C. elegans according to their HSR capacity.Lifespan analysis of the TJ375 C. elegans subpopulations demonstrated that the heterogeneity of HSR dynamics in early adulthood translates to variance of lifespan in these isogenic C. elegans reared in the same environment.More specifically, nematodes displaying a higher HSR capacity, presumably by having a more robust proteostasis, have a longer mean lifespan (Figure 6).Thus, our results point to the variability in protein turnover-the production and degradation of proteins-or more generally in protein homeostasis (proteostasis) fidelity, as an important driver of heterogeneity in the length of the C. elegans adult lifespan.
Our model-based finding, showing that differences in protein expression and turnover capacity among individual worms can explain the overt variability of HSR, agrees well with recent studies in C. elegans and HeLa cells. [18]Beyond confirming the results of past studies, our characterization of HSR dynamics of day-1 and day-2 adult TJ375 indicated a rapid age-dependent decline in the activity of the HSR resulting from a significant decrease in the protein turnover parameters (Figure 5c,d).Consistent with this trend, measurements of TJ375 C. elegans embryos show a much stronger HSR activity than in adult C. elegans (Figure 8c), suggesting a resetting of the proteostasis-tabula rasa-in embryos.Mathematical modeling of HSR dynamics in day-1 and day-2 C. elegans further showed that the degree of age-related proteostasis decline varies among individual C. elegans and even across compartments within an animal.The stochasticity in the age-related proteostasis decline has the effect of increasing animal-to-animal heterogeneity in the protein turnover and HSR dynamics across the whole-body as well as within individual compartments of the nematodes.Since the C. elegans in this study shared the same genome and were maintained in a uniform environment, the stochastic-onset of proteostasis collapse in early adulthood may have an epigenetic origin. [45]In contrast, the protein turnover parameters and HSR capacity along the body of a single animal become more alike with age.Taken together, the results point to a higher inter-animal heterogeneity, but a lower intraanimal variability of proteostasis and stress response capacity with aging.
[46] This collapse may be due to the allocation of resources to the growth of embryos as postulated by the disposable soma theory of aging.For instance, in early adulthood, C. elegans intestine starts producing yolk proteins in high amounts that are secreted and picked up by the developing embryos. [56]Hence, intestinal cells can re-allocate their resources for maintaining the intestinal proteome, which includes the expression of molecular chaperones, to production of yolk proteins.Interestingly, we observed the intestine compartment of individual TJ375 C. elegans experiencing the most prominent decline in HSR among the compartments as well as the highest increase in worm-to-worm variability with age (Tables S2 and S4, Supporting Information).Of note, intestinal hsp-16.2expression has previously been shown to correlate most with TJ375 lifespan. [20]We speculate that the investment in protein maintenance to produce yolk has a trade-off of increasing variability in other cellular functions, such as stress response.Interestingly, inhibiting genes in the yolk production and conversely overexpressing chaperone HSP-16 is both sufficient to increase lifespan of C. elegans, [57,58] suggesting that interfering with the expression of these genes might tip the scale towards somatic protein maintenance, resulting in systemic effects on organismal longevity.
In the context of evolution, it has been suggested that there exists a trade-off between developing and reproducing rapidly and the quality of the embryos. [53]For instance, the hatching size, growth rate, developmental speed, starvation resistance, and fecundity show heterogeneity stemming from the maternal age and packaging of the oocytes. [53]Our data and modeling of HSR in TJ375 embryos indicate that day-1 and day-4 mothers produce lower quality embryos than day-2 mothers, and that day-2 mothers give embryos with the highest average HSR capacity and the lowest HSR variability.The results thus give evidence for a maternal (transgenerational) contribution to the heterogeneity of HSR dynamics that is related to the mother's age, similar to other heterogeneous phenotypes of C. elegans. [53]While proteostasis appears to be reset in embryos, there still exists heterogeneity among the embryos, which may arise from the packaging of nutrients, mRNA, and proteins from the mother to the egg. [53]

Experimental Section
C. elegans Strains and Culture: C. elegans populations were cultured and maintained at 20 °C on nematode growth media (NGM) plates and fed with E. coli OP50 using standard protocols [59] before and between microfluidic experiments.Heat-shock was delivered to the animals by placing the NGM plates with C. elegans inside an incubator at 37 o C for 1 h.Age-synchronized day-1 and day-2 adults, via bleaching, were used in experiments.The timing between synchronization and experiments was tightly controlled for comparable results.For population and individual assays in microfluidic devices, animals were suspended in M9 solution before each experiment, repeatedly washed to ensure the complete removal of progeny and debris.
Microfluidic Device Fabrication: Microfluidic devices were fabricated using standard soft lithographic techniques and consisted of one structured polydimethylsiloxane layer bonded to a planar glass slide.Briefly, microfluidic circuits were designed using AutoCAD 2014 (Autodesk, San Rafael, USA) and printed onto a high-resolution film photomask (Micro Lithography Services Ltd, Chelmsford, UK).Master structures were fabricated on SU-8 (Microchem Corporation, Westborough, USA) coated silicon wafers using conventional photolithographic methods, [60] with feature heights being 95 µm for the "high-throughput" devices (for adult C. elegans) and 65 µm for the "longterm immobilization" devices (for adult C. elegans).Microfluidic devices were subsequently manufactured using standard soft-lithographic techniques. [61]Specifically, a 20:1 wt/wt mixture of polydimethylsiloxane base and Sylgard 184 curing agent (Dow Corning, Midland, USA) [62] was poured over the SU8 mold structure and cured in oven at 70 °C for 8 h.Such a composition ensured that the PDMS was flexible enough for efficient PDMS-valve actuation.The cured PDMS layer was then peeled off, cut to shape, and inlet/outlet was formed using a hole puncher.Finally, the structured PDMS layer was bonded to a planar glass slide (Menzel Gläser, Braunschweig, Germany) in an oxygen plasma, and the bonding was reinforced by a post-bake at 70 o C for 3 h.
High-Throughput C. elegans Screening and Sorting Microfluidic Devices: Two microfluidic devices for high-throughput screening and sorting of C. elegans were developed based on their fluorescent protein reporter expression.The layouts of the devices are depicted in Figures 1e and 2a, respectively.They consist of a central straight channel, where C. elegans (suspended in M9 buffer) flow in a continuous and non-overlapping manner.The nematodes enter the device through an inlet, which has wall-embedded PDMS blades and a central pillar to aid disaggregation of animals.Two additional buffer inlets enable appropriate spacing between single C. elegans flowing through the main channel.Whereas the screening device has only one outlet (Figure 1e), the sorting device has three outlets (Figure 2a), emanating from two consecutive bifurcations, subsequent to the main channel.In the sorting device (Figure 2a), during normal operation, one outlet was kept open to allow fluid and non-sorted C. elegans (i.e., animals expressing GFP levels lower than the minimum threshold value) to exit the device.Upon production of a sorting signal, the total flow was relocated through a specific outlet via actuation of lateral on-chip valves.These on-chip valves were formed by dead-end fluidic channels (Figure 2a, in pink) filled with deionized water and connected to high-pressure sources.Upon actuation, these valves acutely reduce the diameter (by approximately 95%) of the fluidic channel adjacent to the valve, increasing its hydraulic resistance and deflecting flow towards the remaining "open" channel.Finally, it should be noted that a reference channel was used for calibration, where fluorescence intensity originating from a standard dilution of fluorescein sodium salt (Sigma Aldrich, Buchs, Switzerland) was measured at a constant flow rate.
Efficient operation and performance of the high-throughput microfluidic devices were enabled by the functional integration of peripheral components of the system.These components are depicted in Figure S2a,b, Supporting Information.Continuous fluid flow was provided by precision syringe pumps (neMESYS, CETONI GmbH, Korbussen, Germany).Two dosing units delivered M9 buffer to the designated inlets of the device at volumetric flow rates of 50 and 70 μLmin −1 respectively, while the suspension of C. elegans was controlled by other two-fluid units (which delivered them at 80 μLmin −1 ).The optical system (Figure S2a,b, Supporting Information) provided a rapid detection of fluorescent protein reporters within animals.The optical system comprised a laser light source, an inverted confocal microscope, and a photomultiplier tube (PMT) detector (Figure S2b, Supporting Information).For accurate C. elegans GFP protein quantification, a cylindrical lens was used to shape the laser beam into a light sheet orthogonal to the direction of C. elegans flow (an example of fluorescent signal readout is shown in Figure S2e, Supporting Information).A 488 nm laser beam (60 mW, Omicron Laser, Rodgau, Germany) was used as the excitation source and expanded using a telescope consisting of two plano-convex lenses and shaped into a light sheet using a cylindrical lens (LJ1558RM, f = 300 mm, Thorlabs, Newton NJ, USA).The light sheet was then directed into the back aperture of the microscope, using a dichroic mirror (AT DC 505, AHF, Tübingen, Germany) and subsequently focused into the microfluidic channel using a 20× objective (S Plan Fluor ELWD ADM, Nikon, Egg, Switzerland).Fluorescence emission was collected by the objective, passed through the same dichroic mirror, and focused onto a 40 ± 3 µm pinhole (P40H, Thorlabs, Newton NJ, USA) through a tube lens (f = 200, Nikon, Egg, Switzerland).Subsequently, fluorescence was passed through an emission filter (525/45 nm, HC Bright Line Bandpass Filter, AHF, Tübingen, Germany) and onto a PMT (H10722520, Hamamatsu Photonics, Solothurn, Switzerland).A secondary detection system was used for visual inspection of device performance.Specifically, white light, bright field illumination filtered with a 610 nm emission filter (ET Longpass, AHF, Tübingen, Germany) was imaged using a high-speed camera (Phantom Miro, Vision Research, Baden-Baden, Germany).Using such a system, real-time images could be acquired in parallel without this (longer) wavelength light interfering with the fluorescence measurement.Collected light passed through the microscope beam splitter 80/20, with 20% of the beam being directed to the high-speed camera, which was protected from the laser beam with a 532 nm longpass filter (LP Edge Basic, AHF, Tübingen, Germany) and 80% being directed to the PMT.Sorting decisions were made by an FPGA (PCI NI-RIO, National Instruments, Austin, USA) (Figure S2a, Supporting Information) enabled system, and animal sorting was enabled by valve actuation, which promptly changes the fluidic resistance at given positions within the fluidic path, and directs the flow into a specific outlet.On-chip valves were connected to solenoid valves (MHA1, Festo, Esslingen, Germany), which were controlled by sorting impulses delivered by the FPGA.Valve-actuation times were defined by decision-making and the signal-output of the FPGA, the switch of configuration of the solenoid valve (approximately 4 ms), and the inflation time of the PDMS valve (between 8 and 15 ms).The total processing time for one C. elegans (from the fluorescence detection to the sorting in the proper outlet, enabled by actuation of valves) is 75 ms.
Long-Term C. elegans Immobilization Microfluidic Device: The longterm immobilization microfluidic device consists of a main channel, whose size and shape allow efficient, yet gentle, immobilization of C. elegans (Figure 3a, in red). [21]Two on-chip hydraulic valves positioned alongside the main channel allow control over nematode loading and immobilization (Figure 3a, in black).Hydraulic valves were formed by filling dead-end channels with water.Once pressurized, the wall separating the fluidic channel and valve channel deforms resulting in partial blockage of the fluidic channel.In devices designed for immobilization of adult hermaphrodite C. elegans, egg laying was facilitated through a set of pillars roughly positioned alongside the animal's vulva.A highly concentrated bacteria suspension was supplied to an immobilized C. elegans as food source (Figure 3a, in yellow).C. elegans immobilized in such a device remain viable for more than 100 h on average and exhibit normal physiological functions, with both feeding and egg laying occurring at rates comparable to free crawling C. elegans on NGM plates (Figure 3b).Furthermore, C. elegans were immobilized such that high-resolution fluorescence images can be acquired.These capabilities allowed quantitative long-term monitoring of the GFP expression after heat-shock in a non-invasive manner for individual animals.
All images were acquired on an inverted microscope (Nikon Ti-S, Nikon Inc.) equipped with an sCMOS camera (Prime 95B, Photometrics Inc.), a fluorescence LED (LedHUB, Omicron-laserage Laserprodukte GmbH) and a piezo objective drive (Nano-F100S, Mad City Labs Inc.), controlled by a custom build Matlab script (Matlab 2016b.Mathworks Inc.).Images of the immobilized C. elegans were acquired at 1-h intervals, at each time-point a fluorescent Z-stack was acquired (26 slices with 2 µm spacing) using an exposure time of 10 ms and 50% excitation intensity.All images were acquired using a 10× objective (CFI Plan Fluor 10X, Nikon Inc.).The mean fluorescence intensity was quantified from each C. elegans at each time point-measuring a defined region of interest which comprises the whole animal-with the open source software package ImageJ (NIH, Bethesda, USA).
Long-Term C. elegans Embryo Microfluidic Device: The long-term C. elegans embryo microfluidic device consisted of a simple straight channel within which embryos were confined for long-term imaging (Figure 8a).Such a channel was designed with an initial width of approximately 50 µm, so that single embryos could flow into the device and would remain on-chip once loading was completed.Following the loading region, the channel expanded to 200 µm width, here the channel width was chosen such that the channel would fit into the microscope field of view at 40-60× magnification.An array of small pillars was placed into the wide channel, trapping embryos within the channel.Channel height was chosen at 30 µm, such that most of the embryos would be oriented for optimal imaging (Figure 8b).
Embryos were harvested from adult hermaphrodite animals, at different ages (Adult day-1, day-2, and day-4).Following heat-shock at 37 °C for 1-h, adult animals were synchronized via bleaching (200 µl of 1 M NaOH and 400 µl of 10% NaClO), the embryos collected and washed twice with M9 buffer.Embryos were passed through a 40 µm cell strainer to remove any debris, loaded into a short piece of tubing, and connected to the device.Prior to loading, the microfluidic device was filled with M9 buffer and gently pressurized to remove all air from the device.Embryos were finally loaded into the device by gently applying pressure to the syringe attached to the embryo inlet, and passively trapped on the pillar array placed into the wide imaging channel.30 Z-stacks of the embryos throughout the HSR were recorded on an epi-fluorescence microscope with a spacing of 1 µm using a 40× oil immersion objective.Stacks were acquired at 10-min intervals, for a total duration of 12 h.The mean fluorescence intensity was quantified from each C. elegans's embryo at each time point-measuring specific regions of interest which comprised a single embryo-with the open source software package ImageJ (NIH, Bethesda, USA) (Figure 8c).
Western Blot Quantification of GFP and HSP-16.2:Western blot quantification of native GFP and HSP-16.2 was carried out in strain TJ375 in order to study their relative expression at different time points after heat-shock and validate our measurements.For this, samples of C. elegans strain TJ375 were harvested at different time points (4, 8, 12, 18, 24, 48 h) after heat-shock, as well as TJ375 samples without heat-shock.N2 control samples were harvested 8 h after heat-shock and without heat-shock (Figure 1c,d and Figure S1, Supporting Information).Around 2000 C. elegans per sample were placed in M9, frozen immediately in liquid nitrogen and maintained at −80 °C.Protein extraction from C. elegans samples was done by adding 150 µl lysis buffer (RIPA buffer (ThermoFisher #89900, Reinach, Switzerland), 20 mM sodium fluoride (Sigma #67414, Buchs, Switzerland), 2 mM sodium orthovanadate (Sigma #450243, Buchs, Switzerland), and protease inhibitor (Roche #04693116001, Basel, Switzerland)) to each C. elegans pellet and subsequently sonicating them repetitively until tissues were dissolved.Aliquots were collected and processed.For equal loading, the protein concentration of the supernatant was determined with BioRad DC protein assay kit II (BioRad #5000116, Cressier, Switzerland) and standard curve with Albumin (ThermoFisher #23210, Reinach, Switzerland).After 5 min of boiling at 95 o C, 10 µg of protein in 20 µl of loading buffer (42.5 mM Tris-HCl pH 6.8, 1.7% SDS, 8.5% glycerol, 15% 2-mercaptoethanol, and 1% bromophenol blue, final concentrations) were loaded per sample in each well of a NuPAGE 10% Bis-Tris protein gel (ThermoFisher #NP0301BOX, Reinach, Switzerland).Gel electrophoresis was done for 1.5 h at 130 V using a Mini Gel Tank filled with MOPS 1× buffer.Right after, the gel was electro-blotted onto a nitrocellulose membrane (Sigma #GE10600002, Buchs, Switzerland) using a Mini Trans-Blot Cell filled with transfer buffer (25 mM Tris-Base, 0.2 M glycine, and 20% methanol, final concentrations) and applying 100 V for 1 h.Subsequently, the membrane was blocked for 1 h using 5% milk powder in TBS-t buffer (50 mM Tris-Base, 0.15 M NaCl, and 0.1% Tween-20, pH 7.2-7.6)at room temperature.Right after, the membrane was cut (according to bands of molecular weight indicated by the protein ladder) in 3 pieces for identification of HSP-16.2,GFP, and Tubulin respectively.The membrane for identification of HSP-16.2 was incubated in 5% milk powder in TBS-t buffer with a polyclonal anti-HSP-16.2primary antibody (1:8000) courtesy of Prof. Gordon Lithgow.The membrane for identification of GFP was incubated in 5% milk powder in PBS-t buffer (0.13 M NaCl, 2.7 mM KCl, 10 mM Na 2 HPO 4 , 1.8 mM KH 2 PO 4 , and 0.2% Tween-20, pH 7.2 -7.6) with anti-GFP primary antibody (1:32 000, Roche #11814460001, Basel, Switzerland).Lastly, the membrane for identification of Tubulin, as control protein, was incubated in 5% milk powder in TBS-t buffer with anti-tubulin primary antibody (1:500, Sigma #T9026, Buchs, Switzerland)).All membranes with primary antibodies were incubated overnight at 4 o C.Then, after thorough TBS washes, blots were probed with secondary antibodies: Tubulin and GFP with anti-mouse IgG antibody (1:2000) and HSP-16.2 with anti-rabbit IgG antibody (1:1500) in TBS-t suspensions with 5% milk powder for 1 and 1.5 h respectively.Finally, the protein bands were visualized upon few minutes of incubation of the membranes with Bio-Rad Clarity Western ECL Substrate (Bio-Rad #1705061, Cressier, Switzerland).Blots were imaged using ChemiDoc imaging system for chemiluminescence assays with 0.1 s exposure.Density of the bands was quantified using software package ImageJ (NIH, Bethesda, USA), and normalized to Tubulin values (Figure 1c,d and Figure S1, Supporting Information).
HSR Mathematical Modeling: The mathematical model of HSR consists of 12 species (Table 1), 15 reactions (Listed as Equations (1-15) in Figure 4a & Table 2)-five of which were reversible-with 24 parameters.The full model derivation of the Equations (1-15) is provided in Supporting Information S2.The model simulates the cellular response to an increase in intracellular MFP caused by a heat-shock.Prior to the simulation of HSR, a pre-equilibration step was performed by running the model until steady-state (i.e., until the concentrations stabilize).Heat-shock stress was simulated by increasing the misfolding rate of proteins (Prot) from α to β (β > α) for a period of 1 h, emulating the transfer of C. elegans from a 20 to a 37 o C environment.In the model, MFPs were captured by molecular chaperons, specifically the heat-shock protein (HSP), for protein refolding.As a consequence of the heat-shock, Heat-shock Factor 1 (HSF-1) molecules that would otherwise form a complex with HSP (HSP:HSF-1) become free.Free HSF-1 molecules subsequently undergo reversible dimerization (HSF-1 2 ) and trimerization (HSF-1 3 ).HSF-1 trimers were able to translocate to the nucleus, bind to the HSE, and promote the expression of the HSP.The newly produced HSP (in excess of the level needed to refold MFP) will bind to free HSF-1 molecules and thereby repress its own expression.For the TJ375 strain, GFP was encoded downstream of a promoter containing an HSE.Accordingly, the original HSR model was modified to include a balance equation for the GFP (Equation ( 15)), in which the expression rate was the same as that of HSP.As a reporter, the GFP does not interact with any other species in the model, and thus does not influence the HSR dynamics.It should be noted that the model does not include a balance equation for the protein and HSF-1, and thus the total amount of protein (Prot + MFP) and HSF-1 does not vary with time.However, the amount of protein (Prot) and HSF-1 was allowed to differ among individual C. elegans (see below).
Parameter Estimation for the Average HSR Dynamics of Adult C. elegans: In the first stage of parameter estimation, the model parameter estimates were obtained by fitting the simulated GFP concentration to the average fluorescence values amongst all measured C. elegans.When comparing model prediction of GFP expression and the average fluorescence data, the value of the simulated GFP concentration was scaled such that its peak equals the maximum value of the average fluorescence intensity.The parameter estimation was formulated as a minimization of the objective function: where p is the model parameter vector, n d is the number of measurement data (n d = 48), y j is the average GFP fluorescence intensity among the C. elegans at time j, and GFP(p,t j ) is the simulated GFP concentration at time j for the parameter vector p The lower and upper bounds of the parameter p, denoted by lb k and ub k (i.e., lb k ≤ p k ≤ ub k ) for the k-th parameter, were set to 2 orders of magnitude below and above those of the original model parameter values reported in Scheff et al. [34] The optimization was solved using a hybrid scatter search algorithm called MEIGO in MATLAB. [63]The resulting parameter estimates are provided in Table 3. Sensitivity Analysis: A local sensitivity analysis of the model was performed by perturbing each parameter one at a time and simulating the resulting GFP expression curve.The sensitivity coefficients of the HSR dynamics were calculated for the following: the maximum value of GFP and the time to reach the GFP maximum value.In this manner, the parameters that exert a strong influence on the slope of the GFP expression curve were identified.The sensitivity coefficients were computed using a finite difference approximation, as follows: x pe x p i j i j i j j i j (17)   where s i,j is the sensitivity coefficient of x i with respect to perturbation to the j-th parameter, x i is one of the two metrics of HSR dynamics above, Δp j denotes the perturbation to the j-th model parameter-chosen to be 10% of the nominal parameter values-and e j is the unit basis vector in the j-th direction.The sensitivity analysis for the model parameters was performed in Table 3 and in those parameter combinations shown in Table S7, Supporting Information.The average sensitivity coefficients are presented in Figure S6, Supporting Information, which are presented in descendent order according to their absolute magnitude across the top 10 parameter combinations.The sensitivity analysis suggests that the parameters related to the HSP (GFP) transcription and translation, specifically, k 11 , k 6 , and k 4 , the amount of proteins HSF 0 and HFE 0 have a strong influence on the GFP dynamics before its peak.HSR Dynamics in Single C. elegans: The HSR model of the average C. elegans above was further adapted to describe the HSR time course from single C. elegans.The ODE model formulation was kept the same as described above, and further assumed that individual C. elegans differ from the average C. elegans only for a few model parameters.More specifically, for the single C. elegans parameter fitting, parameters related to biologically regulated process were only considered such as the protein turnover, and set the other physicochemical parameters such as protein misfolding rates and HSF-1 dimerization/trimerization to be constant across different animals (see Table 3).The results from the model sensitivity analysis above point to the regulation of protein amount as the most important process that modulates the HSR dynamics.A recent study of HSR in HeLa cells by Guilbert et al. [18] demonstrated that basal protein expression variability regulates the overt HSR heterogeneity.
With these considerations, the possible model parameters that vary among animals were narrowed down to the top 3 parameters based on the sensitivity magnitudes: protein translation rate constant k 6 , protein degradation rate constant k 11 , and the initial amount of proteins HSF 0 .Note that since the model does not include new synthesis and degradation of any forms of HSF by assuming that they exist at pseudo steady-state level, the initial amount of HSF 0 thus reflects the total amount at any time.Among the three parameters, two were related to protein turnover.Thus, in the modeling of HSR dynamics in single animals, the protein translation and degradation rate constants were varied for each C. elegans by scaling the translation and degradation rates of HSP.Note that the initial (total) amount of proteins HSF 0 and Prot 0 were scaled based on protein turnover rate constants.More specifically, the total (initial) amounts of proteins HSF where the subscript l refers to the parameters from the specific C. elegans l and the non-subscripted parameters were taken from the average C. elegans estimates in Table 3.The single C. elegans parameter fitting was performed in the same manner as for the average C. elegans except for the fluorescence data, which now came from individual C. elegans: where n d is the number of measurement data (n d = 48), y l,j is the measured GFP intensity of C. elegans l at time j and GFP l is the simulated GFP concentration at time j for that specific C. elegans l.The upper and lower bound of the parameters were allowed to vary by 1 order of magnitude from the values in Table 3.This optimization was performed using the MATLAB function fmincon with multiple starts (n = 120).The results obtained from the parameter estimation of individual C. elegans fluorescence data are shown in Figure 4d.HSR Dynamics in Embryo C. elegans: Finally, the mathematical model was expanded for describing HSR in embryo C. elegans.Keeping the ODE model formulation previously described, the time-course HSR experimental data obtained from embryo C. elegans was used.The model parameter estimates were first obtained by fitting the simulated GFP concentration to the average fluorescence values amongst all measured C. elegans eggs, shown in Figure 8c.The optimization was as well solved using MEIGO in MATLAB, [63] by iterating 1000 times the calculation, starting from the values shown in Table 3, which describe HSR in adult C. elegans.The result is shown in Table 4.
Subsequently, the individual translation and degradation constant rates (k 6 and k 11 ) were found which would describe the individual time course HSR of individual embryos.For this, the procedure described in the previous section (HSR dynamics in single C. elegans) was reproduced.The resulting parameter estimates for individual embryo C. elegans are shown in Figure 8d and Table S5, Supporting Information.
Lifespan Machine: The automated lifespan analysis was performed using the lifespan machine platform described by. [50]The experiment was prepared identically to the manual lifespan setup with the exception of using tight-fitting petri dishes (BD Falcon Petri Dishes, 50×9mm) and drying the plates without their lids for 30 min prior to starting the experiment.The automated analysis was performed in aircooled Epson V800 scanners using a scanning frequency of one scan per 30 min.The temperature was monitored using temperature probes (Thermoworks, Utah, US) on the scanner flatbed and kept constant at 20 °C.
Manually Performed Lifespan (by Hand): Manual lifespan scoring (by hand) was performed as described in Ref. [64].In brief, about 2000 TJ375 animals were heat-shocked and sorted as described above and then about 100 L4 C. elegans per sorting category were picked on NGM plates containing 50 µM 5-Fluoro-2′deoxyuridine (FUDR), 15 µgml −1 nystatin, and 100 µgml −1 carbenicillin seeded with heat-killed OP50 bacteria [65] and scored at 20 °C.Animals were classified as dead if they failed to respond to prodding.Exploded, bagged, burrowed, or animals that left the agar were censored.Population survival was analyzed with the statistical software R [66] using the survival [67] and survminer [68] packages and with the L4 stage defined as time point zero.
The estimates of survival functions were calculated using the productlimit (Kaplan-Meier) method and the log-rank (Mantel-Cox) method was used to test the null hypothesis.

Figure 1 .
Figure 1.Longitudinal C. elegans experimental HSR characterization at population level.a) Outline of the gene mutation of strain TJ375 [Phsp-16.2::GFP(gpIs1)].[15]b) Schematics of the GFP expression on the TJ375 C. elegans strain following heat-shock.GFP is only apparent on C. elegans hours after heat-shock.c) Schematic of the microfluidic device for C. elegans fluorescent protein screening.The "C. elegans inlet" i) incorporates wallembedded PDMS blades to untangle animals as they enter the device.Additional buffer inlets ii) enable efficient spacing of C. elegans in the main channel.iii) Fluorescence detection volume where the optical detection system is focused onto.Animals exit via the device outlet (iv) and calibration is performed using the reference channel (v).The flow direction is indicated by the black arrow.Scale bars are 2.5 mm.d) Experimental layout of C. elegans population screening experiments.e) Violin plots of the fluorescence mean values of an isogenic population of C. elegans TJ375 (n = 2280) measured at determined time points (2, 4, 8, 12, 24, 32, and 48 h) after heat-shock, as the experimental layout shown in (d).The thicker black lines inside the violins represent 50 percent of the population.
Figure 1.Longitudinal C. elegans experimental HSR characterization at population level.a) Outline of the gene mutation of strain TJ375 [Phsp-16.2::GFP(gpIs1)].[15]b) Schematics of the GFP expression on the TJ375 C. elegans strain following heat-shock.GFP is only apparent on C. elegans hours after heat-shock.c) Schematic of the microfluidic device for C. elegans fluorescent protein screening.The "C. elegans inlet" i) incorporates wallembedded PDMS blades to untangle animals as they enter the device.Additional buffer inlets ii) enable efficient spacing of C. elegans in the main channel.iii) Fluorescence detection volume where the optical detection system is focused onto.Animals exit via the device outlet (iv) and calibration is performed using the reference channel (v).The flow direction is indicated by the black arrow.Scale bars are 2.5 mm.d) Experimental layout of C. elegans population screening experiments.e) Violin plots of the fluorescence mean values of an isogenic population of C. elegans TJ375 (n = 2280) measured at determined time points (2, 4, 8, 12, 24, 32, and 48 h) after heat-shock, as the experimental layout shown in (d).The thicker black lines inside the violins represent 50 percent of the population.

Figure 2 .
Figure 2. C. elegans experimental HSR characterization at sub-population level.a) Schematic of the microfluidic device for C. elegans sorting.The fluidic layers are shown in green and the control layer in pink.The flow direction is indicated by the black arrow.Scale bar is 2.5 mm.The "C. elegans inlet" i) incorporates wall-embedded PDMS blades to untangle animals as they enter the device.Additional buffer inlets ii) enable efficient spacing of animals in the main channel.iii) Fluorescence detection volume where the optical detection system is focused onto.Sorted animals exit via one of the device outlets (iv) and calibration is performed using the reference channel (v).Lateral on-chip valves are labeled (vi).b) Experimental layout of C. elegans population sorting experiments.c) Double sorting and subsequent GFP intensity measurement of a C. elegans population (n = 2600), following layout shown in (b), at different time points after heat-shock.Each measured population is shown with violin plots.Initial sorting is performed 4 h after heatshock, dividing the population into a top and bottom 50%.A secondary sorting step (8 h after heat-shock) further separated the top 50% into its top 16% and bottom 84%.Similarly, the initial bottom 50% of the population is further sorted into a bottom 16% and a top 84%.The four resultant groups are screened at posterior time points (12, 25, and 32 h).The change in the mean value of each group at each sorting or screening time is shown with a distinct color line per group, that is, green: "high-high" group, blue: "high-low", yellow: "low-high", and red: "low-low".

Figure 3 .
Figure 3. C. elegans experimental HSR characterization at individual level.a) Schematic of the microfluidic design that allows for adult C. elegans immobilization.Arrows represent the flow direction.Control layer is represented in black.b) Images of an immobilized animal at determined time points after trapping.Scale bar represents 100 µm.c) Increasing time-lapse fluorescence images from an immobilized day-2 adult TJ375 C. elegans, after 1-h heat-shock induction, in microfluidic device along 48 h, where t = 1 h is at the furthest left and t = 48 h is at the furthest right.Scale bar = 200 µm.d) Measured GFP mean expression over time (t = 48 h) in individual TJ375 C. elegans.Each curve represents a single C. elegans (n = 15) and they are identified with a number next to the curve.The mean response of the dataset is shown in a bold black line.

Figure 4 .
Figure 4. C. elegans in silico HSR mathematical modeling.a) Reaction network describing the C. elegans HSR in mutants with GFP as protein reporter of HSP.Solid arrows indicate reaction paths, whereas dotted arrows indicate an activation where the substrate is not consumed.The number annotations correspond to the listed reactions in Table 2. b) Dynamic expression of the primary species determining HSR: MFP, Trimer HSF-1 bound to HSE (HSF 3 :HSE), mRNA, GFP, and HSP-16.2.Heat-shock takes place between -1 h and 0 h.Positive values on the time axis indicate the time in hours after heat-shock.c) Individual fitting of individual GFP expression for each analyzed animal, by varying only parameters k 6 and k 11 from the average response.Here the solid lines represent the modeled HSR while the circles represent the experimental data.d) k 6 -k 11 parameter space found for the experimental dataset shown in Figure 3d.The number next to each data point identifies the C. elegans analyzed.

Figure 5 .
Figure 5. HSR study for C. elegans with different ages (day-1 versus day-2 adults).a) Violin plots of the fluorescence mean values of individual C. elegans measured in an isogenic cohort at day-1 adulthood (n = 6917) at determined time points (4, 8, 12, 24, and 28 h) after heat-shock.The change in the mean value of each group at each sorting or screening time is shown with the magenta line.b) Homologous violin plots of the fluorescent mean values of individual C. elegans measured in an isogenic cohort at day-2 adulthood (n = 2716).Please note that the arbitrary units used in both (a) and (b) are not equivalent.c) GFP mean expression in single animals over time (t = 30 h): In magenta day-1 adult C. elegans (n = 10) and in cyan day-2 adult C. elegans (n = 15).d) k 6 -k 11 parameter space found for the experimental datasets shown in (c).The number next to each data point identifies the C. elegans analyzed (See Figure3dand FigureS7, Supporting Information), where magenta dots correspond to day-1 adults and cyan dots to day-2 adults.e) Contour plot of GFP expression slope values according to k 6 and k 11 values (translation and degradation rates of GFP), where dark blue represents lower values and yellow higher values.Steeper slope isolines correlate with higher translation rates and lower degradation rates.Overlaid are the measured parameters shown in (d).f) GFP expression curves of LSD2089 (ncl-1, Phsp-16.2::GFP)(dotted lines) and TJ375 (Phsp-16.2::GFP)(solid lines) for adult day-1 (in magenta) and adult day-2 (in cyan) after heat-shock.In each population, for each time point, the mean fluorescence value and its standard deviation are shown, when measured with the high-throughput microfluidic system.

Figure 6 .
Figure 6.Sorting and lifespan studies based on HSR difference.a) In silico k 6 -k 11 parameter space (n = 10 000) created around the experimental data from the individual HSR characterization (n = 25) shown in Figure 5d.b) Simulated profiles (n = 10 000) of GFP expression after heat-shock, where each curve is generated with a unique combination of k 6 and k 11 , from (a). c) Parameter subspaces for groups "high" and "low" which correspond to GFP expression curves with different slopes, when selecting top and bottom 5% of the GFP expression, 5 h after heat-shock.Overlaid it is shown the contour plot of Figure 5e.d) Selected GFP expression profiles with the sorting criteria top and bottom 5% GFP expression at 5 h after heat-shock.e) Lifespan curves resulting from the lifespan assays of the sorted C. elegans groups (combined results from 3 biological repeats).The population statistics are enlisted on the right.

Figure 7 .
Figure 7. Differences in tissue-specific HSR in C. elegans.a) Fluorescence image of an immobilized day-2 adult TJ375 C. elegans in a microfluidic device for long-term immobilization; the segmentation for tissue GFP expression analysis are indicated with yellow boxes, namely head, intestine anterior, intestine central, intestine posterior, and tail.Scale bar = 50 µm.b) GFP mean expression over time (t = 30 h) of different sections of the body of adult day-1 C. elegans after heat-shock.c) GFP mean expression over time (t = 30 h) of different sections of the body of adult day-2 C. elegans after heat-shock.d) Average GFP mean expression over time (t = 30 h) for day-1 (in magenta) and day-2 (in cyan) adult C. elegans of different sections of the body.Error lines represent the standard deviation of the dataset.e) Parameter space of individual sections of C. elegans adult day-1 (in magenta) and day-2 (in cyan).Overlaid are the isolines of GFP expression slope shown in Figure 5e.
Figure 8a,b (described in the Experimental Section section) and recorded their GFP fluorescence intensities every 10 min.The GFP fluorescence time-courses for the embryos are shown

Figure 8 .
Figure 8. HSR in embryonic C. elegans.a) Schematic of the microfluidic design that allows for long-term embryonic C. elegans observation.Arrows represent the flow direction.Scale bar represents 1 mm and 100 µm.b) Bright-field image of a collection of embryos in the microfluidic chamber filled with embryos.Scale bar represents 30 µm.c) GFP expression after heat-shock of embryos coming from different age gravid C. elegans.Magenta lines represent eggs from day-1 adults, cyan lines eggs from day-2 adults, and green lines eggs from day-4 adults; over time (t = 10 h).Bolder lines represent the average response of each dataset.d) k 6 -k 11 parameter space of the embryos shown in (c).Dots in magenta represent embryos from adult day-1, in cyan embryos from adult day-2, and in green embryos from adult day-4.
0 and Prot 0 in individual C. elegans vary with the parameters k 6 and k 11 as follows: (see derivations in Supporting Information S2)

Table 1 .
Molecular species in the HSR model.

Table 2 .
List of reactions in the HSR model.

Table 3 .
Model parameters estimated for average adult GFP expression data.

Table 4 .
Model parameters estimated for average embryo GFP expression data.