Theoretical study of the impact of adaptation on cell-fate heterogeneity and fractional killing

Fractional killing illustrates the cell propensity to display a heterogeneous fate response over a wide range of stimuli. The interplay between the nonlinear and stochastic dynamics of biochemical networks plays a fundamental role in shaping this probabilistic response and in reconciling requirements for heterogeneity and controllability of cell-fate decisions. The stress-induced fate choice between life and death depends on an early adaptation response which may contribute to fractional killing by amplifying small differences between cells. To test this hypothesis, we consider a stochastic modeling framework suited for comprehensive sensitivity analysis of dose response curve through the computation of a fractionality index. Combining bifurcation analysis and Langevin simulation, we show that adaptation dynamics enhances noise-induced cell-fate heterogeneity by shifting from a saddle-node to a saddle-collision transition scenario. The generality of this result is further assessed by a computational analysis of a detailed regulatory network model of apoptosis initiation and by a theoretical analysis of stochastic bifurcation mechanisms. Overall, the present study identifies a cooperative interplay between stochastic, adaptation and decision intracellular processes that could promote cell-fate heterogeneity in many contexts.

In many adaptation and developmental contexts, isogenic cells make stochastic fate decisions to generate diversified cell types and subpopulations 1 . Cell-fate heterogeneity is indeed a key feature of microbial adaptation to adverse environments 2 , of the development and homeostasis of tissues and organs 3 or of tumor resistance to drug therapy 4 . The differential fate outcome of isogenic cells exposed to the same environmental stimuli involves the interplay of stochastic and deterministic mechanisms [5][6][7] , where regulatory mechanisms can determine both the statistics and dynamics of stochastic events and the effect of those stochastic events on molecular trajectories dictating cell fate choices. Several experimental studies have shown that cell-fate decisions are often preceded by a highly fluctuating intracellular dynamics. Pulsatile or oscillatory activities have been observed in signaling pathways operating during the stochastic choice of various differentiation or proliferation fates [8][9][10][11][12][13] . The profile characteristics of those dynamic signaling activities have been proposed to either direct decision outcomes or promote cell-fate heterogeneity 14,15 . Transient dynamics occurring at epigenetics, transcriptome-wide or multicellular levels [16][17][18] have also been proposed to regulate cell-fate heterogeneity and plasticity. All these examples support a key role of transient dynamics in orchestrating fate decisions from diverse signaling and stochastic cues.
An attractive case study is the stochastic fate decision between life and death, commonly termed fractional killing, for which the systematic measure of probabilistic dose-response curves coupled with single-cell analysis of stochastic and dynamical signatures are possible 19 . On this issue, several modelling studies have been devoted to identify which sources of fluctuations and which parts of the apoptotic network could contribute the most to the variability of decision time and outcomes [20][21][22][23] , while the impact of the transient dynamics has been seldomly addressed 24 . Yet, singe-cell analysis of the temporal trajectories of caspase 8 activity in response to TRAIL has revealed a signature of adaptation dynamics whose transient kinetics determines whether a cell survives or dies 25 . Caspase 8 is likely to be part of negative feedback regulation involving for instance the formation of inactive heterodimers of procaspase-8 26 . The importance of transient dynamics of apoptotic inducers has been also emphasized in the case of cisplatin drug exposure 24,27 . The proposed mechanism involves a competition between positive and negative regulation of caspase 8-dependent apoptosis, thereby defining an incoherent feedforward loop. More generally, environmental stressors are prone to upregulate both pro-survival and prodeath pathways [28][29][30][31] through negative feedback or incoherent feedforward loop motifs which ultimately lead to a dynamical adaptation response 32- 34 . These stochastic and deterministic features associated with fractional

