A theoretical basis for cell deaths

Understanding deaths and life-death boundaries of cells is a fundamental challenge in biological sciences. In this study, we present a theoretical framework for investigating cell death. We conceptualize cell death as a controllability problem within dynamical systems, and compute the life-death boundary through the development of “stoichiometric rays”. This method utilizes enzyme activity as control parameters, exploiting the inherent property of enzymes to enhance reaction rates without shifting equilibrium states. This approach facilitates the efficient evaluation of the global controllability of models. We demonstrate the utility of our framework using its application to a toy metabolic model, where we delineate the life-death boundary. The formulation of cell death through mathematical principles provides a foundation for the theoretical study of cellular mortality. SIGNIFICANCE STATEMENT What is death? This fundamental question in biology lacks a clear theoretical framework despite numerous experimental studies. In this study, we present a new way to understand cell death by looking at how cells can or cannot control their states. We define a “dead state” as a state from which a cell cannot return to being alive. Our method, called “Stoichiometric Rays”, helps determine if a cell’s state is dead based on enzymatic reactions. By using this method, we can quantify the life-death boundary in metabolic models. The present framework provides a theoretical basis and a tool for understanding cell death.


I. INTRODUCTION
Death is one of the most fundamental phenomena, and comprehension of the life-death boundary is a pivotal issue in the study of biological systems.Cell death within simple unicellular model organisms such as Escherichia coli and yeast has been extensively studied [1][2][3][4][5][6][7][8][9][10][11].
What do we call "death" and from what features of a given cell do we judge the cell is dead?Despite extensive studies, such characterization of cell death is still under debate [11][12][13][14][15][16].Although experimental researches have explored this phenomenon in depth, theoretical frameworks for defining and characterizing cell death are conspicuously underdeveloped.
The current circumstances of microbial cell death research are in stark contrast to other topics in systems biology, where both experimental and theoretical approaches are integrated to explore cellular phenomena * yhimeoka@g.ecc.u-tokyo.ac.jp quantitatively and to elucidate the underlying biochemical design principles.Microbial adaptation, robustness of living systems, and network motifs are hallmarks of crosstalk between experimental and theoretical approaches [17][18][19][20][21].However, despite substantial experimental data, theoretical explorations of cell death remain underdeveloped.Theoretical studies on microbial death have primarily focused on estimating death rates and assessing their impacts on population dynamics, presupposing a predefined concept of "death" [5][6][7][8]10].Theories for delineating what is death and evaluating cellular viability are indispensable for extracting the quantitative nature of cell death, and accordingly, life.
The necessity for theories of cell death is beyond comprehending experimental observations; it is imperative for the advancement of biological sciences, particularly through the integration of burgeoning computational technologies.Recent developments in computational biology have facilitated the construction of whole-scale kinetic cell models [22] and the application of advanced machine-learning technologies, enabling large-scale simulations of cellular responses to environmental and genetic changes [23,24].This includes modeling scenarios for cell death as well as cell growth [22,25].The development of a robust mathematical framework to define and evaluate cell death within these models is essential to advance our understanding of cellular mortality.
In the present manuscript, we aimed to construct a theoretical foundation for cell death using cell models described by ordinary differential equations.We introduce the definition of the dead state based on the controllability of the cellular state to predefined living states.Next, we develop a tool for judging whether a given state is dead for enzymatic reaction models, the stoichiometric rays.The method leverages the inherent feature of enzymatic reactions in which the catalyst enhances the reaction rate without altering the equilibrium.The sto-ichiometric rays enables the efficient computation of the global controllability of the models.

