Stanislas Larnier
email: stanislas.larnier@laas.fr

Structure tensor based analysis of touching grain kernels for visual grading

Keywords: structure tensor, orientation field, visual grading, image analysis

Automatic visual grading of seed lots with a high density of touching grain kernels is a challenging problem. The structure tensor is a simple and robust descriptor that was developed to analyze textures orientation. Contrarily to segmentation methods which rely on an object based modelling of images, the structure tensor views the sample at a macroscopic scale, like a continuum. Thanks to this tool, it is possible to extract useful information on the orientation of grain kernels even in a bulk. Knowing the rough orientation of a grain kernel could provide an initialization for segmentation techniques. The comparative results with ground truth orientations on four different grain kernels demonstrate the ability of the structure tensor to detect the seed orientations.

Introduction

In agriculture, the global grain harvest reached several billion tons each year. Seed producers exchange their crops at a price determined by the quality of their production. This assessment, called grading, is performed for each set on a representative sample. The difficulty of this assessment is to fully characterize the sample. To do so, it is necessary to qualify each of its elements. Historically, this has been performed manually by an operator. This method is exposed to various problems and the results can vary from one operator to another.

Due to the huge amount of data contained in large images, automatic procedures are indeed essential to assess grain kernel shape properties such as size, orientation, aspect ratio, ... The lack of robust, fast and universal procedures to provide such a geometric description in the case of touching grain kernels is one of the obstacles to exploit the full potential of images.

The visual grading of grain kernels has been widely discussed during the last ten years. A number of studies have been conducted to classify grain based on physical attributes. Some classification studies are done on images with separated grain kernels [START_REF] Kaur | Classification and grading rice using multi-class svm[END_REF][START_REF] Agustin | Automatic milled rice quality analysis[END_REF]. Other studies have been conducted on the grain segmentation with low density of touching grain kernels [START_REF] Zhang | Separation of touching grain kernels in an image by ellipse fitting algorithm[END_REF][START_REF] Yao | An automatic segmentation algorithm for touching rice grains images[END_REF][START_REF] Faessel | Touching grain kernels separation by gap-filling[END_REF][START_REF] Mebatsion | A fourier analysis based algorithm to separate touching kernels in digital images[END_REF]. Segmentation processes that deal with touching grain kernels, with low or high densities, under the constraint of no seeds overlapping, have been lately proposed in [START_REF] Dubosclard | Automatic method for visual grading of seed food products[END_REF] and [8].

The structure tensor is an interesting preprocessing tool which could improve these segmentation methods in the case of high density of touching grain kernels. In the presented experiments, the considered seed samples are placed in bulk, without any arrangement. They are spread over a tray in such way that there is nearly no overlapping between the objects to avoid occlusion.

The structure tensor appeared in the field of image processing in the late 80's [9] for the problem of interest point detection. It was then justified theoretically and popularized in different contexts such as interest point detection [9], [10], texture analysis [11], representation of flow-like images [12], optical flow problems [13] and anisotropic or coherence enhancing diffusion [14]. The results of this paper demonstrate that the structure tensor provides a fast, robust and precise enough tool to retrieve grain kernel orientations. The proposed method has the avantage to be fast and invariant under contrast changes.

In Section 2, the definition of the structure tensor is recalled and some implementation information is provided. Section 3 proposes a way to choose the correct algorithm parameters. Section 4 demonstrates the interest of using the structure tensor on real images with low and high density of touching grain kernels. Its output is compared to the ground-truth orientations obtained by an human operator. These comparisons show that the structure tensor provides a good orientation estimation which could be used in segmentation approaches.

Structure tensor

The structure tensor of an image u is defined as follows

T := G ρ ∇u σ ∇u T σ , (1)
where u σ = G σ u, u denotes a gray-level image G σ and G ρ are Gaussian filters with standard deviation σ ≥ 0 and ρ ≥ 0, which control the amount of pre-and post-smoothing. The parameter σ is called the inner scale, the parameter ρ the outer scale. For a color image, the structure tensor is computed as the sum of the structure tensors for each channel.

In practice, the discrete gradient operator ∇ is defined by

∇ = (∂ 1 , ∂ 2 )
T . In all reported experiments, the partial differential operators ∂ i are defined as convolutions with discrete kernels. From an asymptotic point of view any kernel leading to a consistent discretization should provide good results. However, the kernel design turns out to be crucial to provide good orientation and anisotropy estimates. Key properties of discrete kernels are [15]: i) rotation invariance, ensuring a reliable orientation estimation, ii) separability, ensuring faster computations and iii) no shift, implying the use of centered finite differences. Following these criteria, the authors of [15] suggested to use the following operations: The gradient is computed in the space domain, while the convolution with G σ or G ρ are based on fast Fourier transforms. The overall computational complexity for an image with n pixels is therefore O(nlog(n)). In practice, the structure tensor can be computed in near real-time. In each pixel, the orientation is given by the direction of the eigenvector associated with the highest eigenvalue of the structure tensor.

∂ x u = 1 32   -3 0 3 -10 0 10 -3 0 3   u and ∂ y u = 1 32   -3 -10 -3 0 0 0 3 10 3   u.

Postsmoothing parameter estimation

As seen in the previous section, there is two parameters σ the inner scale and ρ the outer scale for the structure tensor. The experiments show that σ = 0 is sufficient, due to the quality of the acquisitions, and ρ is linked to shape of the grain kernel. Figure 1 presents isolated rice, wheat and barley seeds. The images are binarized and statistics on the features of isolated seeds are computed. Table 1 presents these statistics which will be helpful for the parameter selection in the algorithm. The major axis length (respectively minor axis length) is the scalar specifying the length in pixels of the major axis (respectively minor axis) of the ellipse that has the same normalized second central moments as the grain kernel.

Figure 2 presents the mean maps of the isolated seed shapes. A green ellipse is plotted on it to illustrate the interest of using an ellipse to estimate the statistics in Table 1. The mean shape maps are obtained by performing a partial procrustes superimposition of the isolated seed shapes. The grain kernel shape and their three symmetries (two axial and one central) are optimally translated and rotated before performing the superimposition.

Comparisons with ground truth orientation were performed in order to determine the parameter ρ. The samples for these experiments are presented in Table 1. Geometric statistics on the grain kernels presented in Figure 1.

Rice seeds, Figure 1 sample contains one hundred grain kernels. Table 2 presents the statistics on the absolute value of the angle error for these hundred grain kernels. For all the four studied grain kernels, on average the error is of 2 to 3 degrees, same for the standard deviation. The median is around 2 degrees. There are very few bad orientation estimations, two of them are presented in Figure 8. One can note that the main estimation errors are due to overlapping or broken grain kernels.

Conclusion

The present work justifies the interest of structure tensors to evaluate orientations of grain kernels in images. Comparison with ground truth results proves that even if no arrangement is provided by the user, the grain kernel orientation can be retrieved in the samples. The tests of the present study explain how to find the best parameters for the structure tensor. The amount of postsmoothing is linked to the minor axis of the ellipse which best approximate the grain kernels. The accuracy of the orientation estimation demonstrates the efficiency of the structure tensor as a preprocessing tool for segmentation methods. Implementation on industrial segmentation algorithms is done in [16].

Acknowledgement

The images were provided by Alpha MOS company which develops systems for quality control of food products. The acquisition system is described in [START_REF] Dubosclard | Automatic method for visual grading of seed food products[END_REF]8]. In the barley case, the bad orientation is due to an overlapping.

Fig. 1 .

 1 Fig. 1. Seeds used to extract the geometric statistics presented in Table 1. Top line, rice and wheat seeds. Bottom line, two types of barley seeds.

