Quantitative assessment of cell fate decision between autophagy and apoptosis

Autophagy and apoptosis are cellular processes that regulate cell survival and death, the former by eliminating dysfunctional components in the cell, the latter by programmed cell death. Stress signals can induce either process, and it is unclear how cells ‘assess’ cellular damage and make a ‘life’ or ‘death’ decision upon activating autophagy or apoptosis. A computational model of coupled apoptosis and autophagy is built here to analyze the underlying signaling and regulatory network dynamics. The model explains the experimentally observed differential deployment of autophagy and apoptosis in response to various stress signals. Autophagic response dominates at low-to-moderate stress; whereas the response shifts from autophagy (graded activation) to apoptosis (switch-like activation) with increasing stress intensity. The model reveals that cytoplasmic Ca2+ acts as a rheostat that fine-tunes autophagic and apoptotic responses. A G-protein signaling-mediated feedback loop maintains cytoplasmic Ca2+ level, which in turn governs autophagic response through an AMP-activated protein kinase (AMPK)-mediated feedforward loop. Ca2+/calmodulin-dependent kinase kinase β (CaMKKβ) emerges as a determinant of the competing roles of cytoplasmic Ca2+ in autophagy regulation. The study demonstrates that the proposed model can be advantageously used for interrogating cell regulation events and developing pharmacological strategies for modulating cell decisions.


ODE based Modeling
We describe the dynamics of our biochemical network using a system of ordinary differential equations (ODEs). For each molecular species x i taking part in the pathway there will be an equation of the form dxi dt = f i (x, c). Here f i describes the kinetics of the reactions that produce and consume x i , x are the molecular species taking part in these reactions while the vector of c gives the rate constants governing these reactions. We assume that the kinetic laws governing the individual reactions in the pathway are mass action law. Given that {x 1 , x 2 , . . . , x n } is the set of variables and {c 1 , c 2 , . . . , c m } the set of rate constants. Then each f i will be of the form f i = ri j=1 d j n ij g j , where r i is the number of reactions associated with species x i and d j = −1 (or d j = +1) if x i is a reactant (or product) of the j th reaction. Further, n ij ∈ Z denote the stoichiometric coefficients while g j are rational functions of the form g j = c α x a x b (mass action) with a, b ∈ {1, 2, . . . , n} and α ∈ {1, 2, . . . , m}. Each x i is real-valued function of t with t ∈ R + , where R + denotes the set of non-negative real numbers. We realistically assume that x i (t) takes values in the interval [L i , U i ] where L i and U i are non-negative rationals with L i < U i . Hence the state space of the system will be Thus V will be a bounded subset of R n + . To capture the cell-to-cell variability and uncertainties regarding the initial states we define for each variable The vector form of our system of ODEs is: dx/dt = F (x) with x = (x 1 , x 2 , . . . , x n ) and F (x(i)) = f i . We define the flow Φ : R + × V → V for arbitrary initial vectors v. Intuitively, Φ(t, v) is the state reached under the ODE dynamics if the system starts at v at time 0. The dynamics will be of interest only up to a maximal time point T . We define a trajectory starting from v ∈ V denoted σ v to be the (continuous) function The behavior of our dynamical system is the set of trajectories given by