II. A FRAMEWORK FOR CELL DEATHS
In the present manuscript, we formulate the death judgment problem as a controllability problem of mathematicals models of biochemical systems.Herein, we focus on well-mixed deterministic biochemical reaction systems.The dynamics of these systems are described by the following ordinary differential equations: where x : [0, T ] → R N ≥0 , u : [0, T ] → R P , and J : R P ×R N ≥0 → R R represent the concentration of the chemical species, control parameters, and reaction flux, respectively.Here we denote the number of chemical species, control parameters, and reactions as N, P and R, respectively, and S is N × R stoichiometric matrix.
Definition 1 (Trajectories).We gather all the solutions of Eq. (1) satisfying x(0) = x 0 and x(T ) = x 1 and denote the set by T (T, x 0 , x 1 ).In addition, the union of T (T, x 0 , x 1 ) for all T ≥ 0 is denoted as Note that T (T, x 0 , x 1 ) and T (x 0 , x 1 ) contain all the solutions obtained by varying the control parameter u(t).
Definition 2 (Controllable Set).The controllable set of given state x is defined as the set of states from which the state can be controlled to the target state x.It is given by We define the controllable set of a given set X as a union of the controllable set of points in X, that is, In the present framework, we assume that the representative living states X are given a priori as reference points of the living states (We provide several statemants on the choise of X later.).Using the controllable set and representative living states, we define the dead state as follows: Definition 3 (Dead State).A state z is called the dead state with respect to the representative living states The set of all dead states with respect to X form the dead set and is denoted as D(X), i.e., We term the boundary of the dead set, ∂D(X) as the Separating Alive and Non-life Zone (SANZ) hypersurface.
The SANZ hypersurface is derived from a mythical river in the Japanese Buddhist tradition, the Sanzu River, which separates the world of the living and the afterlife [26].
Note that the "death" is defined here based sorely upon the controllability of the state.The dead state can be a point attractor, limit cycle, infinitely long relaxation process, and etc.As Def. 3 states nothing on the features of the dead state except the controllability, the dead states can have, for instance, metabolic activity comparable to that of the representative living states.In such cases, the "dying state" would be an appropriate word to describe the state, while we use the term "dead state" for uniformity.
The dead state and dead set depend on the choice of the representative living states.The following parts of the section present useful statements for selecting the representative living states.
In the following, we denote x ⇀ y if x, y ∈ R N ≥0 satisfy x ∈ C(y).If x ⇀ y and y ⇀ x holds, we denote x ⇋ y.Proposition 1. ⇀ is preorder and ⇋ is equivalence relation on R N ≥0 .The proof is straightforward from the definition of the controllable set.
From the proposition, we obtain the following corollary.
holds.If x ⇋ y holds, we have Proof.For any z ∈ C(x), z ⇀ x holds.From the assumption, we also have x ⇀ y.Thus, the transitivity of ⇀ results in z ⇀ y.Therefore, C(x) ⊂ C(y) holds.By repeating the same argument for y ⇀ x, we have the second statement.
Thereby, if the states x, y ∈ R N ≥0 are mutually controllable to each other x ⇋ y, then the dead set is the same whether x or y is chosen as the representative living state.This implies that we can use the quotient set R N ≥0 /⇋ for studying the controllable set of the states instead of the original space R N ≥0 .In the following, we denote the equivalence class of x by This corollary states that the "larger" equivalence class in terms of 1 has the larger controllable set.This motivates us to introduce the terminal class of the representative living states X. Definition 4 (Terminal Class).For a given set X ⊂ R N ≥0 , we take the set of the equivalence classes [x] ∈ R N ≥0 /⇋ with nonempty intersection with X and denote it by X ⊂ R N ≥0 /⇋ i.e., We call the maximal element of X with respect to the terminal class of X.We denote the set of all terminal classes of X by Ω(X).
Note that an equivalence class [x * ] ∈ X being the terminal class means that there is no other equivalence class [x] exists.Ω(X) can be an empty, finite, or infinite set.However, if we can assume that every chain in X (a totally ordered subset of X ) has an upper bound with respect to , Zorn's lemma guarantees Ω(X) being nonempty.To satisfy this condition, it is sufficient that the model (Eq.( 1)) has a choice of the constant control parameter u with which the model never shows a divergent behavior.
The This allows us to reduce the representative living states X into a smaller set X * without changing the dead set.For X ⊂ R N ≥0 , suppose that Ω(X) is nonempty.Let us take a point from each terminal class in Ω(X) and denote the set of the chosen points by X * ⊂ R N ≥0 , that is, where the choise of the representative element from [x * ] is arbitrary.We have the following theorem for the dead set with respect to X and X * .
Theorem 1.The dead set with respect to X is identical to the dead set with respect to X * , i.e., As the controllable set of the equivalence class C([x]) and the controllable set of a point x ∈ [x] are identical (Cor.1), we have where X is what defined in Def. 4.
For each ) from Cor. 2. Therefore, we can drop C([y]) from the union in Eq. ( 15) without violating the equality.By repeating the dropping process, only the terminal classes remain and we obtain the following equality: Additionally, since The theorem states a useful property for choosing the representative living states to compute a dead set.Suppose that we have the representative living state X ⊂ R N ≥0 .Without the theorem, we must check the controllability of all points in X to judge the dead state.However, with this theorem, we can take discrete points from each terminal class in Ω(X) and check the controllability of the chosen points to judge the dead state.Note that the set X * can also be finite, which enables further efficient computations of the dead set.
An example of practical applications of this theorem is as follows.Suppose that the model has a point attractor a u with a fixed control parameter value u, and a u is a continuous function of u in U ⊂ R P .In such cases, taking the point attractor with a single parameter choice u 0 ∈ U is equivalent to taking the basin of attraction B(a u ) for all a u , as the representative living states, i.e., This is because a u1 ⇋ a u2 holds for any u 1 , u 2 ∈ U and also we have x ⇀ a u for any x ∈ B(a u ).

