{'fraction_non_alphanumeric': 0.03451248218015449, 'fraction_numerical': 0.00586811656665451, 'mean_word_length': 4.6234153616703955, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Formal characterization and efficient verification of a biological robustness property

Robustness is an observable property for which a chemical reaction network (CRN) can maintain its functionalities despite the influence of different perturbations. In general, to verify whether a network is robust, it is necessary to consider all the possible parameter configurations. This is a process that can entail a massive computational effort. In the work of Rizk et al., the authors propose a definition of robustness in linear temporal logic (LTL) through which, on the basis of multiple numerical timed traces obtained by considering different parameter configurations, they verify the robustness of a reaction network. In this paper, we focus on a notion of initial concentration robustness ($\alpha$-robustness), which is related to the influence of the perturbation of the initial concentration of one species (i.e., the input) on the concentration of another species (i.e., the output) at the steady state. We characterize this notion of robustness in the framework proposed by Rizk et al., and we show that, for monotonic reaction networks, this allows us to drastically reduce the number of traces necessary to verify robustness of the CRN.


I. INTRODUCTION
Two main characteristics define living cells: an intrinsic structural complexity, according to which there is modularity inside the cell itself, and the ability to interact with other networks, working as a system. In addition to this already intricate framework, at various frequencies and timescales, internal and external fluctuations can alter specific functions or traits of biological systems, causing genetic mutations, loss of structural integrity, diseases and so on. Nevertheless, many biological networks can maintain their functionalities despite perturbations: this distinct property is known as robustness [1].
Robust traits are pervasive in biology: they involve various structural levels, such as gene expression, protein folding, metabolic flux, species persistence. For this reason, the study of robustness is essential for biologists, whose aim is to understand the performance and functions of a biological system.
However, it is not easy to investigate biological systems since they often exhibit non-linear and non-intuitive behaviors. They can be studied by performing wet-lab (in vitro) experiments, or through mathematical or computational (in silico) methods on a pathway model. Unfortunately, the applicability of the last approach is often hampered by the complexity of the models to be analyzed (often expressed in terms of ODEs or Markov chains). An alternative way is to infer properties of the system by analyzing only its structure, without studying or simulating its dynamics, as proposed in [2]. However, this approach has been shown to be limited in its applicability because the derived conditions are often very restrictive.
Robustness can be formally studied by applying the methodology proposed by Rizk et al. in [3], [4]. Such a methodology is based on the definition of robustness given by Kitano in [1] as the ability of a system to maintain specific functionalities against perturbations. The robustness of a system is measured as the distance between the system behaviour under perturbations and its reference behaviour expressed as a linear temporal logic (LTL) formula. The distance is computed by using a notion of violation degree that measures how much the temporal logic formula should be changed in order to match traces of perturbed behaviours obtained, for instance, through many simulations.
The approach proposed by Rizk et al. is very general, both in the description of the reference behaviour and in the kind of perturbations considered. On the contrary, we want to focus on a particular notion of robustness, namely the initial concentration robustness that studies the influence of the initial concentrations of species on the concentration of the other species at steady state of the system. For these reasons, we consider the notion of α-robustness, based on continuous Petri nets [5] and interval markings, which extends the notion of absolute concentration robustness considered in [2], [6].
The evolution of a biological system under initial concentration perturbations can be understood by studying the simulations of the systems under all possible combinations of its initial concentrations, a process that generally entails a massive computational effort. In order to drastically reduce the number of simulations needed for this study, we proved that -if the concentration of an output species is monotonic with respect to the concentration of an input species (i.e. the perturbed element) -the number of simulations can be reduced to two by studying the model only on the extreme values of the input concentration range. This result was the subject of the Input-Output theorem on monotonicity [7].
In this paper, we show that α-robustness, proposed by Nasti et al. [8], is a particular instance of the robustness notion proposed framework proposed by Rizk et al., and, therefore, it can be characterized in their framework. The great advantage in case of a CRN showing a monotonic behaviour according to our Input-Output theorem, this allows us to drastically reduce the number of traces necessary to verify robustness.
The paper is organized as follows. We proceed by first introducing the notion of Chemical Reaction Networks (CRN) in Section II-A. In Section II-B, we give the formal definition and the semantics of Linear Temporal Logic, which represents the base of the work done by Rizk et al, described in Section III. In Section IV, we give the formal definitions of initial concentration robustness together with the characterization of α-robustness in terms of the framework proposed by Rizk et al.. We also characterize monotonicity in CRN. We apply our approach to the biological example described in Section III-C. Finally, in Section V we draw some conclusions.

II. BACKGROUND
We introduce some notions that will be assumed in the rest of the paper. In the first part, we focus on the representation of chemical reactions, considering one of the main methods that we can use to describe them: the deterministic approach.
In the second part, we present the Linear Temporal Logic formalism that provides a mathematical notion to express systems behavior, and it is at the basis of the work done by Rizk et al..

A. Chemical Reaction Networks
A chemical reaction is a transformation that involves one or more chemical species, in a specific situation of volume and temperature.
We call reactants the chemical species that are transformed; while those that are the result of the transformation are called products. We can represent a chemical reaction as an equation, showing all the species involved in the process.
A simple example of chemical reaction is the following elementary reaction: In this case, A, B, C, D are the species involved in the process: A and B are the reactants, C and D are the products. The parameters a, b, c, d are called stoichiometric coefficients and represent the number of reactants and products participating in the reaction. The arrow is used to indicate the direction in which a chemical reaction takes place. When we have only one arrow, it means that the reaction is irreversible, that is it is not possible to have the opposite process. To describe the dynamical behaviour of the chemical reaction network, we can use the law of mass action, which states that: the rate of a reaction is proportional to the product of the reactants. Applying the law of mass action to the system, we obtain, for each chemical species, a differential equation describing the production and the consumption of the considered species. Considering the generic chemical equation 1, we obtain: where, in each equation, we isolated the term describing the direct reaction from the one describing the inverse reaction. With these two terms, we implicitly considered, for each element, the processes of consumption and production.

B. Formal definition of LTL
Linear Temporal Logic is a logical formalism. It provides a mathematical notion to express systems behaviors [9], based on a linear-time perspective. The temporal logic is necessary to specify the relative order of events, expressed by elementary modalities, which combined can express complex dynamical properties. Typical properties are oscillations (when a behavior recurs infinitely), reachability (when the system can reach a given state), invariance (when a property is always true), inevitability (when a system has to reach a given state), response (an event causes a specific behavior) [10].
A basic LTL formula φ consists of atomic propositions a ∈ AP , Boolean connectors (∧, ∨, ¬, =⇒ ), and two basic modal operators: • Xφ or • ("next") means that a given formula φ is true in the next state; • φUφ or ∪ ("until") means that given two formulas φ 1 and φ 2 , the formula φ 1 is true, until the formula φ 2 becomes true. The "until" operator ∪ allows to derive other two temporal modalities, defined as follows: • Fφ or "eventually (in the future)" means that a given formula φ is true now or sometime in the future, defined as φ def = true ∪ φ; • Gφ or "globally" means that a given formula φ is true now and forever, defined as φ def = ¬ ¬φ. The LTL formulae are formed according to the following grammar: The atomic proposition a, with a ∈ AP is a state label, representing an assertion about the value of a system variable that has to be evaluated, such as the concentration of a chemical species.
A LTL formula φ represents a property of a trace, which is a infinite path. Given a path and a formula φ, we can formulate precisely when φ holds on the path. For example, the trace T in Figure 1 satisfies the formula φ 1 = F(x ∨ ¬y) because it is true in the first state of the trace T . Instead, the trace does not satisfy the formula φ 2 = G(y) because y is not true in all states.

1) Semantics of LTL:
To precisely formulate when a path satisfies an LTL formula, we define the semantics of LTL formula φ by providing a satisfaction relation |= such that σ |= φ if and only if a property φ is satisfied by a trace σ. [9]: Definition 1: Let φ be an LTL formula over AP and σ ∈ (2 AP ) ω be a trace. The satisfaction relation |= ⊆ (2 AP ) ω × LT L is the smallest relation with the following properties: For the derived operator and the expected result is: In Figure 2, we add an intuitively sketch of the semantics of temporal modalities.