Fig. 2 .

 2 Fig. 2. Mean maps of the shapes: from left to right, rice, wheat, first and second type of barley grain kernels, from top to bottom, with or without a green ellipse.

Fig. 3 .

 3 Fig. 3. Influence of ρ on the mean (left) and the median (right) of the absolute value of the angle errors.

Figures 4 ,

 4 Figures 4, 5, 6 and 7. The best parameter is the one which provide the lowest mean or median of the absolute value of the angle error. The curves in Figure 3 illustrate the influence of ρ on the mean (left) and the median (right) of the absolute value of the angle errors. The rice, wheat, first type and second type of barley results are represented in respectively red star, green plus, blue circle, and black cross. In each case, the best postsmoothing parameter ρ is around one-fourth of the mean pixel length of the minor axis.

Figures 4 , 5 ,

 45 Figures 4, 5, 6 and 7 present a visual comparison between the ground truth orientations (left and plotted in red) and estimated orientations (right and plotted in green) for respectively the rice, wheat and the two barley grain kernels. Each

Fig. 4 .

 4 Fig. 4. Comparison of ground truth orientations (left) and estimated orientations (right) for the rice grain kernels.

Fig. 5 .

 5 Fig. 5. Comparison of ground truth orientations (left) and estimated orientations (right) for the wheat grain kernels.

Fig. 6 .

 6 Fig. 6. Comparison of ground truth orientations (left) and estimated orientations (right) for the barley grain kernels, first type.

Fig. 7 .

 7 Fig. 7. Comparison of ground truth orientations (left) and estimated orientations (right) for the barley grain kernels, second type.

Fig. 8 .

 8 Fig.8. Zooms on two bad orientation estimations. The two images of rice grain kernels at left are extracted from Figure4. The two images of barley grain kernels at right are extracted from Figure6. In the rice case, the bad orientation is due to a broken kernel. In the barley case, the bad orientation is due to an overlapping.

Table 2 .

 2 Statistics on the absolute value of the angle error.

	Seeds	ρ Minimum Maximum Mean Median Standard deviation
	Rice	3	0.150	14.833 2.806 2.266	2.629
	Wheat	5	0.023	15.025 2.628 1.798	2.494
	Barley 1 5	0.043	24.473 2.266 1.572	2.722
	Barley 2 6	0.006	9.274 2.464 1.825	2.202


Franc ¸ois Bertaux

Szymon Stoma

Dirk Drasdo

Gregory Batt
email: gregory.batt@inria.fr

Modeling Dynamics of Cell-to-Cell Variability in TRAIL-Induced Apoptosis Explains Fractional Killing and Predicts Reversible Resistance

Isogenic cells sensing identical external signals can take markedly different decisions. Such decisions often correlate with pre-existing cell-to-cell differences in protein levels. When not neglected in signal transduction models, these differences are accounted for in a static manner, by assuming randomly distributed initial protein levels. However, this approach ignores the a priori non-trivial interplay between signal transduction and the source of this cell-to-cell variability: temporal fluctuations of protein levels in individual cells, driven by noisy synthesis and degradation. Thus, modeling protein fluctuations, rather than their consequences on the initial population heterogeneity, would set the quantitative analysis of signal transduction on firmer grounds. Adopting this dynamical view on cell-to-cell differences amounts to recast extrinsic variability into intrinsic noise. Here, we propose a generic approach to merge, in a systematic and principled manner, signal transduction models with stochastic protein turnover models. When applied to an established kinetic model of TRAILinduced apoptosis, our approach markedly increased model prediction capabilities. One obtains a mechanistic explanation of yet-unexplained observations on fractional killing and non-trivial robust predictions of the temporal evolution of cell resistance to TRAIL in HeLa cells. Our results provide an alternative explanation to survival via induction of survival pathways since no TRAIL-induced regulations are needed and suggest that short-lived anti-apoptotic protein Mcl1 exhibit large and rare fluctuations. More generally, our results highlight the importance of accounting for stochastic protein turnover to quantitatively understand signal transduction over extended durations, and imply that fluctuations of short-lived proteins deserve particular attention.

Introduction

TNF-Related Apoptosis Inducing-Ligand (TRAIL) is a promising therapeutic agent against cancer because it induces apoptosis specifically in tumor cells [START_REF] Wiley | Identification and characterization of a new member of the TNF family that induces apoptosis[END_REF][START_REF] Walczak | Tumoricidal activity of tumor necrosis factor-related apoptosis-inducing ligand in vivo[END_REF][START_REF] Johnstone | The TRAIL apoptotic pathway in cancer onset, progression and therapy[END_REF]. This motivated dozens of clinical trials based on TRAIL-related therapies. However, efficiency was usually limited [START_REF] Dimberg | On the TRAIL to successful cancer therapy? Predicting and counteracting resistance against TRAIL-based therapeutics[END_REF]. Most of the molecular events leading from TRAIL exposure to cell death are known [START_REF] Johnstone | The TRAIL apoptotic pathway in cancer onset, progression and therapy[END_REF]. After TRAIL binding to death receptors, initiator caspases are activated, which in turn promote effector caspases activation either directly or via a mitochondrial pathway (Fig. S1). In most cells, Mitochondrial Outer Membrane Permeabilization (MOMP) is required to efficiently activate effector caspases. Several kinetic models have been proposed to describe a part or all of those biochemical reactions [START_REF] Albeck | Modeling a snap-action, variable-delay switch controlling extrinsic cell death[END_REF][START_REF] Chen | Modeling of the role of a Baxactivation switch in the mitochondrial apoptosis decision[END_REF][START_REF] Fussenegger | A mathematical model of caspase function in apoptosis[END_REF][START_REF] Legewie | Mathematical modeling identifies inhibitors of apoptosis as mediators of positive feedback and bistability[END_REF][START_REF] Neumann | Dynamics within the CD95 death-inducing signaling complex decide life and death of cells[END_REF][START_REF] Rehm | Systems analysis of effector caspase activation and its control by X-linked inhibitor of apoptosis protein[END_REF][START_REF] Schliemann | Heterogeneity reduces sensitivity of cell death for TNF-stimuli[END_REF][START_REF] Bentele | Mathematical modeling reveals threshold mechanism in CD95-induced apoptosis[END_REF].