Results
Stochastic modeling framework for probabilistic fate decisions. Fractional killing can be defined as the population-level property by which isogenic cells exposed to increasing doses of death-inducing stimuli will tend to display a fraction of surviving cells and dying cells, although with increasing probability of death. This stochastic decision process can be studied in a general theoretical framework that applies to cases of fate decisions other than survival and death. Probabilistic cell-fate response commonly involves the interplay between intracellular mechanisms of fate decision and intracellular sources of cell-to-cell variability. Without loss of generality, a possible framework to study such probabilistic fate response to a step stimulus consists in a Langevin differential equation description of the stochastic dynamics of a biochemical reaction network (see Table 1 for notations):  (2) and (9) s c Critical stimulus level without noise s c = s 50 (σ = 0) Fig. 2 x c (t, s c ) Critical trajectory Fig. 2 and Eqs. (7)(8) y(t) , y N Small deviations of x(t) from x c (t) Fig. 5 and Eqs. (7)(8)(9) �(t, t ′ ) Principal fundamental matrix Eq. (7) U(x), �, r K Effective potential, barrier height and Kramers rate Fig. 5 and Eq. (10) Scientific Reports | (2020) 10:17429 | https://doi.org/10.1038/s41598-020-74238-y www.nature.com/scientificreports/ where H(t) and δ(t) are respectively the Heaviside and the Dirac delta functions. In this model, the cell-to-cell variability of stimulus-induced response trajectories x l can originate either from noisy biochemical reactions involving stochastic processes ξ l j (t) or from heterogeneities in stimulus exposure/sensitivity s l or network parameters k l from cell to cell. Although limited or inaccurate to describe some stochastic behaviors of biochemical networks 41 , chemical Langevin equation approach is nevertheless convenient to study asymptotic cases of small noise, large size or separated timescales and, thus, to relate with noise-free dynamical properties 42 .
For the deterministic part of the equation, the fate-decision behavior (e.g., death) can be minimally implemented in the nonlinear dynamics of the biochemical network by the presence of a saddle-node bifurcation mechanism at a critical stimulus level s sn through which the (survival) steady state x st1 is destabilized toward the other (death) steady state x st2 , generally in an irreversible manner. Accordingly, near-identical cells exposed to the same stimulus may display divergent trajectories toward survival or death (Fig. 1a). The population dynamics of noisy or heterogeneous cells described by Eq. (1) can be statistically represented by a probability distribution P( x, t) , which typically follows a Fokker-Planck equation. Stimulus-induced fate decision relates with the establishment of a bimodal distribution such that one can define the decision probability P D (Fig. 1b): where t * is the typical measurement time and W s ( x st ) is the fate attractor basin. It is to emphasize that we consider the typical case of an irreversible fate decision associated with a low or null probability to revert from the state x st2 to x st1 , which is typically the case for death, proliferation or terminal differentiation fate outcomes. From the dose-response curve P D (s) , one can define a fractionality index (Fig. 1c) which quantifies the derivative of this curve around the 50% fate probability ( P D (s 50 ) = 0.5): Based on this simple sensitivity measure of the stochastic fate response, systematic analysis of how η varies with noise strength σ and network parameters k should provide key insights into the interplay of stochastic and nonlinear properties of networks in shaping probabilistic features of fate response.
Adaptation dynamics enforces a saddle-collision mechanism for decision making. To evaluate the impact of adaptation dynamics on the probabilistic properties of stochastic fate decisions, the biochemical reaction model used in Eq. (1) must implement adaptation and switching behaviors. For the ease of mathematical and graphical analysis, we consider a low-dimensional biochemical reaction network 43 , whose interactions between three coarse-grained variables implement a basic setting of a negative feedback-driven adaptation and a positive feedback-driven decision switch (Fig. 2a,b). Starting from a set of biochemical reactions associated with this architecture, a suitable factorization and normalization procedure (see "Methods" section) allows one to derive the following set of differential equation, (2) 2 implement a positive feedback that is strong enough to produce an irreversible transition to death fate. Second, the parameters k 1 = 0.9 and τ 1 = 0.1 < 1 < τ implement a significant and fast stimulusinduced response of x 1 , which gives rise to a marked overshoot dynamics through negative feedback with x 2 . Third, the synthesis rate parameters 1 − β , 1 − k 1 and k 2 = 0.056 satisfy the normalization condition that saddlenode bifurcation occurs for x 1 = x 2 = s = 1 whatever the values of the other model parameters ( k i , β , τ). Scientific Reports | (2020) 10:17429 | https://doi.org/10.1038/s41598-020-74238-y www.nature.com/scientificreports/ In this way, the parameters β ( ∈ [0 : 1] ) and τ can be systematically varied to modulate the overshoot characteristics of adaptation with limited impact on steady-state properties. Indeed, the steady state of x 1 depends on the stimulus according to: which satisfies x 1 (s = 1) = 1 for any β and k 1 values. Furthermore, stability analysis of this steady state establishes the following criteria for the existence of an overdamped overshoot response to the step stimulus s = 1: Accordingly, the adaptation parameters β ∈ [0 : 1] and τ > 1 control the amplitude and timescale of the overshoot response without changing the steady-state value for s = 1 (Fig. 2c). More details regarding the relation between the negative-feedback parameters and the adaptation behavior can be found in a previous modeling study 43 .
Before considering a source of variability, we first need to investigate the main effect of transient adaptation dynamics on the fate decision properties. In this case, probabilistic response and fractional killing do not occur, and the system response to a step stimulus is essentially determined by a threshold s c : a stress amplitude s greater (resp. lower) than s c leads to a fate decision toward death (resp., survival). For an adiabatically-slow increase of the stimulus, the system follows the steady-state branch x st1 (s) of low x 3 values before escaping from it for s > s sn through a saddle-node bifurcation. This is not necessarily the case for a step increase of the stimulus, for which transition to death can occur for a stimulus level s < s sn so that s sn − s c > 0 . The plot of s sn − s c as function of adaptation parameters in Fig. 2d shows that s c = s sn for weak enough or fast enough adaptation, while s c is below s sn for large enough values of both β and τ . These two qualitative regimes manifest, in fact, the existence of distinct instability mechanisms (Fig. 2e). For low values of β or τ and thus a small or no overshoot, the threshold property s = s c = s sn relates with a dynamical trajectory that is destabilized in the vicinity of a saddle-node instability (lower panel of Fig. 2e). In this scenario, x 3 species trigger the fate decision depending on the steadystate value of the adaptive species x 1 , which carries out the role of a bifurcation parameter. In sharp contrast, for high enough value of both β and τ and thus for a significant overshoot response of x 1 , the threshold property s = s c < s sn relates with a dynamical trajectory that collides with a saddle instability (upper panel of Fig. 2e).
In summary, while varying β and τ leads to gradual changes of the amplitude and the timescale of the overshoot adaptation profile, we could identify a threshold boundary in the {β, τ } space which separates between a saddle-node and a saddle-collision instability scenario. In the saddle-collision scenario, the threshold value s c becomes very sensitive to the transient characteristics of the overshoot profile, which suggests that fate decision may also become more sensitive to the sources of cell-to-cell variability that impact transient dynamics.
Adaptation dynamics promotes heterogeneous cell-fate decisions. Based on our comprehensive analysis of the deterministic decision dynamics in the coarse-grained model combining adaptation and bistability, we aim to investigate how adaptation influences cell-fate heterogeneity in a population of noisy cells. We therefore apply the general stochastic modeling framework to this biochemical network model (see "Methods" section) and perform a systematic analysis of the probabilistic properties in the {β, τ } parameter space and for several noise sources and levels ( Fig. 3).
To begin with, we consider the case of cell-to-cell variability arising from molecular noise solely ( s l = s and � k l = � k ) and focus on the two archetypical parameter sets depicted in Fig. 2 that are associated with weak/ fast adaptation ( β = 0.3 and τ = 3 ) and strong/slow adaptation ( β = 1 and τ = 10 ), respectively. Simulation of Langevin equations for 2000 trials shows that, for the same level of noise, strong/slow adaptation leads to probabilistic response associated with a much larger stimulus range and a much smaller derivative at P = 0.5 (Fig. 3a). These differences are quantified by the fractionality index η that is about four-fold larger for strong adaptation ( η ≈ 0.05 ) as compared to weak adaptation ( η ≈ 0.012 ). To illustrate how adaptation may amplify noise to generate more heterogeneous fate response, we plot the noisy single-cell trajectories in the two scenarios associated with the noise level and the same relative change of stimulus level. When adaptation is strong and slow enough, noisy trajectories remain within some neighborhood of the noise-free trajectory until diverging from it toward different fates when approaching the saddle fixed point, with a slight change of respective fate probability when the stimulus increases (Fig. 3b). This is in sharp contrast with the case of weak (or no) adaptation for which noisy trajectories reach first the neighborhood of a stable fixed point, before eventually escaping over the saddle fixed point toward the other fate when the stimulus slightly increases (Fig. 3c). The qualitative difference between these two stochastic decision scenario is confirmed by the distinct scaling laws η ∝ σ b , where b = 1 for strong/slow adaptation and b ≈ 1.2 for weak/fast adaptation (Fig. 3d). This body of evidences strongly suggest that adaptation dynamics promotes cell-fate heterogeneity, mostly by changing the underlying nonlinear mechanism of decision-making. This is confirmed by the plot η = f (β, τ ) (Fig. 3e), which unambiguously shows a qualitative increase of fractionality index η specifically in the parameter domain where the saddle-collision scenario occurs (above the white boundary).
Besides molecular noise, other sources of cell-to-cell variability have been tested, such as stimulus exposure or sensitivity s l (Fig. 3f) or initial conditions x l (t 0 ) (Fig. 3g). Again, a qualitative increase of η is observed in the parameter region associated with a saddle-collision scenario (above the white boundary), though the extent of such increase is much more important for the case of variable initial conditions. This is because variability of initial conditions impacts only transient dynamics, not steady state, while variability of s l impacts steady-state Adaptation dynamics contributes to fractional killing in an apoptosis model. The nonlinear nature of the adaptation-related amplification of noise effect suggests that this mechanism could be effective regardless the complexity of the network model. In other words, we expect to observe a similar noise amplification behavior in more detailed regulatory network model of stress-induced death fate decision as far as the adaptation dynamics leads to a collision to a saddle instability in the state-space of any dimensions. To check this conjecture, we need to replace the effective one-dimensional model of fate decision by a more realistic high-dimensional model of death fate decision. Fractional killing is commonly observed following many types of stress or death ligands, which may trigger death through different pathways 44,45 depending on the involved multi-protein signaling complexes, transcriptional factors and other signaling and metabolic cues (left of Fig. 4a). Among these many possibilities, we consider the canonical case where the stress signal and damage species mainly impact the intrinsic mitochondrial pathway of apoptosis through the control of the Bh3 member of Bcl-2 family 46 . An alternative possibility could have been to consider the case of TRAIL-induced apoptosis involving caspase 8-dependent activation of both extrinsic and mitochondrial pathways 20 .
The choice of Bh3-dependent mitochondrial apoptosis is motivated by a previous biochemical model of apoptosis initiation 24 , which exhibited several interesting features for our study. First, the model focuses on the Scientific Reports | (2020) 10:17429 | https://doi.org/10.1038/s41598-020-74238-y www.nature.com/scientificreports/ initial stage of intrinsic mitochondrial apoptosis, providing a simple picture of the decision-making process by leaving aside the further stages of effector caspase activation and related apoptosis events as well as the complex crosstalk with other programmed death pathways. Second, the topology and parameters of the model are determined in close relation with biological hypothesis and experimental data in a context of chemical stress response. Last, the topological and dynamical properties of the model are featured with a single positive feedback and a bistable behavior which are fully consistent with the minimal set of ingredients that is needed to implement our adaptation-dependent noise-amplification mechanism. In particular, prior knowledge about the bifurcation properties is very helpful to compare s c and s sn thresholds and make the connection between the low-dimensional and high-dimensional models. Therefore, the structure of the detailed model merely consists on the coupling between our above adaptation model (Eq. 4a, b) and the published model of mitochondrial apoptosis initiation 24 (right of Fig. 4a). Specifically, the adaptive species x 1 upregulates the synthesis of a pro-apoptotic BH3-only proteins (e.g., Bad, Bim, Bid), keeping in mind that intracellular stress-signaling pathways impacts the mitochondrial apoptosis pathway at various places 45,46 . Regarding the published apoptosis initiation model, the postranslational interactions between the pro-apoptotic Bh3 and Bax proteins and the anti-apoptotic Bcl-2 proteins implement a positive feedback mechanism. Pro-apoptotic signals are prone to increase the level of free Bh3 proteins with respect to the level of Bh3 proteins bound to Bcl-2. Free Bh3 proteins directly interact with inactive cytosolic Bax proteins, thereby inducing conformational change that leads to their activation and mitochondrial translocation. In turn, the activated mitochondria-localized form of Bax can also bind to Bcl-2, resulting in the release of additional free Bh3 proteins from Bh3-Bcl complexes. For a critical synthesis rate of Bh3 proteins, this positive feedback loop produces a bistable switching behavior via a saddle-node bifurcation from low to high levels of free mitochondrial Bax ( Bax m ) 24 . Then, high enough levels of Bax m would typically induce the release of cytochrome C and mitochondrial outer membrane permeabilization (MOMP) followed by the formation and activation of the apoptosome and the execution of apoptosis. The stochastic dynamics of this regulatory network coupling an adaptation module and an apoptosis module is simulated using again the Langevin formalism of Eq. (1) (see "Methods" section) with fixed σ . Death initiation event is assumed to occur when Bax m reaches the neighborhood of the high-level steady-state branch. For large number of simulation trials, death probability can be measured as a function of the stress stimulus level, s, for distinct adaptation profiles, here with β = 0 or 1 (Fig. 4b). Simulation results reveal that the presence of adaptation leads to a probabilistic response for a broader range of stimulus, which manifests itself by a lower value of the derivative of P death with respect to s/s 50 associated to a four-fold higher value of η . Such significant difference in noise-sensitivity η correlates to well distinct types of dynamical trajectories associated with survival-death fate decisions ( Fig. 4c-d). For β = 1 , the overshoot response of x 1 species leads to a transient increase of Bax m during which trajectories finally diverge from each other toward the survival or death attractor (Fig. 4c). This decision is made when approaching the unstable saddle equilibria along its stable manifold (right panel of Fig. 4c). For β = 0 , the level of Bax m reaches and fluctuates around a steady state of low values, before an eventual noiseinduced switch toward the death state by escaping over the saddle instability along its unstable manifold (Fig. 4d).
The results obtained with a detailed model of apoptosis are thus consistent with those obtained with the model Eq. (4) involving a minimal decision module (Fig. 4b similar to Figs. 3a and 4c-d similar to Fig. 3b-c). To obtain a similar behavior, it should be noted that (i) adaptation and decision timescales had to be adjusted to each other in the detailed model such that the decrease of x 1 during the overshoot profile occurs before Bax m reaches its upper-branch steady state, and that (ii) the molecular noise impacts more the death species than the adaptive species (compare molecular noise of Bax m and x 1 in Fig. 4c, d). A further step would thus be to check whether these two conditions are fullfiled in the detailed modeling of both the specific signaling pathways producing adaptation at the level of damage-repair pathways 47 , stress-response patwhays 24 or death-ligand pathways 25,26 , and of the specific death-regulatory pathways that are triggered by these diverse death-inducing stimuli. In these various cases, adaptation and fate decision processes are prone to be implemented by slightly different regulatory network topologies which may modulate the timescale and stochastic characteristics of the dynamical response and influence the extent of the adaptation-dependent fractional killing.