A. Temporal logic semantics of numerical traces
Numerical simulations are used to obtain the behavior of a biological system, which is described by a numerical timed trace as shown in [3]. A numerical timed trace, expressing the evolution of a system with time, is a finite sequence of By a numerical trace, we can depict different biological phenomena, such as the time evolution of a concentration of a chemical species in a system, as represented in Figure 3. In this example, T = ((0, 2, 0), (1, 6, 4.12), ..., (9, 10, 0)) is the the associated trace, in which each state variable defines a specific concentration level of the species B over time. In [3], the authors use LTL to express dynamical properties of biological systems. LTL operators describe if and when a property φ holds on a trace T . Considering again Figure 3, the formula φ = F([B] > 7) expresses that at some point the concentration of species B is greater than 7.

B. The formal definition of Robustness degree
Since it is interesting to define how much a numerical trace satisfies a formula φ, the authors specify the quantifier-free LTL (QFLTL), which replaces the numerical constants in the atomic propositions of a LTL formula, with free real-valued variables y. Then, in this way, having a formula φ and a vector of real-valued variables y, it is possible to know for which values y the QFLTL formula φ(y) holds on T . At this point, the satisfaction domain is defined as follows: Definition 2 (Satisfaction Domain): Given a QFLTL φ formula, for any trace T , the satisfaction domain of φ(y) is the set of variables y for which φ(y) holds. It is defined as: where q is the number of constants appearing in φ.
Through this approach, the LTL formula becomes an instance of a more general QFLTL formula obtained by variable abstraction. Considering again Example 3, and the formula Moreover, concerning trace T in Example 3, the domain is D T,φ(y1,y2) = {y 1 ≤ 10 ∧ y 2 ≥ 2}, since 2 and 10 are respectively the minimum and the maximum values of the trace.
Given a trace T = (s 0 , s 1 , ..., s n ) and a LTL formula φ, the authors define the notion of violation degree to quantify how much φ must be changed to hold on T . This concept is defined as the Euclidean distance between a formula φ and the domain of the trace T . It is formally defined as follows: Definition 3 (Violation Degree): The violation degree vd(T, φ) of a formula φ with respect to a trace T is the distance between the actual specification and validity domain D T,φ((y)) of the QFLTL formula φ(y) obtained by variable abstraction: Considering Example 3 and the formula , the violation degree is vd(T, φ 2 ) = 2 meaning that φ 2 has to be changed to hold on T .
To define, instead, how much the given LTL formula holds on a given numerical trace, the notion of satisfaction degree is introduced as follows: Definition 4 (Satisfaction Degree): The satisfaction degree sd(T, φ) of a formula φ with respect to a trace T is defined as: where vd(T, φ) represents the violation degree. The value obtained by computation of the satisfaction degree ranges between 0 and 1. The satisfaction degree is equal to 1 when the trace T satisfies a formula φ, otherwise it tends to 0. The concept of satisfaction degree characterizes the definition of robustness.
Definition 5 (Rizk et al. Robustness): The robustness of a system is defined as: where φ is the specification of the functionality in LTL; T p is the numerical trace, representing the system behavior under perturbation p; P is the set of perturbations. The continuous probability distribution characterizes the perturbations, affecting the entire system: each perturbation has its weight, representing how much it can influence the biological behavior under study.
This entire approach is already implemented in BIOCHAM (BIOCHemical Abstract Machine) [11], a software environment for modeling biochemical systems, where the function computing robustness is the following: [variable objective]).
As described in [11], the function computes the robustness degree with respect to the LTL Formula, for the list of [parameters name] and with list of objectives for the free variables of LTL Formula given in [variable objective].
The following two functions can be used respectively to set the coefficient of variation of the parameters and the number of simulations: • option(robustness_coeff_var: number); • option(robustness_samples: integer).