Not all cells of an isogenic population die after TRAIL treatment, even at saturating ligand doses. This fractional killing property is widely shared among cell lines and is critical for therapeutical applications [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF]. In addition, surviving cells were shown to be transiently resistant to a second TRAIL treatment. This reversible resistance property was observed in various cell lines and could also have important implications for therapy [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF][START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF]. Fractional killing is generally thought to result from cross talks between the apoptosis pathway and survival pathways [START_REF] Falschlehner | TRAIL signalling: decisions between life and death[END_REF]. Indeed, several studies reported that TRAIL increases the production of anti-apoptotic proteins, via the activation of survival pathways [START_REF] Chaudhary | Death receptor 5, a new member of the TNFR family, and DR4 induce FADDdependent apoptosis and activate the NF-kappaB pathway[END_REF][START_REF] Son | TRAIL-activated stress kinases suppress apoptosis through transcriptional upregulation of MCL-1[END_REF][START_REF] Sun | TRAIL-induced caspase/p38 activation is responsible for the increased catalytic and invasive activities of Akt[END_REF]. While fractional killing illustrates cell-to-cell variability in the decision between life and death, variability is also observed among cells that die: they commit to death after a highly variable delay from one another [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Rehm | Dynamics of outer mitochondrial membrane permeabilization during apoptosis[END_REF]. This variability cannot be explained by differences in TRAIL-induced gene regulation: it is also observed when cells are co-treated with cycloheximide (CHX), an efficient inhibitor of protein synthesis [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF]. Rather, it was proposed to originate mostly from pre-existing differences in the levels of proteins composing the apoptosis pathway. Indeed, recently divided sister cells die almost synchronously [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Rehm | Dynamics of outer mitochondrial membrane permeabilization during apoptosis[END_REF], as expected if protein content is equally shared between daughters and if noise in signaling reactions play a marginal role. This explanation was supported by modeling: when taken as initial conditions of a deterministic, kinetic model describing signaling reactions, differences in protein levels are sufficient to explain observed variability in death times [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF].

Another observation indicated that this cell-to-cell variability in protein levels is not frozen but result from a dynamical equilibrium driven by fluctuations in individual cells: death synchrony between sister cells was weaker as the duration between division and treatment increased [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Rehm | Dynamics of outer mitochondrial membrane permeabilization during apoptosis[END_REF]. To quantitatively assess this effect, the previously mentioned modeling approach, where only the consequences of protein fluctuations (cell-to-cell variability at a fixed time) are accounted, is inadequate. Instead, protein fluctuations themselves should be modeled. The need to account for protein fluctuations is even more stringent when considering observations after treatment with TRAIL alone. In that case, protein synthesis is not blocked and thus can impact the decision between life and death. Even if TRAIL does not change protein synthesis (via induction of survival pathways), significant differences with TRAIL and CHX treatments are expected, as the constitutively noisy protein synthesis will interact with signaling reactions. Thus, before complexifying the model to account for eventual regulations via survival pathways induction, it is critical to assess how much can be explained when protein synthesis is not altered by TRAIL signaling. Here, we investigate this question by enabling a kinetic model of TRAIL-induced apoptosis with stochastic protein turnover for all proteins, following a generic and principled approach. To our knowledge, this is the first attempt to systematically include gene expression noise in a signal transduction model. It enriches the model with a fundamental property as the dynamics of cell-to-cell variability is represented, allowing disentangling the effects of constitutive protein fluctuations, signaling protein-protein reactions and potentially induced changes in protein synthesis.

Results

Extended vision of signal transduction pathways

Protein synthesis and degradation are subjected to noise, resulting in fluctuations of protein concentrations in individual cells and in cell-to-cell variability at the population level [START_REF] Mcadams | Stochastic mechanisms in gene expression[END_REF][START_REF] Raser | Noise in gene expression: origins, consequences, and control[END_REF]. Such variability could have consequences on signal transduction: aside of conventional epigenetic differences [START_REF] Rando | Timescales of genetic and epigenetic inheritance[END_REF], unequal access to ligand molecules or simply noise in signaling reactions, it often contributes importantly to heterogeneous behavior within an isogenic population [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF]. One approach to account for those differences is to incorporate protein level variability as random initial conditions of an ODE model describing the signaling reactions (''extrinsic noise approach'') [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Gaudet | Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis[END_REF]. However, variability is imposed at time zero and then behavior is deterministic: it is thus not appropriate to study transduction on long time scales, during which protein levels dynamically fluctuate [START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF]. A more natural manner to account for protein level variability is to represent their stochastic synthesis and degradation (''intrinsic noise approach''). Although several studies did account for cell-to-cell differences in protein levels in an extrinsic, static manner via random initial conditions [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Aldridge | Lyapunov exponents and phase diagrams reveal multi-factorial control over TRAILinduced apoptosis[END_REF][START_REF] Stoma | STL-based analysis of TRAIL-induced apoptosis challenges the notion of type I/type II cell line classification[END_REF], and many models of signal transduction considered the effect of noise in protein-protein reactions [START_REF] Shibata | Noisy signal amplification in ultrasensitive signal transduction[END_REF][START_REF] Lapidus | Intrinsic noise, dissipation cost, and robustness of cellular networks: the underlying energy landscape of MAPK signal transduction[END_REF] or in the expression of signal transduction target genes [START_REF] Tay | Single-cell NF-kappaB dynamics reveal digital activation and analogue information processing[END_REF][START_REF] Lipniacki | Single TNFalpha trimers mediating NF-kappaB activation: Stochastic robustness of NF-kappaB signaling[END_REF][START_REF] Lipniacki | Stochastic regulation in early immune response[END_REF], no kinetic model of signal transduction pathways considering systematically noise in protein synthesis and degradation has been developed so far. Here by systematically we mean for all the proteins acting in the pathway.

We propose a modeling approach to account for gene expression noise within kinetic models of signal transduction pathways. Following Singh et al. [START_REF] Singh | Dynamics of protein noise can distinguish between alternate sources of gene-expression variability[END_REF], we model protein turnover with stochastic processes describing mRNA level fluctuations, and deterministic processes for protein translation and degradation (random telegraph model, see also [START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Suter | Mammalian genes are transcribed with widely different bursting kinetics[END_REF][START_REF] Dar | Transcriptional burst frequency and burst size are equally modulated across the human genome[END_REF][START_REF] Vin ˜uelas | Quantifying the contribution of chromatin dynamics to stochastic gene expression reveals long, locus-dependent periods between transcriptional bursts[END_REF][START_REF] Paszek | Modeling stochasticity in gene regulation: characterization in the terms of the underlying distribution function[END_REF][START_REF] Munsky | Listening to the noise: random fluctuations reveal gene network parameters[END_REF][START_REF] Sherman | A Computational Framework for Analyzing Stochasticity in Gene Expression[END_REF]). These processes are integrated into a kinetic model of protein-protein reactions (Fig. 1). While the rates of such stochastic protein turnover models (Fig. 2A) are rarely directly measurable, their value can be constrained by using experimentally measurable data and analytical results (Fig. S2). Recently, significant progress has been made on both experimental and theoretical sides to enable this inference approach [START_REF] Singh | Dynamics of protein noise can distinguish between alternate sources of gene-expression variability[END_REF][START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Suter | Mammalian genes are transcribed with widely different bursting kinetics[END_REF][START_REF] Dar | Transcriptional burst frequency and burst size are equally modulated across the human genome[END_REF][START_REF] Vin ˜uelas | Quantifying the contribution of chromatin dynamics to stochastic gene expression reveals long, locus-dependent periods between transcriptional bursts[END_REF][START_REF] Paszek | Modeling stochasticity in gene regulation: characterization in the terms of the underlying distribution function[END_REF][START_REF] Munsky | Listening to the noise: random fluctuations reveal gene network parameters[END_REF][START_REF] Sherman | A Computational Framework for Analyzing Stochasticity in Gene Expression[END_REF]. Importantly, we found that for typical protein and mRNA half-lives, a large set of promoter rate combinations lead to similar fluctuations at the protein level, as characterized by protein level variability (coefficient of variation) and mixing time (half-autocorrelation time). Interestingly, the obtained mixing time is around 40 hours (Fig. 2B), in the middle of the range of experimentally estimated values for twenty endogenous proteins in human cells [START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF]. Thus, standard stochastic protein turnover models (Fig. 2C) can provide a good approximation of protein fluctuations for most proteins. Only short-lived proteins necessitate particular attention (Fig. S3). This finding is a cornerstone of our approach.

Modeling stochastic protein turnover in TRAIL-induced apoptosis

We applied this approach to TRAIL-induced apoptosis, using EARM kinetic model [START_REF] Albeck | Modeling a snap-action, variable-delay switch controlling extrinsic cell death[END_REF][START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] to describe protein-protein reactions taking place between TRAIL exposure and cell death commitment. We equipped all native proteins with a default model of stochastic protein turnover, with the exception of proteins with fast turnover, here Flip and Mcl1 [START_REF] Nijhawan | Elimination of Mcl-1 is required for the initiation of apoptosis following ultraviolet irradiation[END_REF][START_REF] Poukkula | Rapid turnover of c-FLIPshort is determined by its unique C-terminal tail[END_REF][START_REF] Laussmann | Proteasome Inhibition Can Impair Caspase-8 Activation upon Submaximal Stimulation of Apoptotic Tumor Necrosis Factor-related Apoptosis Inducing Ligand (TRAIL) Signaling[END_REF]. Details of model construction are given in Text S1 (see dedicated section, Fig. S4 and Tables S1-3). Importantly, all parameters have been constrained based on experimental data and analytical results, with the exception of four parameters (''ON'' and ''OFF'' promoter switching rates, mRNA and protein half-lives for Flip and Mcl1). Note that because of their similar protein and mRNA half-life, we use the same couple of promoter switching rates for both proteins in our exploration of parameter space. This drastically limited the number of introduced degrees of freedom and made it possible to systematically explore realistic ranges for remaining parameters. To study the influence of stochastic protein turnover on fractional killing and reversible resistance, we sought to confront our model with existing quantitative data about TRAIL-induced apoptosis in HeLa cells. Those experiments, described in detail later, can be classified into two groups based on the type of information they contain: 1) quantification of the variability in cell fate, 2) characterization of the transient memory

Author Summary

TRAIL induces apoptosis selectively in cancer cells and is currently tested in clinics. Having a mechanistic understanding of TRAIL resistance could help to limit its apparition. Several observations suggested that protein level fluctuations play an important role in TRAIL resistance and its acquisition. However, quantitative, systems-level approaches to investigate their role in cellular decisionmaking processes are lacking. We propose a generic and principled approach to extend signal transduction models with protein fluctuation models for all proteins in the pathway. The key aspect is to use standard protein fluctuation models for long-lived proteins. We show that its application to TRAIL-induced apoptosis provide a quantitative, mechanistic explanation to previously published but yet unexplained critical observations. in cell state. While previous approaches using ODE models with distributions for initial protein levels (capturing a static description of cell-to-cell variability) [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Gaudet | Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis[END_REF] are potentially able to reproduce the first type of data, a dynamic view on cell-to-cell variability as proposed in our model is needed to account for both types of data. We adopted the following strategy: first, search for models able to reproduce observations on cell fate variability; and second test whether valid models can robustly predict observed behaviors where transient memory matters.

Stochastic protein turnover models predict transient memory in cell sensitivity to TRAIL and CHX

Using live-cell microscopy, Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] investigated the fate of hundreds of cells after exposure to TRAIL and CHX (10 ng/mL and 2.5 mg/mL, Fig. 3A). All cells undergo MOMP with a highly variable delay (from 2 to 8 hours, Fig. 3B). To study cell fate inheritance, the authors also recorded 20 hours before treatment to identify sister cells (Fig. 4A). They were found to have highly correlated MOMP times (correlation coefficient close to 1 for recently divided cells, about 0.5 for older sisters -Fig. 4B, black curve). Here, the MOMP time distribution provides a quantification of the cell fate variability, while MOMP time correlations between sister cells also give information on the transient memory in cell state. Within our framework, in-silico reproduction of those experiments is straightforward (Figs. 3D and4C), enabling us to investigate possible origins of transient cell fate inheritance. We first asked if the observed cell fate variability could be reproduced. In the model, it is only determined by protein levels at treatment time (behavior is deterministic as synthesis is assumed to be fully blocked by CHX and noise in signaling reactions is neglected), and differences between sister cells are only caused by protein synthesis noise occurring between division and treatment (in agreement with the fact that recently divided sisters died almost synchronously, we assumed an equal repartition of protein content at division). We found that excellent agreement with observed MOMP time variability can be obtained (Fig. 3E). Further analysis revealed that such agreement requires Flip and Mcl1 protein half-life to be short and to fall within a narrow range (between 0.3 and 0.6 hours, Fig. S5-C). This model prediction is consistent with previous measurements in HeLa cells (30 and 40 minutes for Flip short isoform and Mcl1 respectively, [START_REF] Nijhawan | Elimination of Mcl-1 is required for the initiation of apoptosis following ultraviolet irradiation[END_REF][START_REF] Poukkula | Rapid turnover of c-FLIPshort is determined by its unique C-terminal tail[END_REF]). In contrast, Flip and Mcl1 mRNA halflife and promoter switching rates are not strongly constrained, probably because their influence on cell fate is limited by the rapid protein level decrease caused by synthesis blockade. We then asked whether our extended model also capture transient cell fate inheritance (Fig. 4B). It is the case: fitted models accurately predict the MOMP time correlation between sister cells (Fig. 4D -black curve, Fig. S5-C). Of note, assuming standard promoter switching rates for Flip and Mcl1 (but accounting for their short mRNA and protein half-life -this parameterization will later be referred as the ''non-fitted'' model) already provides a good agreement for both MOMP time distribution and MOMP time correlation between sister cells (Fig. S6). This non-trivial result shows that the speed at which the sensitivity to TRAIL and CHX fluctuates in single cells is well captured and thus suggests that our generic approach permits to describe fluctuations of protein levels with sufficient accuracy.

Large, rare fluctuations of short-lived anti-apoptotic proteins explain fractional killing and predict transient cell fate inheritance

Spencer and colleagues repeated this experiment but treated cells with TRAIL alone (250 ng/mL). In this condition, an important fraction of cells died fast (MOMP in ,2 hours) but 40% were still alive after 8 hours (Fig. 3C), illustrating the fractional killing property. Also, cell fate inheritance between sister cells was markedly changed: only young sister cells that underwent MOMP rapidly were importantly correlated (Fig. 4B, grey curve). We asked whether the observed cell fate variability, including fractional killing, could be reproduced in-silico. Within our modeling assumptions, absence of co-treatment with CHX makes a fundamental difference: as synthesis continues, the effect of gene expression noise during TRAIL-induced apoptosis could be investigated, and comparison with the TRAIL and CHX condition is insightful. Strikingly, we found that quantitative agreement for both MOMP time distribution and surviving fraction could be obtained (Fig. 3F). Robustness analysis showed that rates of the Flip and Mcl1 stochastic protein turnover model, and particularly promoter switching rates, are in this case strongly constrained. Interestingly, MOMP time distribution and surviving fraction constrain those values differently (Fig. 5A-B), resulting in an narrow ranges for their values: agreement for both observations together is obtained only when promoter switching rates are both low (Fig. 5C). Such low switching rates lead to large, rare fluctuations of protein levels (Fig. 5D). Those atypical fluctuations phenotypes are expected to leave a signature at the population level: the shape of the protein level distribution would be bimodal rather than resembling a lognormal distribution (Fig. S7). This property is thus a model prediction. Those fluctuations are likely to impact how the fate of sister cells diverge with time. Thus, we asked whether the model could also account for the observed fast loss of cell fate inheritance. Remarkably, the fitted models accurately and robustly predict MOMP time correlations between sister cells (Figs. 4D and5E). As mentioned earlier, the same couple of promoter switching rates was used for Flip and Mcl1 during exploration, but further analysis showed that assuming low promoter switching rates for Mcl1 alone was sufficient to obtain quantitative agreement for MOMP distributions surviving fractions, and that sister cells MOMP time correlations were still correctly predicted (Fig. S8). Thus, comparison with transient cell fate inheritance data supports that large, rare fluctuations of Mcl1 could be responsible for the observed cell fate variability.

Accounting for stochastic protein turnover predicts reversible resistance

Recently, reversible resistance was observed among various cell lines [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF]. Cell populations were submitted to two consecutive TRAIL treatments. The duration between treatments was varied from 1 day to 1 week (Fig. 6A). One-day survivors were significantly more resistant than the initial population, but such resistance was significantly decreased or even lost in one-week survivors. Thus, cells surviving a first TRAIL treatment are transiently resistant. Remarkably, in-silico reproduction of those (Fig. 6B) showed that our model predicts the presence of reversible resistance (Fig. 6C): one-day survivors exhibit a dose-dependent increase of resistance to a second TRAIL treatment, which disappears after 3 to 5 days. This is surprising since our model does not include induced regulation mediated by survival pathways. Moreover, the presence of reversible resistance is a robust property of the model as it is also obtained when assuming standard promoter switching rates for Mcl1 and Flip (''non-fitted'' model, Fig. S6). However, agreement with related experimental data [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF][START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF] is only qualitative. While we cannot exclude that model parameterizations allowing a quantitative agreement exist, it might be needed to include additional mechanisms such as survival pathways induction to explain the observed sustained resistance gain after one week when treating cells with a high TRAIL dose [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF].

Molecular determinants of fractional killing and reversible resistance

What are the mechanisms behind cell escape to TRAILinduced apoptosis, either on the short-term (fractional killing) or the long-term (reversible resistance)? Using the fact that in-silico, all protein, mRNA levels and gene activity states can be monitored in single cells, we investigated those questions at the molecular level. To study the influence of pre-existing differences on cell fate, we compared at the time of stimulation the sub-population of 'future survivors' with the whole population (Fig. 7A-B). Future survivors strongly stood out by their Mcl1 protein level and gene activity state (Fig. 7B). Flip also appeared to play an important role in determining cell decision, and smaller but significant effect was also seen for Bid, Bax, Bcl2 and XIAP. Although it is a good predictor of cell fate, initial Mcl1 gene activity status does not completely determine survival: neither all Mcl1 ''ON'' cells survived nor all Mcl1 ''OFF'' cells died. Thus, pre-existing differences in protein levels and promoter activities are major determinant of cell fate but stochastic events in gene expression occurring during signal transduction also play a role. While timing of death for cells treated with TRAIL and CHX appeared to be multi-factorial [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF], our results suggest that cell survival is predominantly determined by Mcl1 (Fig. 7B). This important role of Mcl1 is robustly predicted. Indeed, it also holds for the ''nonfitted'' model, which assume standard promoter switching rates for all proteins, including Mcl1 and Flip (Fig. S6).

To investigate the determinants of reversible resistance, we tracked the temporal evolution of protein levels in surviving cells (Fig. 7 A&C). The protein level composition of one day survivors contrasts with the protein content observed in future survivors: almost all protein levels differ importantly from the naı ¨ve population composition, while that was the case only for Mcl1 and Flip in future survivors. This is expected as all proteins are partly activated during signal transduction, leading to a higher degradation (active forms have a shorter half-life). Therefore, the distinction between the causes of cell survival and the consequences of partial apoptosis induction cannot be easily resolved by the sole observation of protein levels in survivors. When signaling stops, recovery of protein levels is expected to follow exponential kinetics governed by the turnover rate (Fig. 7C, inset -see the death receptor (R), pro-caspase 8, Bar and Bid). Deviation from such kinetics indicates either the persistence of signaling reactions that continue to consume proteins (it is the case for Apaf, procaspase 9 and XIAP, as further analysis confirmed) or is a consequence of important selection. Indeed, while Mcl1 and Flip should recover normal levels in a few hours in absence of selection because of their high turnover rate, Mcl1 levels (but not Flip levels) are still strongly higher than in naive cells one day after TRAIL treatment, consistently with previous observations on the relative selection strength that operated on them. Together, those results indicate that recovery phenotypes in surviving cells result from a complex interplay of three distinct effects: selection during apoptosis, transcriptional noise and protein turnover as a driving force tending to reset protein levels to their initial, pre-stimulus distribution, and long-term residual signaling activity. This explains why it is difficult to understand the recovery process and justifies the use of modeling to disentangle the various contributions.

Discussion

Rehm et al. and Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Rehm | Dynamics of outer mitochondrial membrane permeabilization during apoptosis[END_REF] made two insightful observations about TRAIL-induced apoptosis. First, recently born sister cells died almost synchronously when treated with TRAIL and a protein synthesis inhibitor, while in contrast, unrelated cells died after highly variable durations. This demonstrated that TRAIL signaling is mostly deterministic when protein synthesis is blocked and that the timing of death is determined by the cell internal state at the time of treatment. Second, they observed that such synchrony in sister cells death is gradually lost as the time between division and treatment increases. This showed that the cell 'TRAIL sensitivity state' (the part of cell internal state involved in death timing determination) naturally fluctuates over a dozen of hours. In addition, the modeling results in ref. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] highly suggested that such state is mainly composed by the various levels of the proteins acting in the extrinsic apoptosis pathway. In parallel, important progress on the characterization of the stochasticity in gene expression has been made: the two-state transcriptional bursting model was shown to permit high accuracy and several approaches to infer its parameters were proposed, enabling the quantitative modeling of protein fluctuations in single cells [START_REF] Singh | Dynamics of protein noise can distinguish between alternate sources of gene-expression variability[END_REF][START_REF] Raj | Stochastic mRNA synthesis in mammalian cells[END_REF][START_REF] Suter | Mammalian genes are transcribed with widely different bursting kinetics[END_REF][START_REF] Dar | Transcriptional burst frequency and burst size are equally modulated across the human genome[END_REF][START_REF] Vin ˜uelas | Quantifying the contribution of chromatin dynamics to stochastic gene expression reveals long, locus-dependent periods between transcriptional bursts[END_REF][START_REF] Paszek | Modeling stochasticity in gene regulation: characterization in the terms of the underlying distribution function[END_REF][START_REF] Munsky | Listening to the noise: random fluctuations reveal gene network parameters[END_REF].

Modeling protein fluctuations in TRAIL-induced apoptosis

In this study, we merged those two approaches by integrating such stochastic models of gene expression within an existing kinetic model of TRAIL-induced apoptosis [START_REF] Albeck | Modeling a snap-action, variable-delay switch controlling extrinsic cell death[END_REF] in a systematic and principled manner. Doing so provides advantages compared to previous approaches to account for cell-to-cell variability in protein levels [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF][START_REF] Gaudet | Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis[END_REF]. First, variability is not considered as an ''input'' parameter but arises naturally from stochastic fluctuations. The dynamics of this variability is thus intrinsically represented within the system, allowing investigating the effects [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF][START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF] and model predictions (our study). Comparison of resistance gain experimental measurements between the two TRAIL doses (500 and 50 ng/mL for [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF] and [START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF] respectively) should be done with care, as the measurement and quantification method differed. doi:10.1371/journal.pcbi.1003893.g006 of transient memory in protein levels. Second, the influence of protein synthesis noise during TRAIL-induced apoptosis could also be investigated. Importantly, we followed a parsimonious parameterization strategy, motivated by the fact that fluctuations of long-lived proteins are rather insensitive to the precise kinetics of transcriptional bursting, enabling us to equip most proteins (long-lived proteins) with reasonably accurate fluctuation models even in absence of gene expression data for each and every promoter. The sister cells experiment for which cells were treated with TRAIL and CHX provided ideal data to validate our modeling approach: in that case, behavior is mostly deterministic as soon as treatment starts and only fluctuations occurring before treatment are responsible for death time variability and decorrelation between sister cells. Moreover, gene regulation via survival pathways induction is ineffective as protein synthesis is blocked. Because our model was able to quantitatively reproduce the MOMP time distribution and then accurately predicted sister cells correlation, our modeling approach appears as a promising tool to investigate the effect of protein fluctuations on signal transduction, despite the limitations inherent to its simplicity (for Figure 7. Molecular determinants of cell fate and resistance to repeated TRAIL treatments. (A) Cartoon illustrating that the determinants of cell fate and resistance can be studied by analyzing the over-representation of protein levels in 'future survivors' (cells that will still be alive after treatment) at the time of treatment, and in surviving cells at day X, respectively. (B) Cell fate determinants analysis: over-representation (compared to initial population) of protein level (blue) and promoter activity (red) at the time of treatment in 'future survivors'. Asterisks mark differences that passed a 5% significance test. (C) Resistance determinants analysis: over-representation of protein levels in surviving cells at day X. Inset illustrates the recovery kinetics expected from protein turnover only (i.e. in absence of significant selection effect or residual signaling activity). Therefore, deviation from such kinetics indicates the presence of a selection effect or residual signaling activity. doi:10.1371/journal.pcbi.1003893.g007 example, stochastic gene expression events were assumed to be independent between proteins, neglecting the fact that levels of different proteins can be partially correlated [START_REF] Gaudet | Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis[END_REF][START_REF] Sigal | Variability and memory of protein levels in human cells[END_REF], possibly because of common transcription factors or coordinated chromatin-state transitions). Of note, good agreement was readily obtained when assuming for Flip and Mcl1 standard promoter switching rates, but short protein and mRNA half-lives -in agreement with available knowledge (''non-fitted'' model, Fig. S6). Finally, the transposition of fluctuation timescales from individual proteins into 'TRAIL sensitivity states' is not trivial: while in our model, stable proteins levels are mixed in about 40 hours, cells were switching between 'fast dying' and 'slow dying' phenotypes more rapidly (about 10-15 hours). As combinatorial effects are at play, mechanistic models of protein-proteins reactions are needed to link protein-level timescales with more high-level phenotypic transitions [START_REF] Gupta | Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells[END_REF].

Questioning the role of survival pathways

Several studies reported that TRAIL can induce survival pathways [START_REF] Chaudhary | Death receptor 5, a new member of the TNFR family, and DR4 induce FADDdependent apoptosis and activate the NF-kappaB pathway[END_REF][START_REF] Son | TRAIL-activated stress kinases suppress apoptosis through transcriptional upregulation of MCL-1[END_REF][START_REF] Sun | TRAIL-induced caspase/p38 activation is responsible for the increased catalytic and invasive activities of Akt[END_REF]. How such induced changes affect signal transduction and eventually stop apoptotic signaling remains unclear. On the other hand, the contribution of constitutive protein synthesis noise, which is responsible for pre-existing differences between cells, has not been evaluated. Although it does not exclude the existence of other mechanisms, an important result of our study is that fractional killing can be obtained without assuming any TRAIL-induced regulation. Alternatively, we find that because of its fast turnover, constitutive expression of the Mcl1 protein has the potential to rescue cells from TRAIL apoptotic signaling. In this context, solely accounting for protein fluctuations within the TRAIL apoptosis pathway predicts the fractional killing property (Figs. 3 andS6). While our results challenge current opinion on the role of survival pathways in TRAIL-induced apoptosis, they are consistent with observations made on wild type HeLa cells that neither blocking NF-kB response nor inhibiting the Akt pathway do significantly change the surviving cell fraction after TRAIL treatment [START_REF] Braeuer | Constitutively activated nuclear factor-kappaB, but not induced NF-kappaB, leads to TRAIL resistance by up-regulation of X-linked inhibitor of apoptosis protein in human cancer cells[END_REF][START_REF] Lalaoui | TRAIL-R4 promotes tumor growth and resistance to apoptosis in cervical carcinoma HeLa cells through AKT[END_REF]. The pivotal role for Mcl1 in TRAIL-induced apoptosis predicted by our model is consistent with the recent finding that Mcl1 silencing by shRNA in HeLa cells completely sensitize cells to TRAIL [START_REF] Lemke | Selective CDK9 inhibition overcomes TRAIL resistance by concomitant suppression of cFlip and Mcl-1[END_REF]. While moderate fluctuations of Mcl1 levels were sufficient to obtain fractional killing, a quantitative agreement with the Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] single-cell data (MOMP time distribution and surviving fraction) required large and rare Mcl1 fluctuations, caused by rare switches between long periods of gene activity or inactivity. Interestingly, in that case, the observed rapid loss of MOMP time correlation between sister cells quantitatively emerged from model simulations. Flip is often mentioned as a key factor in cell resistance to TRAIL [START_REF] Shirley | Targeting c-FLIP in cancer[END_REF], but in our model Flip has less impact on cell survival than Mcl1. Consistently, Lemke et al. [START_REF] Lemke | Selective CDK9 inhibition overcomes TRAIL resistance by concomitant suppression of cFlip and Mcl-1[END_REF] silencing experiments demonstrated a dominant role for Mcl1 and a synergy with Flip. However, our model might underestimate the role of Flip: the representation of DISC-related events in EARM is simple and thus does not account for recent biological findings, including the stoichiometry between its components [START_REF] Dickens | A death effector domain chain DISC model reveals a crucial role for caspase-8 chain assembly in mediating apoptotic cell death[END_REF][START_REF] Schleich | Stoichiometry of the CD95 death-inducing signaling complex: experimental and modeling evidence for a death effector domain chain model[END_REF]. Improving how DISC assembly is modeled might thus be needed to elucidate the precise role of Flip in fractional killing and reversible resistance, especially for cell lines that express higher Flip amounts than HeLa.

Origins of reversible resistance: joint effect of selection and stochastic protein turnover

A second significant result reported here is that our model predicts the phenomenon of reversible resistance, showing that constitutively noisy protein synthesis, protein-protein interactions and protein degradation are by themselves sufficient to explain a dose-dependent, significant increase of resistance in recent survivors and its gradual loss within 3-5 days. This result is consistent with the observation that Nf-kB blockade does not change resistance acquisition after TRAIL treatment [START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF] (in MCF10A cells; HeLa cells have not been tested). In-silico analysis at the molecular level revealed that reversible resistance as predicted by the model was shaped by a complex interplay between 1) selection based on protein levels and transcriptional activity, 2) protein turnover and 3) residual signaling activity. As opposed to the death process, which involves a sharp and complete activation of effector caspases, our results suggest that recovery in cells that did not commit to death is a slow and complex process. While one should not conclude from our results that parallel activation of survival pathways by TRAIL plays no role in reversible resistance, our results show that the sole contribution of protein level fluctuations occurring within the extrinsic apoptosis pathway can partly lead to reversible resistance. Thus, protein fluctuations should be accounted for to gain quantitative insights into reversible resistance.