III. THE STOICHIOMETRIC RAYS: A TOOL FOR THE DEATH JUDGMENT
Thus far, we have proposed a possible definition of the dead state and presented several statements on the dead set.This definition raises the question of whether such global controllability can be computed.In this section, we show that the controllable set is efficiently computable for models with enzymatic activity and the external concentration of chemicals as the control parameter.Before moving to the main body, we briefly describe previous studies related to the controllability of chemical reaction models.
Most biochemical reactions can take place within reasonable timescales only by the catalytic aid of enzymes [27], and thereby, the control of intracellular states relies on the modulation of enzyme activities and concentrations (hereafter collectively referred to as activities).Thus, it is indispendable to develop a mathematical technique to compute the controllability of intracellular states by modulating enzyme activities for constructing a framework of cellular viability in terms of Def. 3.
Generic frameworks for the controllability of chemical reaction systems have been studied in the control theory field [28][29][30][31].However, there are critical limitations in previous studies when applied to biological systems; Saperstone et.al. [28] studied the controllability of linear reaction rate functions, with which we could not model many-body reactions.In Farkas et.al. [29] and Dochain et.al. [30], nonlinear models were dealt with, but only local controllability was discussed.In contrast, global controllability of a nonlinear model was studied by Drexler et.al. [31], although they allowed the kinetic rate constants to be negative-valued, to which no physical reality corresponds.
In the following, we present a simple method to compute the global controllability of enzymatic reaction models, the stoichiometric rays.The concept is hinted at the inherent feature of enzymatic reactions; enzymes enhance reactions without altering the equilibrium.We deal with the phase space to the positive orthant, R N >0 2 .
2 In typical cases, the method works for R N ≥0 , including the two Hereafter, we restrict our attention to the cases in which the reaction rate function in Eq, ( 1) is expressed by the following form We refer to f (x) and p(x) as the kinetic, and thermodynamic parts, respectively.f r (x) is a strictly positive function.n ± cr ≥ 0 represents the reaction degree of chemical c in the forward (+) or backward (−) reactions of the rth reaction.k r ≥ 0 corresponds to the Boltzmann factor of the reaction r.u ∈ R R ≥0 is the control parameter, representing the activities of the enzymes.In the following, we focus on the modulations of enzymatic activities where the forward and backward reactions cannot be modulated independently.However, when the forward and backward reactions are independently modulated, for example by controlling the external concentrations of chemicals, we can extend our method to deal with such situations by treating the forward and backward reactions as independent reactions.
Note that the popular (bio)chemical reaction kinetics have the form shown in Eq, ( 19)-( 21), such as the massaction kinetics, (generalized) Michaelis-Menten kinetics, ordered-and random many-body reaction kinetics, and the ping-pong kinetics [33].An important feature of the reaction kinetics of this form is that the direction of the reaction is set by the thermodynamic part as p(x) is related to the thermodynamic force of the reaction.On the other hand, the remaining part f (x) is purely kinetic, and it modulates only the absolute value of the reaction rate, but not the direction.For reversible Michaelis-Menten kinetics, v max ([S]−k[P ])/(K+[S]+[P ]), the thermodynamic part corresponds to its numerator, [S]−k [P ].Now, we allow the enzyme activities u to temporally vary in R R ≥0 and compute the controllable set of a given target state x * , C(x * ).The model equation is given by where ⊙ denotes the element-wise product.First, we divide the phase space R N >0 into subsets based on the directionalities of the reactions.Note that the kinetic part is strictly positive and u ∈ R R ≥0 holds.Consequently, the directionalities of the reactions are fully set only by the thermodynamic part p(x) 3 .Let σ be a examples presented in the latter half of the manuscript.However, it is known that special care is necessary for mathematical claims at the boundary (x i = 0 for one or more chemicals) [32].In the present manuscript, we avoid dealing with the boundary. 3For the directionality of the rth reaction at ur = 0, we define it by sgn Jr(u, x)| ur =0 = lim ur →+0 sgn Jr(u, x).
binary vector in {1, −1} R , and define the direction subset W σ as Next, we introduce the null-reaction manifold of reaction r given by Now the phase space R N >0 is compartmentalized by the null-reaction manifolds {M r } R−1 r=0 into the direction subsets W σ4 .Let path ξ(t) be the solution of Eq. ( 22) with a given u(t) reaching the target state x * at t = T from the source state x 0 , i.e., ξ(T ) = x * and ξ(0) = x 0 , meaning that x 0 ∈ C(x * ).Additionally, we assume that if ξ(t) intersects with any M i , the intersection is transversal 5 .As ξ(t) is the solution of Eq. ( 22), ξ(t) is given by Next, we fragment the path.Note that the union of all the closures of the direction subsets, W σ , covers R N >0 , and ξ(t) intersects transversally to the null-reaction manifolds, if any.Thus, we can divide the interval [0, T ] into L + 1 segments 25) is simplified as Recall that all the functions inside the integral are nonnegative.Thus, ũ(t) := u(t)⊙f (ξ(t))⊙|p(ξ(t))| satisfies ũ(t) ∈ R R ≥0 .This means that we can consider the ramped function ũ(t) as a control parameter.By introducing U i (t) := t τi ũ(s)ds, Eq, ( 26) is further simplified as Therefore, for the solution of Eq.( 22), ξ(t) with ξ(0) = x 0 and ξ(T ) = x * , the following holds; There exist a set of time-intervals {I i } L i=0 , reaction directionalities {σ (i) } L i=0 , and non-negative-valued functions {U i (t)} L i=0 such that path in the interval t ∈ I i is represented in the form of Eq, (27).
The converse of the above argument holds.Suppose that there is a continuous path ξ ⊂ R N >0 with x 0 and x * as its endpoints.Since ξ is a continuous path, it is continuously parameterizable by t ∈ [0, T ] so that ξ(0) = x 0 and ξ(T ) = x * hold.By following the above argument conversely, we can see that if there exist {I i } L i=0 , {σ (i) } L i=0 , and {U i (t) ≥ 0} L i=0 satisfying Eq, ( 27) and dU i (t)/dt ≥ 0, there exists a control u(t) given by and u(τ i ) = lim t→τi−0 u(t), where the division is performed element-wise.With this u(t), ξ(t) satisfies Eq, ( 22).Thus, x 0 is an element of the controllable set C(x * ).The summary of the above conditions leads to the definition of the single stoichiometric path.