C. Example: Application to the ERK signaling pathway
As example, we show how to analyze the dynamics of the ERK signalling pathway.
A signalling pathway consists of enzymatic cascades, having a starting species that triggers the other connected reactions. In general, there is a particular species, namely the transductor that perceives an initial stimulus, which activates the cascade amplifying the signal for the next enzymatic process.
One of the most important examples of such processes is the ERK pathway, involved in many biological phenomena such as cell's growth and differentiation. We consider a particular portion of the mathematical model of the ERK pathway (denoted as ERK * ) implemented by Schilling et al. [12] and available to the public on the BioModels Database (BIOMD0000000270). We indicate the species and the kinetics rates as originally denoted in the model in [12]. For simplicity, we refer to the reaction using the notation R i , where i is the kinetics rate index. The reactions involved are the following: in Table I we reported the coefficient rates and the initial conditions of ERK * system. The species PRaf, involved in the reaction R 21 , acts as catalyst promoter, which means that its concentration positively influences the production of the species PMek1. A catalytic species increases the reaction rate of the reaction in which it is involved, and during this process, its concentration is not consumed. To verify the robustness of the system, we need to perform many simulations, considering different kinds of perturbations. For instance, if we want to test whether the species PPMek1 is robust with respect to oscillations on the initial concentration of Raf, we need to compute many simulations, one for each possible (continuous) value of the initial concentration of Raf. In Figure 4, we show how the concentration of the species PPMek1 varies at the steady state with respect to the initial concentration of the species Raf. We notice that PPMek1 does not exhibit absolute robustness.