Theoretical description of stochastic decision properties.
We have shown that the sensitivity of cell-fate decision to molecular noise depends on the state-space paths taken to reach a saddle instability, along, either, its stable manifold or its unstable manifold (Fig. 5a). In order to get further insights into the stochastic nonlinear dynamics involved in this process, we develop a perturbation approach in the limit of small noise and small stimulus changes for which specific scaling laws η = aσ b have been obtained (Fig. 3d). Scaling analysis near instabilities is a common approach to characterize qualitative dynamical behaviors as function of noise, timescales and bifurcation parameters (see textbook 42 or some case study 48,49 ).
In the case of a saddle-collision scenario, perturbed trajectories evolve in the neighborhood of the deterministic trajectory x c (t, s c ) that connects the initial condition � x c (t = 0) = � x st1 (s = 0) and the saddle fixed point � x c (t → ∞) = � x sad (s c ) (Figs. 2d and 3d) and, thus, live on the stable manifold of this saddle W s ( x sad ) that separates the different fate attractors. Along this singular deterministic trajectory, some local Lyapunov stability exponents (i.e., time-dependent eigenvalues of the Jacobian matrix J( x c (t)) ) become positive such as to amplify transverse perturbations due to molecular noise or heterogeneous initial conditions. Mathematically speaking, linearization of Eq. (1a) about x c (t) defines a class of Langevin equations for the perturbed trajectories � y(t) = � x(t, s c + δs, σ ) − � x c (t) whose solution can be decomposed as � y(t) = δs � y δs (t) + σ � y σ (t) where �(t, t ′ ) is the principal fundamental matrix and b δs/σ are the normalized stimuli and noise perturbation vectors given by: To compute the fractionality index η , the key statement is that fate decisions (resp., probability) are determined by the deviation (resp., distribution) of trajectories onto the normal direction y N of the (N − 1)-dimensional stable manifold of the saddle. The mean and variance of the normal distribution P(y N , t) evolves in time until the decision time t * at which the distribution splits and many trajectories leave the neighborhood of x c (t) (Fig. 5b). Rewriting Eq. (2) as P D = R + P(y N , t * )dy N and decomposing dP D ds = ( dP D d�y N � )( d�y N � ds ) = (2πσ 2 �y 2 N,σ �) −1/2 (�y N,δs �) , we can derive the following expression for η: This formula shows an asymptotic scaling law η ∝ σ (Figs. 3d and 5d) while the prefactor depends in a sophisticated manner on the local stability properties (via ) and sensitivity properties (via b ) of the transient trajectory.
A similar derivation in the 1D case has been previously performed to show that this scaling law also depends on the speed of the trajectory toward the saddle instability 48 .
The sensitivity to noise is very different in the other scenario of noise-induced escape from a stable state ( x st1 ) over a saddle barrier ( x sad ), which is a very common behavior associated with the escape from metastability 49,50 . For this decision-making regime, an effective potential U, a potential barrier �(s) = U(� x sad (s)) − U(� x st1 (s)) and a Kramers escape rate r K (s) ∝ exp (2�(s)/σ 2 ) can be usually defined, even for multi-dimensional systems and multiplicative noise 51 (Fig. 5c). Given a fate-decision probability P D (t) ≈ 1 − exp (−r K t) , the fractionality index can be derived and approximated as : For the one-dimensional model of bistability used in Eq. (4c), the particular scaling relation s 50 ∂ s � ∼ σ 0.8 (as the threshold s 50 depends on σ ) leads to the scaling law η ∝ σ 1.2 obtained in Fig. 3d.
To conclude, these very distinct formulas for η highlight that the conversion of intracellular fluctuations into heterogeneous cellular fate response sharply differ depending on the transition scenario. The saddle-collision scenario is characterized with the amplification of small perturbations due to the local instability of trajectories when approaching the saddle state during the overshoot of decision variables (e.g., x 3 or Bax m ). In contrast, the more common scenario of a noise-induced escape from a metastable state does not display this amplification mechanism, while the transition rate r K is very sensitive to stimulus level due to the exponential-like dependency on the saddle barrier height.