Accounting for gene expression noise appears necessary to investigate signal transduction

While here we focused on TRAIL-induced apoptosis, our modeling approach is generic and can be applied to other signal transduction pathways. Our results showed that even in absence of induced gene regulation, gene expression noise interacts with signaling dynamics on a non-trivial manner. Thus, even in contexts where the influence of induced gene-regulation is indisputable, its sound quantification probably requires to investigate first the role of constitutive gene expression noise. Only then models could be enriched parsimoniously with wellcharacterized regulatory links until all observations are successfully explained. Significant advances to allow such detailed characterization of gene regulation occurred recently [START_REF] Tay | Single-cell NF-kappaB dynamics reveal digital activation and analogue information processing[END_REF][START_REF] Molina | Stimulusinduced modulation of transcriptional bursting in a single mammalian gene[END_REF][START_REF] Neuert | Systematic Identification of Signal-Activated Stochastic Gene Regulation[END_REF]. Following such approaches could significantly extend the reach of models of signal transduction towards accurate, single-cell level description of populations submitted to varying signaling contexts over multiple cell generations. We interpret the four first reactions as stochastic reactions and the translation and protein degradation reactions as deterministic reactions. Thus, the model state is given by a Boolean G~G on ~1{G off , an integer m and a continuous variable P. The statistical properties of this system can be investigated analytically by writing the corresponding chemical master equation, which describes the temporal evolution of the state joint probability distribution. While the obtained equations cannot be solved analytically, it can be used to derive moments of the underlying distribution by applying generating functions tech-niques, as was done by Paszek [START_REF] Paszek | Modeling stochasticity in gene regulation: characterization in the terms of the underlying distribution function[END_REF]. We next give the analytical expressions of several moments of interest:

Methods

Modeling stochastic protein turnover

E G ½ ~k on k on zk off ~kon c g c g ~kon zk off À Á CV G ½ 2 ~1{E G ½ E G ½ E m ½ ~E G ½ k sm c m , CV m ½ 2 ~cm c m zc g CV G ½ 2 z 1 E m ½ E P ½ ~E m ½ k sp c p , CV P ½ 2 ~cp c m c p zc m zc g À Á c m zc p À Á c m zc g À Á c p zc g À Á CV G ½ 2 z c p c m zc p 1 E m ½
While those expressions characterize the steady-state distribution of the model, and thus can be compared to snapshot measurements of a cell population, they do not provide information about how fast single cells move with time within those distributions. Such information can be obtained by computing the autocorrelation function. Its expression and derivation are given in Supplementary Results (Text S1).

Modeling TRAIL-induced apoptosis and protein fluctuations

We used the EARM model described in Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] to represent protein-protein reactions. No kinetic rates were changed, except the rate of pC6 cleavage by C3, set to zero to remove the feedback loop. Degradation rates of non-native forms were modified according to biological knowledge (Table S3). The effect of those changes is discussed in SI Appendix. Rates of the standard stochastic protein turnover models and rates of Flip&Mcl1 models were determined following algorithm to incorporate constraints based on biological knowledge (Fig. S4). Their values are given in Table S1-2. Notably, for comparability and consistency with the findings of Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF], we kept their mean protein values, even if for a few proteins new data were available, as for Mcl1 for example [START_REF] Lindner | Systems Analysis of BCL2 Protein Family Interactions Establishes a Model to Predict Responses to Chemotherapy[END_REF], provided that the order of magnitude is the same.