Statistical Model Checking
We define a bounded linear time temporal logic (BLTL) as follows. An atomic proposition in our logic will be of the form (i, , u) with L i ≤ < u ≤ U i . Such a proposition will be interpreted as "the current concentration level of x i falls in the interval [ , u]". We fix a finite set of such atomic propositions AP = {A 1 , ......A k }. The formulas of BLTL are: (i) Every atomic proposition as well as the constants true, false are BLTL formulas. (ii) If ψ, ψ are BLTL formulas then ¬ψ and ψ ∨ ψ are BLTL formulas. (iii) If ψ, ψ are BLTL formulas and t ≤ T is a positive integer then ψU ≤t ψ and ψU t ψ are BLTL formulas. The derived propositional operators such as ∧, ⊃, ≡ and the temporal operators G ≤t , G t , F ≤t and F t are defined in the usual way. We will interpret the formulas of our logic at the finite set of time points T = {0, 1, . . . , T }. We assume T has been chosen such that it exceeds the last time for which experimental data is available. The semantics of the logic is defined in terms of the relation σ, t |= ϕ where σ is a trajectory in BEH and t ∈ T .
• ¬ and ∨ are interpreted in the usual way.
We define models(ψ) = {σ | σ, 0 |= ψ, σ ∈ BEH}. We shall make statements of the form P >r (ψ), where the intended meaning is that the probability that a trajectory in BEH belongs to models(ψ) exceeds r. To assign precise meaning to such a statement, we define a probability measure over sets of trajectories. We identify BEH with IN IT , and define the set M odels(ψ) ⊆ IN IT as: v ∈ M odels(ψ) iff there exists σ ∈ models(ψ) with σ(0) = v. To assign a probability to M odels(ψ), we construct a probability measure over the standard σ-algebra B generated by the open It is a standard fact that P extends in a unique way to the probability measure P : B → [0, 1] such that P (IN IT ) = 1 and P (∅) = 0. We define the formulas as P ≥r ψ and P ≤r ψ are probabilistic BLTL formulas provided r ∈ [0, 1) , r ∈ (0, 1] and ψ is a BLTL formula. We shall say that S, the system of ODEs meets the specification P ≥r ψ -and this is denoted S |= P ≥r ψ -iff P (M odels(ψ)) ≥ r, while S |= P ≤r ψ iff P (M odels(ψ)) ≤ r .
We formulate whether S |= P ≥r ψ, as a sequential hypothesis test between the null hypothesis H0 : p ≥ r + δ and the alternative hypothesis H1 : p ≤ r − δ, where p = P (M odels(ψ)). Here, δ signifies the indifference region supplied by the user. The strength of the test is decided by parameters α and β which bound the Type-I (false positive) and Type-II (false negative) errors respectively. Thus the verification is carried out approximately but with guaranteed confidence levels and error bounds. The test proceeds by generating a sequence of sample trajectories σ 1 , σ 2 , . . . by randomly sampling an initial state from IN IT . One assumes a corresponding sequence of Bernoulli random variables y 1 , y 2 . . ., where each y k is assigned the value 1 if σ k , 0 |= ψ. Otherwise y k is assigned the value 0. For each m ≥ 1, after drawing m samples, we compute a quantity q m as: Hypothesis H0 is accepted if q m ≥ A, and Hypothesis H1 is accepted if q m ≤ B. If neither is the case then another sample is drawn. The constants A and B are chosen such that it results in a test of strength (α, β). In practice, a good approximation is A = 1−β α and B = β 1−α .

Parameter Estimation using Statistical Model Checking
We encode experimental data as a BLTL formula ψ exp . Let O ⊆ {x 1 , x 2 , . . . , x n } be the set of variables for which experimental data is available and which has been fixed as training data to be used for parameter estimation. Assume . . , τ i Ti } are the time points at which the concentration level of x i has been measured and reported as is so chosen that it reflects the noisiness, the limited precision and the cell-population-based nature of the experimental data.
We then set ψ exp = i∈O ψ i exp . We encode qualitative dynamic constraints as a BLTL formula ψ qlty . We fix the probabilistic BLTL formula P ≥r (ψ exp ∧ ψ qlty ), where r will capture the confidence level with which we wish to assess the goodness of the fit of the current set of parameters to experimental data and qualitative trends. We also fix an indifference region δ and the strength of the test (α, β). We use the constants r = 0.9, δ = 0.05, α = 0.05 and β = 0.05.
Let θ = {c 1 , c 2 , . . . , c K } be the set of unknown rate constants whose values we wish to estimate. The outer loop of our parameter estimation procedure will run as follows. We shall assume for convenience that the search strategy uses a single set of parameter values (one for each unknown rate constant) in each round.
(i) Fix θ 0 , which assigns a value to each unknown rate constant. This represents the initial guess. Set = 0.
(ii) With θ as the current set of rate constant values, run the statistical model checking procedure to verify the individual conjuncts of ψ exp ∧ ψ qlty with the chosen strengths.
(iii) Based on the answers returned by these tests compute F (θ ), where F is the objective function.
(iv) Check if the value of the objective function is sufficiently high or has reached a predetermined bound.
(v) If yes, return θ as the estimated value.
(vi) Else fix a new set of rate constant values θ +1 as dictated by the search strategy. Increment to + 1 and return to step (ii).
The objective function is formed as follows. Let θ be an assignment of values to the unknown rate constants. Let J i exp (= T i ) be the number of conjuncts in ψ exp and J qlty the number of conjuncts in ψ qlty . Let J i,+ exp (θ) be the number of formulas of the form ψ t i (a conjunct in ψ i exp ) such that the statistical test for P ≥r (ψ t i ) accepts the null hypothesis (that is, P ≥r (ψ t i ) holds) with the strength ( α J , β), where J = i∈O J i exp . Similarly, let J + qlty (θ) be the number of conjuncts in ψ qlty of the form ψ ,qlty that pass the statistical test P ≥r (ψ ,qlty ) with the strength ( α J , β). Then F (θ) is computed via: Thus the goodness to fit of θ is measured by how well it agrees with the qualitative properties as well as the number of experimental data points with which there is acceptable agreement. To avoid over-training the model, we do not insist that every qualitative property and every data point must fit well with the dynamics predicted by θ.
The search strategy deployed in step (vi) above will use the values F (θ ) to traverse the space of candidate parameter vectors. We employ a global method called Stochastic Ranking Evolutionary Strategy (SRES) [2], since it is known to perform well in the context of pathway models [1].