Definition 5 (Single Stoichiometric Path
Let SP L (x * ) denotes the source points of all the single stoichiometric paths to x * with all possible choices of the sign sequence with a length less than or equal to L, and SP(x * ) := lim L→∞ SP L (x * ).We call SP(x * ) the stoichiometric paths while SP L (x * ) is termed the finite stoichiometric paths with length L. Note that the stoichiometric paths equals to the controllable set.
The stoichiometric paths is a useful equivalent of the controllable set, whereas the conditions 2 and 4 in the definition are laborious to check if a given path and U i (t) satisfies them.Thus, we introduce the single stoichiometric ray.It is easy to compute, and the collection of them gives exactly the stoichiometric paths if the thermodynamic part p(x) of a model equation is linear.
Suppose that there is a pair of points and satisfying Then, we can construct a ray connecting the two points, ξ(τ i ) and ξ(τ i+1 ), as where 0 ≤ s ≤ 1.
Here, U i (τ i+1 )s ≥ 0 and d(U i (τ i+1 )s)/ds ≥ 0 follow.Additionally, the points between ξ(τ i ) and ξ(τ i+1 ) are in the direction subset W σ (i) if the thermodynamic part is linear.Thus, ξ(t) and U i satisfy the conditions 2 to 4 in Def. 5 in t ∈ I i (the condition 3 holds only for the linear thermodynamic part).Since the timescales are arbitrarily chosen by modulating the magnitude of u(t), we can rewrite ξ(τ i ) and U i (τ i+1 ) as ξ i and U i , respectively, and obtain the following definition.

Definition 6 (Single Stoichiometric Ray). A set of points {ξ i } L+1
i=0 is called a single stoichiometric ray from x 0 to x * with signs {σ (i) } L i=0 if there exists a.

ξ
Let us define SR L (x * ) and SR(x * ), termed the finite stoichiometric rays with length L and the stoichiometric rays6 in the same manner as the stoichiometric path.Note that only the evaluations of L discrete points are required to determine the existence of a single stoichiometric ray.
Let us remark that SR L (x * ) = SP L (x * ) holds if the thermodynamic part p(x) is linear, and thus, the stoichiometric rays is equivalent to the controllable set.If the thermodynamic part is nonlinear, SP L (x * ) ⊆ SR L (x * ) holds in general.This is because the existence of control parameters in a given direction subset, U i (t), is evaluated based only on points on the null-reaction manifolds.
As shown in the illustrative example in Fig. 1A, the point x on M i is judged to be reachable to M j even though all possible paths are blocked by another direction subset (the region surrounded by M k in the figure).However, the stoichiometric rays correctly judge the reachability in cases illustrated in Fig. 1B.
Let us now show how the stoichiometric rays is computed using the reversible Brusselator model as an illustrative example.The model equation is given by d dt Note that the thermodynamic part vector p in the reversible Brusselator is a linear function of a and b whereas the reaction rate function vector f ⊙p is a nonlinear function.
The finite stoichiometric rays for the two choices of the target state x * are shown in Fig. 2

IV. DEATH OF A TOY METABOLIC MODEL A. Computing the dead set
Next, we apply the stoichiometric rays to a toy model of metabolism and compute its dead set.The model is an abstraction of the glycolytic pathway consisting of the metabolites X, Y, ATP, and ADP (a network schematic and reaction list are shown in Fig. 3A).
The rate equation is given by7 d dt , where x, y, z, and w represent the concentrations of metabolites X, Y, ATP, and ADP, respectively.In the second line of the above equation, f (x) = 1 is omitted.The thermodynamic part of the model is nonlinear.As the total concentration of ATP and ADP is a conserved quantity in the model, we replace w by z tot − z where z tot is the total concentration of the two (we set z tot to unity in the following).As shown in Fig. 3B, the model exhibits bistability with a given constant u.On one attractor (the endpoint of the red line), X molecules taken up from the environment are converted into Y and secreted into the external environment resulting in the net production of a single ATP molecule per single X molecule.In contrast at the other attractor (the endpoint of the blue line), almost all X molecules taken up are secreted to the external environment via reaction R 1 .Y is also taken up from the external environment while consuming 3ATPs.Note that reaction R 0 consumes a single ATP molecule to take up a single X molecule, and the reaction R 1 produces a single ATP molecule.For the model to obtain a net gain of ATP, X must be converted into Y by consuming one more ATP molecule via reaction R 2 .However, once the model falls into the attractor where ATP is depleted (endpoint of the blue line in Fig. 3B), it cannot invest ATP via R 2 , and the model is "dead-locked" in the state.Hereafter, we refer to the former and latter attractors as the active attractor x A and inactive attractor x I , respectively.Now, we compute the dead states of the model according to Def. 3. In this framework, we first need to select the representative living states.Here, we suppose that the active attractor x A is the representative living state.According to Def. 3, we compute the complement set of the stoichiometric paths of the active attractor, D(x A ) = R N >0 \SP(x A ).This complement set corresponds to the dead set with x A as the representative living state.Note that choosing x A as the repre-sentative living state is equivalent to choosing mutually reachable attractors from/to x A by continuously changing the enzyme activities and their basin of attraction (Recall Eq. ( 18)).By selecting a single-point attractor as the representative living state, we can effectively select a large region in the phase space.
Due to the challenges associated with computing the stoichiometric paths, we utilize the stoichiometric rays for the computation.Let us term D(x A ) = R N >0 \SR(x A ) the underestimated dead set.As the stoichiometric rays is a superset of the stoichiometric paths, the underestimated dead set is a subset of the dead set.Thus the stoichiometric rays does not show a false negative of returnability; If x ∈ D(x A ) holds, then x is guaranteed to be dead with respect to the active attractor.In Fig. 4A, the underestimated dead set with a finite L, ) is depicted in the phase space (For the computational procedure, see SI text section S3 and SI Codes).The points inside the underestimated dead set cannot return to the active attractor regardless of how the enzyme activities are modulated with a given number of possible flips.The inactive attractor is contained in the underestimated dead set, and thus the inactive attractor is indeed the "death" attractor at least with the L = 4 flips 8 .Also, the boundary of the under- 8 We computed the underestimated dead set for L = 6 with a lower estimated dead set ∂ DL (x A ) now gives an estimate of the SANZ surface.
In Fig. 4B, we plot the points in the underestimated dead set with an analytically estimated SANZ surface of the set on the 2-dimensional slice of the phase space where z value is set to a single value.The SANZ surface is calculated from the thermodynamic and stoichiometric constraints that each single stoichiometric ray should satisfy.The curves match well with the boundary of the underestimated dead set (for the detailed calculation, see SI Text section S4).