In-silico ''sister cells'' and ''repeated TRAIL'' experiments

The simulation procedure is detailed in SI Appendix. Briefly, mother cell states (promoter activity, mRNA and protein levels) were chosen by Monte-Carlo sampling and the two sister cells were constructed by duplication of the mother cell state. Promoter activity and mRNA fluctuations were simulated using an implementation of the Gillespie algorithm. For the ODEs governing evolution of all protein levels, the Semi-Implicit Extrapolation method was used. MOMP was considered to have occurred when half of mitochondrial Smac has been released. For all simulations, at least 10 4 pairs of sister cells were simulated. Quantitative comparison of simulation results with related experimental data (Fig. 5) was achieved by computing a measure of model-data agreement described in Supplementary Methods (Text S1). For repeated TRAIL simulations, a naı ¨ve population of 10 4 cells was obtained as in the sister cells experiment. Cells divided after a gaussian cell cycle duration (mean 27 hours, standard deviation 3 hours). Death was detected via cPARP levels as in Gaudet et al. [START_REF] Gaudet | Exploring the contextual sensitivity of factors that determine cell-to-cell variability in receptor-mediated apoptosis[END_REF].

Supporting Information

Figure S1 Simplified view of the TRAIL-induced apoptosis pathway. Only the type of each protein-protein interaction (activation or inhibition) is represented. A diagram displaying the detail of protein complexes and activated forms can be find in [START_REF] Albeck | Modeling a snap-action, variable-delay switch controlling extrinsic cell death[END_REF]. (PDF)