Discussion
The present modeling study deciphers the role of adaptation dynamics in promoting cell-fate heterogeneity associated for instance with the fractional killing behavior. A common property of adaptation is the transient overshoot of some cellular variables above its steady state value, which can be implemented by diverse circuit topologies 32 and which is subjected to tradeoffs associated with homeostatic or sensory process 33,34,52,53 . In addition, we propose that this transient overshoot dynamics can also significantly impact fate-switching behaviors, so as to extend the stimulus range of fate heterogeneity and to allow for tunable fate probability. This adaptationdependent fate stochasticity relies on the manner how the overshoot of some intracellular species drive cell state in the neighborhoohd of a saddle instability, rather than a saddle-node instability, along a path where molecular noise are more prone to promote divergent decisions. This noise-amplification behavior illustrates how molecular noise and instability mechanisms can cooperate to shape cellular dynamics, like genetic timers 54 , boundary formation 55 or versatile sensory processing 56 .
The biological relevance of the proposed mechanism is most likely in a context of fractional killing for which the choice between life and death depends on adaptation processes. The timescales of those adaptation responses range from half an hour to few hours depending on stress type and regulatory mechanisms 47,57,58 which is of the range of magnitude of the initiator caspase rise time and death onset timing. Moreover, noise-induced fate heterogeneity is the most effective when fluctuating variables are those involved in the positive feedback that triggers death initiation. This requirement is consistent with modeling evidences that variability in diverse regulatory molecules can contribute in very different ways to variability in cell death outcomes 20 and that the main contributions seem to occur in the initial decision commitement phase, whether it is at the level of the fluctuations of short-lived antiapoptotic proteins 22 or the stochastic assembly of DISC/RIPoptosome platform 23 . The manner how cell fate is determined by the impact of these fluctuations at the level of concentration trajectories has been also investigated 24,25,27 . In relation to these studies, our study presents a broad and comprehensive view of this cooperative process and, thus, provides strategies, by monitoring transient characteristics, to either increase or reduce fractional killing.
The profile characteristics of adaptation dynamics, such as the ratio between its maximal and steady-state values, are highly sensitive to the temporal profile of the stimulus. Ramp increase of a stimulus or a preconditioning  59 , the osmotic stress response of yeast 60 or ethanol stress in Bacillus 61 . In case of stress-induced fate response, monitoring the stress stimulus profile would therefore be expected to modulate not only threshold stimulus level ( s 50 ), but also the degree of heterogeneity of the response ( η ). This provides a practical mean to test the role of adaptation for cell-fate heterogeneity, and to design dose delivery protocols of treatment to cope with fractional killing of cancer cells or microbial organisms. It is tempting to extrapolate the biological relevance of such adaptation-dependent mechanism beyond the scope of fractional killing and transient adaptation dynamics. The mechanism itself only requires a regulatory network featured with an upsteam overshoot response and a downstream switching response, which could be implemented by diverse network topologies and in diverse cell-fate contexts. For instance, overshoot dynamics can also occur in regulatory systems comprising incoherent feedforward loops 32 , but also in excitable or pulsatile systems combining negative and positive feedback loops. For the latter case, numerous signaling pathways, such as P53, Erk or NF-κ B, display a versatile pulsatile dynamics, which has been proposed to expand signal-processing capabilities and determine cell fate accordingly 14,15 . In relation to our result, the transient and stochastic characteristics of these signaling dynamics may also suggest a role for promoting cell-fate heterogeneity. This is supported by some experimental evidences that have mapped the cell-cell variability of the pulsing dynamics of Erk 10,11 , p53 62,63 , β-catenin 13 and NF-κB 12 with the heterogeneity of cell-fate outcomes. Data-driven and fine-grained modeling of specific dynamic signaling and fate-regulatory pathways 11,20,62,64 are definitively the step further to evaluate on a case-by-case basis to which extent transient adaptation or pulsing dynamics may contribute, fortuitously or functionally, to cell-fate heterogeneity.

Methods
Coarse-grained model. The set of regulatory reactions depicted in Fig. 3a consists in the following basal/ regulated synthesis terms and basal/regulated degradation terms: which can be translated into a system of differential equations using the law of mass action : To obtain the set of Eq. (4), we perform a nondimenzionalization procedure to reduce the number of parameters and to define effective parameters that control separately different features of the dynamics such as response timescales, transient nonlinear response and steady states. Accordingly, we have introduced dimensionless time t , concentration x i and stimulus s and defined rescaling variables ( X i,0 , S 0 ) and aggregate parameters ( τ i , β and k i ), as the following: These changes of variables and parameters simplify Eq. (12) into Eq. (4) where dimensionless time t is noted t again for simplicity. The chemical Langevin equation associated to Eq. (4) is also characterised with a rescaled noise σ = (�d 3 X 0 ) −1/2 where is the system size and X i,0 ≡ X 0 ∀i . Finally, reaction rates and stoichiometry matrices are given by: where k 1 = 0.9 , k 2 = 0.056 , k 3 = 1 , k 4 = 0.2 , τ 1 = 0.1 , while β and τ are varied.