B. The global transitivity of the toy metabolic model
The application of the stoichiometric rays is not limited only to computations of the dead set with pre-defined representative living states.The stoichiometric rays enables us to investigate the mutual transitivity among states without a priori definition of living states.In this section, we present a transition diagram of the toy metabolic model.The computation of the transition diagram corresponds to the computation of the equivalent classes and partial order ( , R 3 >0 /⇌) defined in the framework section (Sec.II).The transition diagram thus dictates the global transitivity of the model.This provides insights for choosing the representative living states, and in addition, criteria for the model selection suitable for investigating cell death.
For the computation of the transition diagram, we sampled several points in the phase space (Fig. 5A) in addition to the active and inactive attractors, and computed controllability for all pairs of the sampled points.Note that we computed only the controllability to the active attractor in the previous part, but here we compute the controllability of all pairs of points mutually.With the computation, we construct the directed graph G L where the nodes are the sampled points and the directed edges represent the controllability from the source point to the target point, i.e, we put the directed edge from x to y if x ∈ SR L (y) holds.
Next, we computed the condensation of G L which, denoted GL .Condensation is a directed graph in which strongly connected components (SCCs) of G L are contracted to a single node.The nodes of G L within the same SCC are mutually reachable, and thus, each node in the condensation graph GL corresponds to equivalent class (see Sec. II).Hereafter, we refer to the condensation graph GL as a transition diagram.
The transition diagram is shown in Fig. 5B, describing the partial order relation on R 3 >0 /⇌.In addition, because the condensation graph is a directed acyclic graph, resolution of the point sampling due to the computational cost.We confirmed that the increase of L from 4 to 6 does not make the underestimated non-returnable set larger.(see Fig. S3).
the diagram highlights the hierarchy of the state transitions or the "potential" of the states.The active and inactive attractors are in the SCCs represented by the yellow-and pink squares in Fig. 5B, respectively.Let us refer to the pink SCC as the terminal SCC, C * , since it has only in-edges, and thus, the states in C * = Ω(R 3 >0 ) (cf.Def. 4) cannot escape from the set regardless of the control u(t).The transitivity among states with a given number of flips is fully captured by the transition diagram while we carried out the computations of the transitivity for a limited number of points in the phase space due to the computational cost.
The computation of the transition diagram provides a guide for selecting the representative living states, and in addition, a criterion for the model selection suitable for investigating cell death.It is commonly believed that controlling to the living states is difficult and possible only from limited states.Also, all living and non-living systems can be non-living (dead state), i.e., control to the dead state should always be possible regardless of the source state.If we accept these beliefs, the representative living state should be taken from the SCCs except for the terminal SCC.When we additionally require that the representative living state be the attractor (i.e., the state is stable without temporal modulation of enzymatic activities), the representative living state should be taken from the SCC containing the active attractor (the yellow square in Fig. 5B).
In the toy metabolic model, the active attractor and the inactive attractors are in the non-terminal SCC and the terminal SCC, respectively.However, this is not always the case for any biochemical model.In order to study cell death employing mathematical models, we may require the models to satisfy conditions such as the ones above.The transition diagram is useful method for checking if a given model satisfies requirements as the diagram captures the global transitivity in the model.