Figure S2 Numerical simulation and analytical characterization of stochastic protein turnover models. A stochastic protein turnover model is defined by six rates: the promoter activity switching rates, the mRNA production and degradation rates, the protein per mRNA synthesis rate and the protein degradation rate. Numerical simulations can be used to simulate temporal fluctuations in single cells. When a population of cells is simulated, the cell-to-cell variability can be studied. After some time, cell-to-cell variability reaches a steady state. Analytical calculations on the stochastic protein turnover model provide expressions characterizing the steady-state variability (moments of the steady-state distributions), but also fluctuations (autocorrelation functions). Complete expressions are given in Text S1. (PDF)

Figure S3 Fluctuations of protein levels caused by transcriptional bursting are smoothed out for long-lived proteins. Comparison of protein level coefficient of variation (A) and half-autocorrelation time (C) as a function of transcriptional bursting rates for two situations: a short-lived protein and mRNA (half-lives of 2 and 1 hours, resp.) and a long-lived protein and mRNA (27 and 9 hours, resp.). Other rates of the stochastic protein turnover model are chosen such that mean protein and mRNA level are the same (1000 and 17).

Figure 1 .

 1 Figure 1. Accounting for stochastic protein turnover in signal transduction pathways. Scheme of the modeling approach. Protein-protein interactions mediating signal transduction (signal transduction layer) are modeled by ordinary differential equations. In parallel, promoter activity changes, mRNA production and degradation (gene expression layer) are seen as stochastic events and generate fluctuations in mRNA levels. This impacts the synthesis rates of the corresponding proteins. Together with protein degradation, it generates fluctuations in protein levels (here shown in absence of transduction). Only a fragment of the extrinsic apoptosis pathway is shown. Deterministic/stochastic interpretation of chemical reactions is represented with black/red arrows respectively. doi:10.1371/journal.pcbi.1003893.g001