IV. THE DEFINITION OF THE α-ROBUSTNESS AND ITS CHARACTERIZATION IN THE RIZK FRAMEWORK
The definition of robustness degree given by Rizk et al. represents the average satisfaction degree of the property of interest over all admissible perturbations, possibly weighted by their probabilities [4].
In our case, as already described by Nasti et al. in [8], we want to focus on the evaluation of the initial concentration robustness. We want to vary the concentration of at least one chemical species (namely the input) and verifying, at the equilibrium, if the concentration of another species (namely the output) is included in an interval of possible values.
In order to do that, we recall the definition of the initial concentration robustness, using continuous Petri nets.
Definition 6 (Continuous Petri net): A continuous Petri net N can be defined as a quintuple P, T, F, W, m 0 where: • P is the set of continuous places, conceptually one for each considered kind of system resource; • T is the set of continuous transitions that consume and produce resources; represents the set of arcs in terms of a function giving the weight of the arc as result: a weight equal to 0 means that the arc is not present; • W : F → R ≥0 is a function, which associates each transition with a rate; • m 0 is the initial marking, that is the initial distribution of tokens (representing resource instances) among places. A marking is defined formally as m : P → R ≥0 .
Tokens are movable objects, assigned to places, that are consumed by transitions in the input places and produced in the output places. Graphically, a Petri net is drawn as a graph with nodes representing places and transitions. Circles are used for places and rectangles for transitions. Tokens are drawn as black dots inside places. Graph edges represent arcs and are labeled with their weights. To faithfully model biochemical networks, the marking of a place is not an integer (the number of tokens) but a positive real number (called token value representing the concentration of a chemical species. Each transition is associated with a kinetic constant, that determines the rate of (continuous) flow of tokens from the input to the output places of the transition.
In order to give the definition, we recall some definitions introduced by Nasti et al. in [8]. The initial marking is defined as an assignment of a fixed value to each place p. Now, it is possible to generalize the idea of initial marking by considering a marking as an assignment of a interval of values to each place p of the Petri net.
We first recall the definition of the domain of intervals. An interval marking in which at least one interval is nontrivial represents an infinite set of markings, one for each possible combination of values of the non-trivial intervals. Therefore, given an interval marking, we relate it with the markings as in the original Petri nets formalism in the following way: In a Petri net we assume that there exists at least one input place and exactly one output place representing input and output species of the modeled biochemical network, respectively. Under this assumption, we can give the formal definition of robustness.
Definition 9 (α-Robustness): A Petri net N with output place O is α-robust with respect to a given interval marking m [ ] iff ∃k ∈ R such that ∀m ∈ m [ ] , the marking m corresponding to the steady state reachable from m, is such that A. The initial concentration robustness in the general Rizk's framework The specific notion of the initial concentration robustness is analyzable in the general framework proposed by Rizk et al..
To match the two frameworks, we restrict the analysis of the system behavior at the equilibrium. We show that Definition 9 is an instance of Definition 5, expressed using the Rizk et al. framework shown in [3].
Theorem 1 (Definition 9 in the context of Definition 5): Given a Petri Net P N with output place O, an initial marking m [ ] ∈ M [ ] and a continuous probability distribution prob(m), defined on m [ ] such that the integral of the pdf is normalized to 1. P N is α-robust with respect to m [ ] iff there exists an interval [min, max] ∈ R such that R s φ,P = 1, Hence, according to Theorem 1, we define the initial concentration robustness as: We can implement the Formula 7 in BIOCHAM: where: • robustness is the function, implemented in BIOCHAM [11], computing the robustness measure; • In is the concentration of the input species, which is perturbed; • O is the concentration of the output species; • F and G are two specific temporal operators, respectively meaning future and globally; • the expression [min -> x, max -> y] represents the assignment of the interval limits; <= max)) represents the system behavior that has to be evaluated. The formula expresses that, when the system reaches the steady state, the concentration of output species is within a range. The definition of α-robustness is simpler and much less general than the one considered by Rizk et al.. However, it is conceived with the aim of enabling further studies on sufficient conditions that could allow robustness to be assesses by avoiding (or significantly reducing) the number of simulations to be performed. This could be obtained, for instance, by adapting conditions already considered in the context of monotonicity analysis [13], as we will describe in the next Section.