V. DISCUSSION
Questions such as "what is life?" and "what is death?"relate to how we can characterize living and dead states.
Providing a straightforward answer for such a characterization is challenging even for mathematical models of cells.However, quantitative evaluation of the "life-death boundary" (SANZ hypersurface), or in other words, characterization of the difference between life and death can be possible.In the present manuscript, we proposed a theoretical framework to study the dead set and SANZ hypersurface with the development of a tool for computing it.
In the present framework, the dead state is defined based on the controllability to the representative living state.Thanks to the useful mathematical features of controllability relationship, x ⇀ y, the death judgment is robust on the choice of the representative living state; The choice of representative living state is arbitrary, as long as they are chosen from the same equivalence class.Further, Thm. 1 states that even if there are uncountably many living states, the representative living states can be countably many, or even, a finite number.
Our strategy for defining the dead state is to use the law of exclusion, of the form "what is not living is dead".This allows us to avoid direct characterization of dead states.In this sense, our characterization of death is similar to that of the regrow experiment, rather than the staining assay in which the dead state is characterized by stainability.For such indirect characterizations, it is necessary to determine whether a cell will eventually regrow.This is an extremely hard challenge in experimental systems.However, by focusing only on mathematical models, we can determine, in principle, whether the cell will regrow.The stoichiometric rays provides a reasonable method for computing this.
The heart of the stoichiometric rays is that the modulation of enzymatic activity cannot directly change the directionalities of the reactions.As a consequence of this feature, the null-reaction manifolds and direction subsets are invariant to enzymatic activities.This allows us to efficiently compute the stoichiometric rays.The usefulness of the stoichiometric rays was demonstrated by using two simple nonlinear models of chemical reactions.Especially, in the toy metabolic model, we showed that the stoichiometric rays is a powerful tool for computing the dead set in which no control can return the state to the active attractor with a given number of flips.
Thermodynamic (ir)reversibility is believed to be key to grasping the biological (ir)reversibility such as the ir-reversibility of some diseases and death.However, there is no system-level link between thermodynamic and biological reversibilities except in a few cases where simple causality leads to biological irreversibility such as in transmissible spongiform encephalopathy [34].The stoichiometric rays enables us to bridge the (ir)reversibility of each reaction to the system-level reversibility; Recall that the SANZ surface of the dead set is calculated from the thermodynamic constraint, i.e., directionality of the reactions, and the stoichiometric constraint.Further studies should pave the way to unveil the relationship between thermodynamic irreversibility and biological irreversibility.This may reveal the energy (or "negentropy" [35]) required to maintain cellular states away from the SANZ hypersurface [36].
Finally, we would like to mention an extension of the present framework to the probabilistic formulation and its implications.In the present deterministic version, the transition from the dead state to the living state is irreversible.In stochastic processes, however, the transition from the state considered dead in the deterministic model to the living states is possible with some probability.Thus, in probabilistic formulation, life and death are not binary categories, but are continuously linked states.
Given the small size of cells, the probabilistic formulation is more suitable for describing cellular state transitions.Death is commonly believed to be irreversible.However, if the theory of death reveals the probabilistic nature of life-death transitions, it may have multiple implications for reforming the concept of "death".
(For the computational procedure, see SI text section S1 and SI Codes).In the figures, SR L (x * )\SR L−1 (x * ) is filled with distinct colors.As shown in Fig 2A, it is possible to reach the point x * = (0.6, 0.6) from anywhere in the phase space, SR L (x * ) = R 2 >0 (L ≥ 3), while we need to choose the starting point from the specific region for the other target point x * = (1.5, 1.5) as in Fig. 2B because SR L (x * ) R 2 >0 (L ≥ 0) holds for the point.