Figure 2 .

 2 Figure 2. Standard stochastic protein turnover models. (A) Schematic description of the reactions constituting the stochastic protein turnover model. Gene activity switches, mRNA production and degradation (red arrows) are stochastic reactions. Protein synthesis and degradation reactions (black arrows) are deterministic. (B) For typical mRNA and protein half-lives, promoter switching rates have a limited influence on the protein level half-autocorrelation time: using an analytical derivation of the protein level autocorrelation function (see Supplementary Results, Text S1), the protein half-autocorrelation time is plotted against mean promoter switching times. (C) Behavior of a standard stochastic protein turnover model. Promoter switching rates respect typical ranges observed by Suter et al. [36] and lead to a protein level coefficient of variation (CV) of 0.25. See Fig. S3 for more details. Upper plots show three representative single-cell time courses of protein and mRNA levels. Histogram at the bottom displays the corresponding distribution of protein level obtained when simulating a large number of cells for a long duration, corresponding to a snapshot of the cell-to-cell variability expected in a population. doi:10.1371/journal.pcbi.1003893.g002

Figure 3 .

 3 Figure 3. Cell fate variability in TRAIL-induced apoptosis. (A-C) Cell fate variability experiments performed in [13]. (A) HeLa cell populations were treated with either 10 ng/mL of TRAIL and 2.5 mg/mL of cycloheximide (CHX) or 250 ng/mL of TRAIL alone. Cells were tracked during 8 hours by live-cell microscopy and MOMP time was detected via mitochondrial release of a fluorescent reporter. (B-C) Histograms of MOMP times and surviving fractions observed for treatment with (B) TRAIL and CHX or (C) TRAIL alone. (D-F) In-silico reproduction of those experiments with our ''fitted'' model (i.e. the parameterization in the explored parameter space region giving the best agreement for cell fate variability data, described in Tables S1-3). (D) Simulations (see Supplementary Methods in Text S1 for details). (E-F) Results for the (E) TRAIL and CHX or (F) TRAIL alone treatments. For the latter case, representative model trajectories are given in Fig. S9. doi:10.1371/journal.pcbi.1003893.g003