B. Input-Output Monotonicity in CRN
In general, the verification of robustness requires a huge number of simulations [2], being necessary to test the system behavior for all the possible combinations of initial concentrations of the involved chemical species.
In order to reduce the computational effort, we apply a sufficient condition that allows us to analyze the monotonicity between the output of the system (the under study chemical species at the steady state) with respect to the initial concentration of the input (the perturbed chemical species) [7]. In this context, the input and output are in a monotonicity relation if the concentration of output at any time either increases or decreases due to an increase in the initial concentration of the input, within an interval [min, max] ∈ R: Definition 10 (Positive (respectively Negative) Input-Output Monotonicity): Given a set of reactions R, species S O is positively monotonic (resp. negatively monotonic) with respect to S I in R if, for any two initial states S 0 , S 0 as above, If the monotonicity is established, we can consider only the extreme values of the input concentration interval, as described in [14] by Gori et al.. Indeed, as we will see in detail in Section IV-C, if the output is monotonic with respect to the input, then just two simulations are necessary: one with input = min and one with input = max. The dynamics of the output in all the other (intermediate) cases is included in the results we obtained from these two configurations.
The proposed condition is based on a constraint on the structure of the chemical reaction network that can be evaluated efficiently, without the need of performing simulations. Following the guidelines of [13], it is based on a graph representation of the chemical reaction network, namely the R-graph, enriched with information about cooperation and competition among reactions (i.e., the Labelled R-graph).
The R-graph is formally defined as follows. Definition 11 (R-graph): Given a finite set of reactions R over a set of species S, the R-graph of R is the signed graph R, E + , E − , where E + ⊆ (R × R) and E − ⊆ (R × R) are defined as follows: there is a species which is a reactant (or product) in both reactions R i and R j . Intuitively, given two reactions, we draw a positive edge between them if they cooperate each other, hence, for example, the product of a reaction is among the reactants of the other reaction. Instead, we draw a negative edge between two reactions if both share the same reactants, hence they compete.
In the R-graph, essentially, each edge is associated with one of the symbols + or −. We can label in the same way also its vertices, assigning to each one a sign + or −: Definition 12 (Consistent labeling): A map σ : R → {+, −} is called a consistent labeling of the R-graph if the following two properties hold: The condition is then expressed as a set of constraints on the graph structure. In addition, studying the signs of the stoichiometric matrix of the network, it is possible to predict if a variation on the concentration of the input affects the output positively or not, as described in [7], where we find the main result: Theorem 2: Given a set of chemical reactions R, and two species S j I (the input species) and S j O (the output species, with with j O = j I ). If the following three conditions hold: 1) the R-graph of R admits a consistent labeling σ; 2) the input species S j I is involved in exactly one reaction R i I ∈ R; 3) the output species S j O is involved in exactly one reaction R i O ∈ R; then, the species S j O is positively monotonic with respect to have opposite signs, and negatively monotonic if they have the same sign.
We refer to [15] for the proof of Theorem 2.
C. The Input-Output monotonicity applied to the ERK signaling pathway Consider again the example of the ERK signaling pathway, described in Section III-C.
We now apply the definition of the α-robustness and we identify Raf and PPMek1 as the input and the output of the network. We assume that the initial concentrations of all the species are fixed, but that the initial concentration of species Raf can vary from 1 to 100.
In order to reduce the computational effort of simulation, we apply the sufficient condition of the Input-Output monotonicity. Indeed, if species PPMek1 is monotonic with respect to Raf, then just two simulations are necessary: one with Raf = 1 and one with Raf = 100. The dynamics of PPMek1 in all the other (intermediate) cases is included in the results we obtained from these two simulations.
As shown in detail in [7], [15], we notice that the species PRaf is the link between the two groups of reactions: in the first reaction, PRaf acts as the output of the network, while in the second one it acts as a kinetic constant. Since the reactions follow one another in a chain, we assume that PRaf reaches the steady-state before the PPMek1 species, which represents the final product of the CRN, and we justify the assumption by simulations, as shown in Figure 5.
For the second sub-network it would be natural to select PRaf as the input of the network, since it represents the link between the first and the second block. However, we cannot choose PRaf as input, since it does not appear in the role of reactant or product. Therefore, we choose Mek1 as input. Indeed, we know that the concentration of [P Raf ] increases the rate of the reaction in which Mek1 is the only reactant. Consequently, to verify if PPMek1 is monotonic with respect to PRaf, we can verify if this species is monotonic with respect to the Mek1. Therefore, we choose Mek1 and PPMek1, respectively, as the input and the output. Then, we proceed by building the labelled R-graph, represented in Figure 6, and computing the following stoichiometric matrix: By calculating the products Γ j I i I σ(i I ) and Γ j O i O sσ(i O ), which have opposite signs, we find that the species PPMek1 is positively monotonic with respect to the species Mek1. As consequence, since PRaf is positively monotonic with respect to Raf and PPMek1 is positively monotonic with respect to Mek1, PPMek1 is positively monotonic also with respect to Raf, as we already noticed in Figure 4.
This result enables us to use just two simulations (i.e., Raf= 1 and Raf= 100) to verify if the PPMek1 is robust with respect to the variation of the initial concentration of Raf. V. CONCLUSIONS Robustness is a crucial feature of many biological systems because this property allows their correct functioning in presence of molecular noise and environmental fluctuations.
In order to verify the system's robustness, many strategies have been proposed. Among them, one of the most common approaches is to simulate the system with all possible combinations of initial concentrations of chemical species. However, this requires a considerable number of simulations (in general an infinite number).
In [3], Rizk et al. propose a general and computational framework for the definition of the robustness of biological functions with respect to a set of perturbations, based on LTL, an expressive language for specifying dynamical behaviors widely used in computer science and engineering. Using this framework, implemented in BIOCHAM, they are able to describe on average how the system behaves under perturbations.
In this paper, we show that we can cast our formal definition of the initial concentration robustness in the general framework proposed by Rizk et al.. Then, we can apply the sufficient condition of the Input-Output monotonicity in chemical reaction networks to substantially reduce the computational effort required by the huge number of simulations needed to prove the robustness property.