FIG. 1 .
FIG. 1.An illustrative example of the stoichiometric rays overestimates the controllable set.(A) All the paths from the pink point are blocked by the null-reaction manifold M k .However, Def. 6 judges that there are rays from the pink point to the highlighted region on Mj. (B) An example case where a part of the single stoichiometric rays are blocked, but the stoichiometric rays correctly judge the reachability.

FIG. 4 .
FIG. 4. (A) The active attractor (red point), the inactive attractor (blue box), and the points in the underestimated dead set (black dots) are plotted.The gray rectangle is the convex hull of the dead states.The null-reaction manifolds M0, M2, and M4 are drawn as the salmon, green-blue, and olivecolored surfaces.The null-reaction manifold M5 is not illustrated because it is not necessary for the argument.The intersections of the null-manifold pairs (M0, M2), (M2, M4), and (M4, M0) are represented by the solid, dash, and dotdash curves, respectively.(B) An analytic estimate of the boundary of the underestimated dead set (red curve) and the points in the underestimated dead set (green-blue points) are plotted.log 10 z value is fixed to −1.7.We set L = 4.The same k values with Fig. 3 are used.

FIG. 5 .
FIG. 5.The transition diagram of the discrete states of the toy metabolic model.(A) The points for which the mutual transitivity is computed.The points are grouped and labeled by the strongly connected components (SCCs).The dashed lines are eye guides.The active attractor (yellow square) and the inactive attractor (pink square) are shown together.Note that the green triangle is on the grid, but not the active attractor.(B) The transition diagram of the SCCs.Only the two SCCs (pink-and yellow squares) contain more than one state.The terminal SCC (pink square, indicated by the reddashed rectangle), C * , has only in-edges.The same k values with Fig. 3 are used, and L = 4.