Figure 4 .

 4 Figure 4. Transient cell fate inheritance in TRAIL-induced apoptosis. (A-B) Experiments measuring correlation of MOMP times between sister cells performed in [13]. (A) HeLa cells were recorded from 20 hours before treatment as in Fig. 3A. Sister cells were identified to permit comparison of their fate. (B) Quantification of cell fate inheritance was realized by computing the correlation between sister cells MOMP time as a function of the duration between division and MOMP (averaged between sisters). (C-D) In-silico reproduction of those experiments with the model of Fig. 3. (C) Description. (D) Quantification of cell fate inheritance was applied to simulation results as in (B). See Supplementary Methods in Text S1 for details. doi:10.1371/journal.pcbi.1003893.g004

Figure 5 .

 5 Figure 5. Model fitting to cell fate variability data predicts large, rare fluctuations of Flip/Mcl1 and transient cell fate inheritance. (A-C) Agreement between model prediction and experimental data for (A) death (i.e. MOMP) time distribution, (B) surviving fraction after 8 hours, and (C) both together, for treatment by TRAIL alone (250 ng/mL), as a function of Flip/Mcl1 promoter switching times (other parameters as in Table S1). (D) Representative protein level fluctuations of Mcl1 described by a stochastic protein turnover model allowing good agreement for both MOMP time distribution and surviving fraction. This model has been used for Figs. 3,4,6 and 7. (E) Model-data agreement for MOMP time correlation between sister cells. For (A), (B), (C) and (E), agreement quality increases from red to green. The quantification algorithm is detailed in Supplementary Methods (Text S1). doi:10.1371/journal.pcbi.1003893.g005

Figure 6 .

 6 Figure 6. Reversible resistance in repeated TRAIL treatments. (A) Schematic description of the 'repeated TRAIL' experiments performed in [14,15] to characterize reversible resistance in HeLa cells. (B) In-silico reproduction of these experiments with our model (details in Supplementary Methods, Text S1). (C) Resistance gains in surviving cells relative to naı ¨ve cells as a function of time between the two TRAIL treatments. Data are shown for experimental observations[START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF][START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF] and model predictions (our study). Comparison of resistance gain experimental measurements between the two TRAIL doses (500 and 50 ng/mL for[START_REF] Flusberg | Cells surviving fractional killing by TRAIL exhibit transient but sustainable resistance and inflammatory phenotypes[END_REF] and[START_REF] Flusberg | Modulating cell-to-cell variability and sensitivity to death ligands by co-drugging[END_REF] respectively) should be done with care, as the measurement and quantification method differed. doi:10.1371/journal.pcbi.1003893.g006

Denoting

  G on and G off the two states of the promoter, m the number of mRNAs and P the protein level, the stochastic protein turnover

1 (

 1 Translation and protein degradation).

  Figure S3 Fluctuations of protein levels caused by transcriptional bursting are smoothed out for long-lived proteins. Comparison of protein level coefficient of variation (A) and half-autocorrelation time (C) as a function of transcriptional bursting rates for two situations: a short-lived protein and mRNA (half-lives of 2 and 1 hours, resp.) and a long-lived protein and mRNA (27 and 9 hours, resp.). Other rates of the stochastic protein turnover model are chosen such that mean protein and mRNA level are the same (1000 and 17). Combinations of Ton and Toff values ranging from 0.1 to 5 hours and 0.1 to 10 hours respectively were tested. Ton and Toff are mean ON and OFF time of the gene. (B and D) Representation of the range of values obtained for all models tested in (A) and (C). (PDF) Figure S4 Building stochastic protein turnover models for TRAIL-induced apoptosis. Routine followed to choose rates of all 17 native proteins in the EARM kinetic model of TRAIL-induced apoptosis. Typical values from multi-genes studies in mammalian cells are used to constrain rate values. Specific attention is given to Flip and Mcl1 because they are known to be short-lived, and thus more prone to exhibit large variations. (PDF) Figure S5 Stochastic protein turnover models captures fluctuations of cell sensitivity to TRAIL and CHX. (A) Best found agreement between model and data for MOMP times distribution in the +CHX condition. Obtained for Flip and Mcl1 model rates such that protein/mRNA half-life and mean ON/ OFF promoter activity duration equaled 0.4/1 and 1.9/3.1 hours respectively. See Supplementary Methods (Text S1) for quantification of model data agreement. (B) Best found agreement between model and data for MOMP time correlation between sisters in the +CHX condition. Obtained for Flip and Mcl1 model rates such that protein/mRNA half-life and mean ON/OFF promoter activity duration equaled 0.3/1 and 0.35/24 hours respectively. (C) Influence of Flip and Mcl1 model rates on Model-Data agreement in the +CHX condition. For each parameter, we

PLOS Computational Biology | www.ploscompbiol.org

October 2014 | Volume 10 | Issue 10 | e1003893

Acknowledgments

We thank Laurence Calzone, Deborah Flusberg, Olivier Gandrillon, Olivier Micheau, Peter Sorger and Andrei Zinovyev.

This work was supported by the research grants Syne2arti ANR-10-COSINUS-007 and Iceberg ANR-IABI-3096 from the French National Research Agency. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

plot the model to data distance corresponding to the best model when all other three parameters are varied. (PDF)

Figure S6 The ''non-fitted'' model quantitatively predicts TRAIL+CHX single-cell data and lead to fractional killing and reversible resistance for TRAIL alone treatments. In the non-fitted model, Flip and Mcl1 promoter switching rates are standard (Ton = 0.1 hours and Toff = 2.6 hours) but the short-half life of their mRNA and protein is accounted for (2 hours and 0.5 hours respectively). We reproduce here for this model all the results presented in the main text for the ''fitted'' model. Quantitative agreement is obtained for TRAIL + CHX single cell data from Spencer et al. [START_REF] Spencer | Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF] (MOMP time distribution and sister cell MOMP time correlations). No quantitative agreement is obtained in the case of TRAIL alone treatments, but the existence of fractional killing and reversible resistance is nevertheless predicted. Note that because fewer cells were simulated compared to main text figures (50000 instead of 10 5 for sister cell experiments), sister correlation curves appears slightly noisier. (PDF)

Figure S7 Mcl1 and Flip fluctuations for standard or ''fitted'' promoter switching rates. For the top frame, promoter switching rates are standard (as in Fig. 2). Because mRNA and protein half-lives are short, protein level fluctuates more rapidly and the steady-state distribution is changed (it is wider and the mode is in 0) compared to the standard stochastic protein turnover model (Fig. 2). On the bottom frame, the steady-state distribution becomes bimodal because the promoter switching rates are low compared to mRNA and protein degradation. Figure S11 Resistance gain in one day survivors is robust regarding the presence/absence of feedback loop and the degradation of active forms. In-silico repeated TRAIL experiment (as in Fig. 6) was repeated for variants of the ''fitted'' model regarding presence/absence of the C3-.C6-.C8 feedback loop and the default active forms half-life. Resistance gain in one-day survivors is shown. Simulations were repeated 4 times with 10 4 cells, error bars indicate standard deviation of estimated resistance gain between replicates. (PDF)

Figure S12 Long-term population survival is not possible with the feedback loop and stable active forms. Evolution of alive cell number in populations treated in-silico as in Fig. 6, for model several variants regarding presence/absence of the C3-.C6-.C8 feedback loop and the default active forms halflife. (PDF)

Author Contributions

Conceived and designed the experiments: FB SS DD GB. Performed the experiments: FB SS. Analyzed the data: FB SS DD GB. Wrote the paper: FB GB. Developed numerical simulation code: FB SS.