Incorporating Time Delays in Process Hitting Framework for Dynamical Modeling of Large Biological Regulatory Networks

Modeling and simulation of molecular systems helps in understanding the behavioral mechanism of biological regulation. Time delays in production and degradation of expressions are important parameters in biological regulation. Constraints on time delays provide insight into the dynamical behavior of a Biological Regulatory Network (BRN). A recently introduced Process Hitting (PH) Framework has been found efficient in static analysis of large BRNs, however, it lacks the inference of time delays and thus determination of their constraints associated with the evolution of the expression levels of biological entities of BRN is not possible. In this paper we propose a Hybrid Process Hitting scheme for introducing time delays in Process Hitting Framework for dynamical modeling and analysis of Large Biological Regulatory Networks. It provides valuable insights into the time delays corresponding to the changes in the expression levels of biological entities thus possibly helping in identification of therapeutic targets. The proposed framework is applied to a well-known BRNs of Bacteriophage λ and ERBB Receptor-regulated G1/S transition involved in the breast cancer to demonstrate the viability of our approach. Using the proposed approach, we are able to perform goal-oriented reduction of the BRN and also determine the constraints on time delays characterizing the evolution (dynamics) of the reduced BRN.

1. Introduction

Biological Regulatory Networks (BRNs) are used to describe almost all biological functions to cover the interactions taking place in various chemical reactions in living organisms. The nature of these interactions is continuous and stochastic in nature and the rate of change in the underlying dynamics of these interactions is variable. A large number of formal approaches have been devised to model the topology of BRNs and to analyse their dynamical behaviors (De Jong, 2002; Sheikh et al., 2016). BRNs are traditionally modeled by differential equations but they are non-linear and thus difficult to solve analytically. Moreover, the available experimental data is mostly qualitative in nature and do not aid in precise determination of quantitative parameters for the differential equations. A powerful tool involving Lie algebra has been reported to solve a large class of non-linear epidemic spreading systems analytically (Shang, 2012, 2015a), however it has not been applied to the dynamical modeling of BRNs. René Thomas Discrete Modeling technique (Thomas, 1978) proposed a qualitative method for representing the change in expression level of a gene or biological entity with a logical function having discrete values. The information on whether the gene expression would increase or decay is contained in the K-parameters associated with each entity of the BRN. The values of these parameters determine the dynamical behavior of the system that is represented by a State Transition Graph (STG). The state graph grows exponentially with the increase in the number of entities thus it becomes difficult to analyse medium or large scale BRNs (containing more than 10 or more genes). Moreover, this approach ignores the time taken for the gene to reach the expression level required for production or degradation and therefore, uses an asynchronous evolution of the dynamics of the system. Hybrid Modeling Technique (Ahmad et al., 2006) addressed this shortcoming by associating time delays in BRNs by considering the change in qualitative levels as piece-wise linear functions. This made it possible to perform model checking and obtain constraints in terms of production and/or degradation delays. The hybrid approach works well for BRNs consisting up to 7–8 genes, however, it is unable to compute the time delays for larger networks due to the state space explosion.

A Process Hitting framework for analysis of large BRNs has been proposed (Paulevé et al., 2011; Folschette et al., 2012) which can handle networks comprising of thousands of genes. The very high scalability of this approach is due to the following factors:

• It takes into consideration only the most permissible (generalized) dynamics possible in the interaction graph of a BRN instead of dealing with the whole state space.

• The generalized dynamics is refined with the help of cooperativity between two or more genes, which have a combined influence on any other gene in the network.

The Process Hitting approach is thus able to quickly perform static analysis like determination of stable states, successive reachability and inferring the K-parameters of the BRN. This technique is based on Stochastic π-Calculus and allows for synthesis of temporal and stochastic parameters, which enables the simulation of dynamical behavior to some extent. However, there is no mechanism for incorporating the time delays in the proposed framework. Lately, time parameters have been introduced into Process Hitting with classes of priorities (Folschette et al., 2015); however, it still does not infer time delays of the interactions.

The Process Hitting is also able to perform goal-oriented reduction of the large BRNs by taking into consideration the minimal traces that are necessary to achieve the desired reachability goal. It applies the cutsets to preserve the minimal traces through Gene Knock-out/in/down technique (Paulevé et al., 2013).

Previously, enhanced modeling of BRNs based on timed automata has been performed in Siebert and Bockmayr (2008) which introduced time delays in the René Thomas Discrete Modelling Framework and it is possible to simulate the dynamics of BRN. However, this approach does not infer the constraints on time delays rather the numerical values have to be adjusted manually to obtain a certain behavior. Moreover, it is also limited to small BRNs as it introduces intermediate states, which results in an even greater number of possible states in the state transition graph.

Recently, new methodology has been proposed for the inference of BRNs through a time extension of Automata Networks using the time series data and known influences among the genes (Ben Abdallah et al., 2017). This modeling approach requires the observed experimental time series data for deducing the BRN model which satisfies the dynamical behavior depicted by the observed data.

Our approach is to extend the Process Hitting Framework by introducing Time Delays in it so that it becomes possible to determine the constraints on the activation and degradation delays associated with the evolution of a gene in the dynamical model of the BRN. In order to achieve it, we represent the biological entities of the BRN with Biological Linear Hybrid Automata (Bio-LHA) (Ahmad et al., 2009; Aslam et al., 2014) and enrich them with time delays while the logical rules for gene activation or degradation are provided by the cooperative hits as determined by Process Hitting. Since PH only takes into consideration the permissible (possible) dynamics of the BRN as given in its Interaction Graph (IG), therefore, we have a reduced STG to be represented by the Hybrid Automata.

The proposed approach, Hybrid Process Hitting, has been applied to a toy BRN to fully explain the steps involved and subsequently applied to the biologically well-known network of Bacteriophage λ and ERBB Receptor-regulated G1/S transition involved in the breast cancer to successfully determine the constraints on Time Delays in less computational time. Moreover, due to its reduced complexity, the proposed approach would be able to handle BRNs with more number of genes as compared to the Hybrid Modeling approach that takes into account the complete STG while determining the constraints on Time Delays.

The rest of the paper is organized as follows: section 2 discusses the Process Hitting Framework while the basics of René Thomas Discrete Modelling Framework and Hybrid Modeling are described in section 3. The proposed methodology for incorporation of Time Delays in Process Hitting Framework is given in section 4 and the proposed Hybrid Process Hitting approach is applied to the BRNs of Bacteriophage λ and ERBB-receptor regulated G1/S transition in sections 5, 6 followed by concluding remarks and future work in sections 7, 8, respectively. A glossary of mathematical notations used in the definitions in this paper is provided as supplementary file (Data Sheet 2).

2. Process Hitting (PH) Framework

This section describes the Process Hitting Framework which is used for the static analysis of large scale biological regulatory networks (Paulevé et al., 2011; Folschette et al., 2012) and Goal-oriented model reduction.

2.1. Biological Regulatory Networks (BRNs)

BRNs are used to model the interactions between various biological entities (Proteins or RNA). Quite a few complex processes are involved in the regulations between them, however, these are simplified to only two actions, namely, activation and inhibition (Thomas and d'Ari, 1990). BRNs are often described as an Interaction Graph in which genes or other biological entities are shown as nodes and their activating or inhibiting influences on the other elements are represented by positive or negative edges, respectively. The simple representation of BRN was extended to automatically derive its behavioral model by using model-checking tools (Bernot et al., 2007). The discrete levels are defined for a particular concentration level of a gene after which its influence on other entities in the BRN changes. The activating (resp. inhibiting) influence of a gene is exerted once it is above (resp. below) the threshold level.

Definition 1. (Interaction Graph). An Interaction Graph IG = (N, E) of a Biological Regulatory Network is defined where N is the set of all nodes and E is the set of all directed edges in which an edge from node p to q is $e=\left(p\underset{}{\overset{t,s}{\to }}q\right)$ and is labelled by the pair (t, s) such that t is a positive integer and represents the qualitative threshold level required for interaction, and s ∈ {+, −} describes the type of interaction, i.e., “+” for activation and “−” for inhibition.

The interaction graph described above is quite similar with the discrete version of Deffuant model in opinion evolution in a network (Shang, 2015b).

Definition 2. (Effective Levels). For a regulation $p\underset{}{\overset{t,s}{\to }}q$ in a BRN and s being + (resp. −), the effective levels $LE{V}_p^{+}$ are those threshold levels {t; l_p} (resp.{0;t − 1}) of gene p at and above (resp. below) which it will activate (resp. inhibit) gene q, where l_p is the highest threshold level of gene p. Conversely, effective levels $LE{V}_p^{-}$ would be {0;t − 1}(resp.{t; l_p}) for which gene q would be activated by gene p.

Example 1. In the BRN shown in Figure 1B, the edge $p\underset{}{\overset{1,+}{\to }}q$ depicts the activating influence of gene p on gene q at threshold level 1. So effective levels $LE{V}_p^{+}$ are {1;1} and $LE{V}_p^{-}$ are {0;0} where gene p would be inhibiting gene q. Similarly, for some other interaction $p\underset{}{\overset{2,-}{\to }}q$ with l_p = 3, the effective levels $LE{V}_p^{+}$ are {0;1} and $LE{V}_p^{-}$ are {2;3}.

FIGURE 1

Figure 1. (A) An example of Automata Network (AN). Labeled boxes represent the Genes (automaton) whereas the circles represent the Levels (local states) within an automata and labeled with numbers, e.g., p₀ is the local state labeled 0 in the box p. Also the local transitions (e.g., $q_0\underset{}{\overset{\left\{p_1\right\}}{\to }}{q}_1$) are represented by arrows labeled by their necessary conditions. The network state 〈p₁, q₀, r₀〉 is depicted by grayed local states. (B) Interaction graph of a toy BRN.

2.2. Automata Networks

Starting from any arbitrary state of the system, Process Hitting takes into account all the possible evolutions (Generalized Dynamics) of the BRN and keeps evolving the network till it reaches a steady state where no more interactions can take place. In Process Hitting Framework, each gene is represented as automata and its current level as local state. A Process Hitting state is formed by gathering one local state of each automata present in the BRN. A local state i of a automata p is denoted as p_i. At any given time for each automaton only one local state is present which represents the current level of activity of that gene. Accordingly, the set of all local states is the global state of the BRN. The local transition is the local possible evolutions inside an automaton. The transition $p_i\underset{}{\overset{l}{\to }}{p}_j$ is the change of local state p_i to local state p_j made possible by the presence of a set l of local states belonging to other automata which should be active in the current global state. Thus, all local states in l cooperate to switch the active local state of p from p_i to p_j. It is noted that l contains at most one local state of another automaton and none from automaton p; it can also be empty. Automata Network is formally defined as:

Definition 3. (Automata Network). An Automata Network (AN) is a triple $\left(\Sigma ,L,{H}\right)$ where:

• Σ = {p, q, …} is the finite set of automata,

• L = Π_p∈ΣL_p is the finite set of global states with L_p = {p₀, …, p_lp} the finite set of local states of automata p ∈ Σ and l_p a positive integer. Also p ≠ q ⇒ ∀(p_i, q_j) ∈ L_p × L_q, p_i ≠ q_j, and LS = ∪_p∈ΣL_p denotes the set of all the local states.

• for each $p\in \Sigma ,{{H}}_p=\{p_i\underset{}{\overset{l}{\to }}{p}_j\in L_p\times \wp \left(\text{LS}\backslash L_p\right)\times L_p|p_i\ne p_j\}$, is the set of local transitions on automata p; ${H}={\cup }_{p\in \Sigma }{{H}}_p$ is the set of all local transitions in the network.

Definition 4. (Step). For an AN $\left(\Sigma ,L,{H}\right)$, a step $\theta \subseteq {H}$ is a subset of local transitions ${H}$ where, for each automaton a ∈ Σ, there is at most one local change in a.

For any local transition $\theta =p_i\underset{}{\overset{\ell }{\to }}{p}_j$, p_i is the origin, ℓ is the condition and p_j is the destination and, respectively, noted as ^•θ, cond(θ) and θ^•.

Definition 5. (Trace). For an AN $\left(\Sigma ,L,{H}\right)$ and a state l ∈ L, a trace π is a sequence of successive steps θ which leads to a defined goal state l_T. The trace π is denoted as θ₁::θ₂::….

Given an Automata Network $\left(\Sigma ,L,{H}\right)$ and a state l ∈ L, the goal state l_T ∈ L is reachable from i if there exists a trace π with ^•π ⊆ l and l_T ∈ π^•.

Example 2. Figure 1A shows the Automata Network of a toy BRN depicted in Figure 1B where:

• Σ = {p, q, r},

• L_p = {p₀, p₁}, L_q = {q₀, q₁}, L_r = {r₀, r₁} and

• the set of transitions as:

$\begin{array}{l}{H}=\{q_0\stackrel{\left\{p_1\right\}}{\to }{q}_1,q_1\stackrel{\left\{p_0\right\}}{\to }{q}_0,r_0\stackrel{\left\{p_1,q_1\right\}}{\to }{r}_1,r_1\stackrel{\left\{p_0\right\}}{\to }{r}_0,\\ r_1\stackrel{\left\{q_0\right\}}{\to }{r}_0\}.\end{array}$

Starting from the state l = 〈p₁, q₀, r₀〉, only one transition is playable: $\theta =q_0\underset{}{\overset{\left\{p_1\right\}}{\to }}{q}_1$ which updates the automaton state to l · θ = 〈p₁, q₁, r₀〉. From this state, the next possible transition is $r_0\underset{}{\overset{\left\{p_1,q_1\right\}}{\to }}{r}_1$ which leads to the state 〈p₁, q₁, r₁〉. Thus, from the state l = 〈p₁, q₀, r₀〉 to 〈p₁, q₁, r₁〉, the trace π is given by: $\pi =\{q_0\underset{}{\overset{\left\{p_1\right\}}{\to }}{q}_1\}::\{r_0\underset{}{\overset{\left\{p_1,q_1\right\}}{\to }}{r}_1\}.$

2.3. Goal Oriented Reduction and Cutsets

In the Automata Networks representing the BRNs, a goal is typically the activation (resp. inhibition) of some gene or transcription factor or kinase, etc. From an initial given state for each entity in BRN, the goal is reachable if a sequence of steps is present which leads to the state which contains the goal (Paulevé, 2018). For example, the goal in a signaling network is usually the activation of the downstream transcription factor starting from the initial inactive state.

A minimal trace corresponds to the path or sequence of steps (containing no cycle) starting from the initial state to the defined goal state containing the gene activation / inhibition. In a cycle, the initial global state is visited twice, i.e., the path leads back to its starting point. There can be multiple traces which satisfy the goal reachability in the network. So, the goal oriented reduction aims to preserve all minimal traces to the defined goal.

For a global AN $\left(\Sigma ,L,{H}\right)$, an initial state l ∈ L and a reachability goal state l_T, the goal oriented reduction identifies a subset of local transitions ${H}$ that are sufficient for producing all the minimal traces leading to l_T from l.

2.3.1. Cutsets for Goal Reachability

Cutsets are the sets of local states such that each trace reaching the goal includes a transition involving one of these local states. Also, these sets contain the necessary processes which when disabled would prevent the considered reachability goal. Hence, disabling of all transitions having precondition intersecting with cutset will remove all the traces leading to the goal. The algorithm for calculating the Cutset is described in detail in Paulevé et al. (2013) alongwith the complete methodology. An explanatory image describing the algorithm is given in Figure 2.

FIGURE 2

Figure 2. Explanatory image of the algorithm for computation of Goal-Oriented Reachability and Cut-Sets in Process Hitting (Paulevé, 2018).

The Cutsets provide information on potential therapeutic targets if the considered reachability goal represents a diseased condition by preventing all the local states of a cutset to act using other Gene Knock-out/in/down techniques (Paulevé et al., 2013).

2.4. Regulation Conditions for Cooperative Hits

The combined effect or influence of multiple genes in a BRN is governed by logical functions (Richard et al., 2006; Bernot et al., 2008) unless there are biological observations which provide evidence that the genes are exerting an influence other than that given by the logical function. In that case the biological observation take precedence over the logical function.

Here we give all the possible cases for two genes p and q to cooperate and exert influence on a third gene r at threshold levels t_p and t_q, respectively, and define the rules for logical functions for activation (upregulation) of gene r when both p and q are acting simultaneously. Figure 3 depicts the four cases which cover all possible combinations of two entities which are ¬p ∧ ¬q, ¬p ∧ q, p ∧ ¬q and p ∧ q.

FIGURE 3

Figure 3. Four possible cases of cooperation of two genes p and q to influence a third gene r. (A) Both p and q are inhibitors and would upregulate gene r when both are below threshold value t_p and t_q, respectively. (B–D) Other possible combinations of gene p and q.

Definition 6. (Regulation Condition for Individual Hits). The regulation condition Cond_ACT (resp. Cond_INH) for the individual hit contained in the set ${H}$ by a gene p having interaction with gene r as $\left(p\underset{}{\overset{t_p,s_p}{\to }}r\right)$ is defined as:

$\begin{array}{l}con{d}_{ACT}=\{\begin{array}{ll}{p}_L,\hfill & s_p=\{-\}\hfill \\ p_H,\hfill & s_p=\{+\}\hfill \end{array}\text{ and}\\ con{d}_{INH}=\{\begin{array}{ll}{p}_H,\hfill & s_p=\{-\}\hfill \\ p_L,\hfill & s_p=\{+\}\hfill \end{array}\end{array}$

Now we consider the combined influence of two genes p and q on one single entity r and define the rule for this cooperative interaction. If both the individual influencers p and q are the “inhibitors” of target gene r, then it will be activated when individually both p and q are at p_L and q_L represented by pq_LL. For all other combinations of p and q given by pq_LH, pq_HL, pq_HH, the target gene r would be inhibited.

Definition 7. (Regulation Condition for Cooperative Hits). The regulation condition Cond_ACT (resp. Cond_INH) for the cooperative hit ${H}$ by two genes p and q having their individual interactions with gene r as $\left(p\underset{}{\overset{t_p,s_p}{\to }}r\right)$ and $\left(q\underset{}{\overset{t_q,s_q}{\to }}r\right)$ respectively is defined as:

$con{d}_{ACT}=\{\begin{array}{cc}p{q}_{LL},& s_p=\{-\},s_q=\{-\}\\ p{q}_{LH},& s_p=\{-\},s_q=\{+\}\\ p{q}_{HL},& s_p=\{+\},s_q=\{-\}\\ p{q}_{HH},& s_p=\{+\},s_q=\{+\}\end{array}$

and

$con{d}_{INH}=\{\begin{array}{ll}p{q}_{LH},p{q}_{HL},p{q}_{HH},\hfill & s_p=\{-\},s_q=\{-\}\hfill \\ p{q}_{LL},p{q}_{HL},p{q}_{HH},\hfill & s_p=\{-\},s_q=\{+\}\hfill \\ p{q}_{LL},p{q}_{LH},p{q}_{HH},\hfill & s_p=\{+\},s_q=\{-\}\hfill \\ p{q}_{LL},p{q}_{LH},p{q}_{HL},\hfill & s_p=\{+\},s_q=\{+\}\hfill \end{array}$

We see that for Case I (¬p ∧ ¬q) shown in Figure 3A, the logical AND function implies that q would be upregulated when both p and q are below threshold t_p and t_q, respectively, which is represented by pq_LL and it would be downregulated for all other values, i.e., pq_LH, pq_HL, and pq_HH. Once the appropriate regulation condition has been obtained then the corresponding effective levels are determined in a straight forward manner. For Cond_ACT = pq_LL the effective levels are found to be $LE{V}_p^{-}=\left\{0;t_p-1\right\}$ and $LE{V}_q^{-}=\left\{0;t_q-1\right\}$. Similarly, the conditions and effective levels for activation (resp. inhibition) for Cases II, III and IV are given in Table 1.

TABLE 1

Table 1. The regulation conditions and effective levels for Up (resp. Down) Regulation by two genes p and q having a combined influence on another gene in a BRN are governed by the logical function.

2.4.1. Multi-Level Genes

Now we consider the interactions involving multi-level genes which form a cooperative sort. For instance, consider gene p and gene q to have 3 and 4 threshold levels, respectively (i.e., l_p = 1 and l_q = 2), and they interact with gene r at threshold level t_p = 1 and t_q = 2, respectively, as shown in Figure 4A. Since the individual genes activate (resp. inhibit) the target gene for all the discrete levels above (resp. below) the threshold level, therefore, this interaction can be represented by a Boolean function with H containing all the discrete levels required for activation and L having levels for inhibition. Then for the case ¬p ∧ q which is the logical function for upregulation of gene r, the corresponding Cond_ACT is pq_LH which gives us the effective levels $LE{V}_p^{-}<1=\left\{0;0\right\}$ and $LE{V}_q^{+}\ge 2=\left\{2;2\right\}$. Figure 4B shows the AN for gene r on which multi-level genes p and q are acting at the same time. The corresponding condition for upregulation for gene r requires p_L = {p < 1} and q_H = {q ≥ 2}. Similarly, the downregulation takes place when either p_H = {p ≥ 1} or q_L = {q < 2}.

FIGURE 4

Figure 4. (A) Multi-level genes p and q act on gene r. For the case ¬p ∧ q and gene p having l_p = 1 and q having l_q = 2 interacting with gene r at threshold levels (1, −) and (2, +), respectively. (B) The Automata Network for gene r on which multi-level genes p and q are acting at the same time. The upregulation of AN for gene r takes place for the condition p_L = {p < 1} and q_H = {q ≥ 2}.

2.5. Abstraction of Switch Conditions for Gene Regulation

We intend to model the Process Hitting in a BRN with the help of concurrent finite timed automata (Alur, 1999) in which the qualitative levels of each gene would be represented by the states of its corresponding automaton and the gene evolution would be determined by the automaton dynamics. For this we need to abstract the switch conditions under which the automaton would move from one state to another.

Our approach is to directly derive the switch conditions from the effective levels corresponding to each action contained in the set ${H}$. In order to model the timed automaton representing the dynamics of gene r, we abstract the switch conditions which are required for it to move from one state to another. These switch conditions are abstracted directly from the Regulation conditions and effective levels.

Definition 8. (Hit Part). We abstract the hit part μ from the permissible process hits contained in ${H}$ corresponding to up (resp. down) regulation of each gene in the following manner:

$\mu =\{\begin{array}{ll}{p}_i,\hfill & i\in \{L,H\},forhitbysinglegene\hfill \\ p_{xy},\hfill & xy\in \{LL,LH,HL,HH\},forcooperativehit\hfill \end{array}$

where L = [0; t_p − 1] and H = [t_p; l_p].

Corresponding to each hit part μ, we obtain the qualitative threshold levels which permit the change in expression level of the target gene. From these levels we deduce the appropriate switch conditions ∧ which result in the dynamics of the target gene.

Definition 9. (Switch Conditions for Gene Regulation). We define Switch Conditions ∧ for each process hit contained in the set ${H}$ as below:

1. for every hit by a single gene μ = p_i, i ∈ {L, H}, the switch condition is:
$\wedge =\{\begin{array}{ll}p=\left\{0,\dots ,t_p-1\right\},\hfill & \text{ for }i=L\hfill \\ p=\left\{t_p,\dots ,l_p\right\},\hfill & \text{ for }i=H\hfill \end{array}$

2. for every cooperative hit μ = p_xy, xy ∈ {LL, LH, HL, HH}, the switch condition is:
$\wedge =\{\begin{array}{l}p=\left\{0,\mathrm{...},t_p-1\right\},q=\left\{0,\mathrm{...},t_q-1\right\},\hfill \\ \text{ for }x=L,y=L\hfill \\ p=\left\{0,\mathrm{...},t_p-1\right\},q=\left\{t_q,\mathrm{...},l_q\right\},\hfill \\ \text{ for }x=L,y=H\hfill \\ p=\left\{t_p,\mathrm{...},l_p\right\},q=\left\{0,\mathrm{...},t_q-1\right\},\hfill \\ \text{ for }x=H,y=L\hfill \\ p=\left\{t_p,\mathrm{...},l_p\right\},q=\left\{t_q,\mathrm{...},l_q\right\},\hfill \\ \text{ for }x=H,y=H\hfill \end{array}$

Example 3. For the transition $r_0\underset{}{\overset{p_0,q_1}{\to }}{r}_1$ given in Example 3 in which the gene r is upregulated from r₀ to r₁ in the presence of p₀ and q₁, the effective levels are obtained as $LE{V}_p^{-}=0$ and $LE{V}_q^{+}=1$. We see that gene r can change its state from r₀ to r₁ only when gene p and q are at threshold levels 0 and 1, respectively. Therefore, the appropriate switch condition is derived as ∧ :{p = 0, q = 1}. For the multi-level genes shown in Figure 4B, the effective levels for regulation condition Cond_ACT of gene r were obtained as $LE{V}_p^{-}<1=\left\{0;0\right\}$ and $LE{V}_q^{+}\ge 2=\left\{2;3\right\}$. Thus, in this case the switch conditions would be ∧ :{(p = 0, q = 2) ∨ (p = 0, q = 3)} which are all the combinations of effective levels of genes p and q at which gene r is upregulated.

3. Hybrid Modeling

3.1. René Thomas' Logic Formalism

Initially, René Thomas presented a qualitative formalism based on Boolean logic applicable on Biological Regulatory Networks (Thomas, 1978, 2013) which closely approximated the ODE models. Later, the boolean model was found unable to represent various interactions taking place in BRNs at varying gene expression concentrations. This led to the presentation of kinetic logic formalism which allows the modeling of discretely abstracted concentration levels other than boolean as well (Thomas and d'Ari, 1990). It has been further enriched with Parametric Time Delays in Hybrid Modeling (Ahmad et al., 2006).

3.2. Hybrid Modeling Framework

It is difficult to model the changes in concentration levels of a protein through the discrete modeling frameworks. For instance, the protein expression level which often follows a sigmoidal curve (Figure 5A) is represented by positive integers in a discrete model where the change in level from 0 to 1 takes place instantly (Figure 5B) which is not a true representation of the actual increase in concentration of protein in a cellular environment. Therefore, the changes in concentration expression levels were proposed to be represented by piecewise linear curves for the modeling of the sigmoidal nature of protein expression (Ahmad et al., 2008) as shown in Figure 5C, as opposed to the step functions used in discrete models. Figures 5D–F) depict the modeling in case of degradation of protein.

FIGURE 5

Figure 5. Activation and Inhibition delays. The clock variable h_p measures the time of evolution between two discrete levels. Initially the clock is set to zero and the changes in the level occurs in a delay time $d_p^{+/-}$. The change in expression level during Gene Activation is represented by (A) Sigmoidal function (B) Step function and (C) Piece-wise Linear function. Similarly, the change in expression level during Gene Inhibition is represented by (D) Sigmoidal function (E) Step function and (F) Piece-wise Linear function.

A clock variable h_p is associated with each protein to measure the time it takes (delay) to go from a particular expression level to the next one. Initially, the value of each delay h_p is set to zero, which then increase to $d_p^{+}$ (resp. $d_p^{-}$) that signifies the production (resp. degradation) delay, i.e., the time it takes to increase (resp. decrease) the concentration level of the particular protein by 1 threshold level. Evolution rate for each clock h_p is set using first order derivative dh_p/dt = x where x ∈ {0, 1, −1} which represents no change, increase or decrease, respectively (Ahmad et al., 2008). Since the exact value of the delays associated with the production or degradation of proteins is not known in most cases, therefore, unvalued parametric delays are used instead and are determined through hybrid modeling (Ahmad et al., 2012). Hybrid modeling approach has been applied for behavioral modeling of several biological networks including MAL-associated network (Ahmad et al., 2012), dengue virus pathogenesis and clearance mechanism (Aslam et al., 2014), tail-resorption in tadpole's metamorphosis (Khalis et al., 2009), role of O-linked N-acetylglucosamine transferase in oncogenesis and cancer progression in hexosamine biosynthetic pathway (Saeed et al., 2016) and immunity control mechanism in bacteriophage lambda (Richard et al., 2006).

The hybrid modeling method has also been applied in opinion evolution in social networks (Shang, 2018). The Hybrid consensus for averager-copier voter networks with non-rational agents has been modeled and it is shown that there is a distinct influence of voters on the consensus and evolution of the opinion process.

The delays are of two types; Signal Transmission Delay and Signal Processing Delay. While the outcome for both types of delay would be the same, however, the processing delay is likely to cause errors in the consensus behavior whereas it is independent of the transmission delay (Shang, 2017). For study of BRNs we have to deal with the processing delay, hence, keeping track of the transient changes during the delay period is important and essentially required to be taken into account for correctly modeling the evolution of BRN and understanding its dynamical behavior.

3.3. Parametric Bio Linear Hybrid Automaton

We formally define hybrid automaton here by mainly using the notations and definitions given in Alur (1999) and adapted from Ahmad et al. (2006). We use a set of clock variables X = {h₁, …, h_n} to represent time which progress synchronously in accordance with $\dot{h_i}=1$ or −1 and can be reset to zero. The clock constraints φ for the set ϕ(X) is defined by the grammar:

$\phi ::=\text{true}|h\le d|h\ge d|h=d|{\phi }_1\wedge {\phi }_2$

where h ∈ X and d is an un-valued parameter.

Definition 10. (Parametric Biological Linear Hybrid Automaton (Bio-LHA)). A parametric Bio linear hybrid automaton 𝔹 is a tuple (L, l₀, D, X, E, Inv, Dif) where,

• L is a finite set of locations,

• l₀ ∈ L is the initial location,

• D is a finite set of parameters (delays),

• X is a finite set of real-valued variable (clocks),

• E ⊆ L × φ × 2^X × L is a finite set of edges with typical element e = (l, g, R, l′) ∈ E representing an edge from l to l′ with guard g (clock constraint of the form h = d) and the reset set R ⊆ X,

• Inv: L → φ assigns an invariant to any location,

• Dif: L × X → {−1, 0, 1} maps each pair (l, h) to an evolution rate.

The Transition System related semantics of the parametric Bio-LHA is given below according to the time domain 𝕋, where 𝕋* = 𝕋\{0}.

Definition 11. (Semantics of Bio-LHA). Let γ be a valuation for the parameters P and ν represents the values of clocks in a location. The (𝕋, γ)-semantics of a parametric Bio-LHA 𝔹 = (L, l₀, P, X, E, inv, dif) is defined as a timed transition system _B = (𝕊, s₀, 𝕋, →) where:

(1) 𝕊 = {(l, ν)|l ∈ L and ν ⊧ Inv(l)};
(2) s₀ is the initial state and
(3) the relation → ⊆ 𝕊 × 𝕋 × 𝕊 is defined for t ∈ 𝕋 as:

• discrete transitions: $\left(\ell ,\nu \right)\underset{}{\overset{0}{\to }}\left({\ell }^{\prime },{\nu }^{\prime }\right)$ iff ∃(ℓ, g, R, ℓ′) ∈ E such that g(ν) = true, ν′(h) = 0 if h ∈ R and ν′(h) = ν(h) if h ∉ R.

• continuous transitions: For t ∈ T*, $\left(\ell ,\nu \right)\underset{}{\overset{t}{\to }}\left({\ell }^{\prime },{\nu }^{\prime }\right)$ iff ℓ′ = ℓ, ν′(h) = ν(h) + Dif(ℓ, h) × t, and for every t′ ∈ [0, t], (ν(h) + Dif(ℓ, h) × t′⊧Inv(l), where ⊧ represents satisfaction operator.

3.4. Network of Bio-LHA

A BRN is composed of various entities which evolve concurrently depending on the activating (resp. inhibiting) influence of other entities as defined in the BRN. We have represented each entity with a Bio-LHA and its dynamics is also defined above. Next we define the network of Hybrid Automaton 𝔹_N which is built up from composition of individual Bio-LHA 𝔹₁, …, 𝔹_n.

Definition 12. (Network of Bio-LHA). Let 𝔹_i: = (L_i, l_{0, i}, D_i, X_i, E_i, Inv_i, Dif_i), i ∈ 1, …, n be the component Bio-LHA representing the genes of the BRN. The network of Bio-LHA representing n genes is 𝔸 = (L, l₀, D, X, E, Inv, Dif) where,

• L = L₁×… × L_n, is the set of locations,

• l₀ = (l_{0, 1}…, l_{0, n}) ∈ L is the initial location,

• $D={\displaystyle \underset{i=1}{\overset{n}{\cup }}{D}_i}$, is the set of delays,

• $X={\displaystyle \underset{i=1}{\overset{n}{\cup }}{X}_i}$, is the set of clocks,

• The transition relation E is defined by the following asynchronous rule:

   If l_i is any component location of l ∈ L then $\left(\left(\dots l_i,\dots \right),g_i,R_i,\left(\dots l_i^{\prime },\dots \right)\right)$ ∈ E for $\left(l_i,g_i,R_i,l_i^{\prime }\right)\in E_i$ and all other component locations of l do not change. The guard g_i is the conjuction of switch condition (∧) and clock constraint (φ),

• Inv : $L\to \underset{i=1}{\overset{n}{{\displaystyle \cup }}}$, Inv_i(l), is mapping of invariants to locations,

• Dif : L×X → {−1, 0, 1}, is the mapping of clocks to evolution rates.

3.5. Representation of Hybrid Automaton

The modeling of Hybrid Automaton representing evolution of gene q is depicted in Figure 6 which is inhibited by gene r. The set of transitions ${H}$ is given by $\{q_0\underset{}{\overset{r_0}{\to }}{q}_1\},\{q_1\underset{}{\overset{r_1}{\to }}{q}_0\}$. The automaton is modeled with two states q₀ and q₁ which represents the discrete threshold levels of the gene q. Each state is labeled with its evolution rate $d_{q_1}^{+/-}$ and its invariant for which the automaton remains in that state. The edges are labeled with the guard condition on the clock variable as well as the switch condition ∧ on the gene level $LE{V}_q^{+/-}$. When the guard becomes True, the automaton moves to the next state with the clock being reset and also updating the variable q which contains the current value of gene threshold level.

FIGURE 6

Figure 6. A toy example of hybrid automaton. Various labels describe the function of each symbol used in the automaton.

3.6. Implementation of Hybrid Automaton With HyTech

HyTech (Henzinger et al., 1997) is a model checking tool developed for the formal verification of real-time systems. It is suitable for modeling of Linear Hybrid Automata and has a powerful set of analysis commands as well as parametric synthesis capability. Using HyTech for modeling of BRN, each gene is modeled by an automaton. Thus, the current state of the BRN would be given by the current state of all the automata considered together.

Example 4. The code for implementation of toy example of hybrid automaton geneQ shown in Figure 6 in HyTech is given below. Keyword loc is used to specify the state of automaton (e.g., q₀), keyword while is followed by the invariant condition (e.g., hq <= dpq1), keyword wait is used to give the rate of increase or decrease during the transition (e.g., {dhq = 1}), keyword when is followed by the guard condition (e.g., r = 0 & hq = dpq1), keyword do is used for updating of the state variable (e.g., {hq′ = 0, q′ = 1}) and keyword goto is followed by the next state of the automaton (e.g., q₁).

automaton geneQ

synclabs: ;

initially q0 & q=0 & hq=0;

loc q0: while hq < =dpq1 wait {dhq=1}

when r=0 & hq=dpq1 do {hq'=0,q'=1} goto q1;

loc q1: while hq < =dnq1 wait {dhq=-1}

when r=1 & hq=dnq1 do {hq'=0,q'=0} goto q0;

end --geneQ

We use the parametric delay variable hq here with all the transitions between different states of automaton which evolve at some rate till the switch condition ∧ becomes true for respective transitions. In this way we are able to get the range of values for this delay variable for the automaton to move from one state to another. Accordingly, we are able to determine the constraints on the delay variable corresponding to each state of automaton. When the system of automata representing the complete interaction graph of BRN is composed with switch conditions and delay variables for each gene evolution, we get the constraints on delay variables for each of the composed state of automata. The complete methodology is given in next section.

4. Incorporating Time Delays in Process Hitting Framework

The aim of this section is to give the methodology of our proposed Hybrid Process Hitting technique in which we have incorporated the Time Delays from Hybrid Modeling described above into the Process Hitting Framework.

4.1. Obtaining Process Hits for BRN

Starting with the Interaction Graph of a BRN which specifies all the regulations (type and threshold level) of all its genes on each other, all the local transitions in the BRN were captured through application of Process Hitting as a first step. We now consider an example BRN in which two individual genes p and q are interacting with gene r as shown in Figure 7A. By considering the permissible actions in the generalized dynamics in this BRN we get the possible local transitions ${H}$ as given below:

${H}=\{q_0\stackrel{\left\{r_0\right\}}{\to }{q}_1,q_1\stackrel{\left\{r_1\right\}}{\to }{q}_0,r_0\stackrel{\left\{p_1,q_1\right\}}{\to }{r}_1,r_1\stackrel{\left\{p_0\right\},\left\{q_0\right\}}{\to }{r}_0\}.$

FIGURE 7

Figure 7. (A) An example BRN in which three genes p, q, and r interact with each other. Gene p and q are the activators of gene r whereas gene r inhibits gene q. (B) State Transition Graph for the toy BRN is shown in which the system of automata comprising of genes p, q, and r completes a cycle from state 〈1, 0, 0〉 through states 〈1, 1, 0〉, 〈1, 1, 1〉, and 〈1, 0, 1〉.

4.2. Deducing the Switch Conditions and Effective Levels

Now using Table 1, the effective levels of interacting genes p and q are determined from the corresponding condition required by the cooperative sort for upregulation of gene r. For pq₁₁ we get $LE{V}_p^{+}=\left\{1\right\}$ for gene p and $LE{V}_q^{+}=\left\{1\right\}$ for gene q. We interpret the values of effective levels in a straight forward manner in which the levels represent the current state of the individual automaton of gene p and q. In this case, the automaton geneP is at state P₁ and automaton geneQ is at state Q₁. For automaton geneR which describes the dynamics of gene r, the effective levels, therefore, specifies the guard condition in terms of exact combination of state variables of other automata that would induce the transition from state R₀ to R₁.

4.3. State Transition Graph

All the automata when considered together would represent the complete model of BRN. Hence, the transitions from one state to another would be governed by the dynamics of the modeled BRN. At any given time, the locations of all the automata give us the current state of the BRN. A change in location of any one automata result in change of state of BRN. For the BRN shown in Figure 7A consisting of genes p, q and r with expression levels as 1, 0 and), respectively, the current state is denoted as 〈1, 0, 0〉 with corresponding locations of automata as loc p₁, loc q₀ and loc r₀. Here, we note that the switch condition for evolution of automaton geneQ is satisfied (i.e., $LE{V}_r^{+}=\left\{0\right\}$) for it to move from location loc q₀ to loc q₁, thus the next state of the system would be 〈1, 1, 0〉.

In the state 〈1, 0, 0〉, it may be noted that the required switch conditions for evolution of automaton geneR were not satisfied, so this automaton could not evolve. However, we see that in state 〈1, 1, 0〉 the switch conditions for evolution of geneR are now satisfied (i.e., p = 1 and q = 1) therefore, the next transition lead the system to state 〈1, 1, 1〉. Subsequently, the system evolves to state 〈1, 0, 1〉 from where it comes back to state 〈1, 0, 0〉. The corresponding state transition graph for the dynamics of the system model is shown in Figure 7B.

4.4. Enriching the Automaton With Time Delays

The transition from one state of the system to another takes place when the concentration level increases (resp. decreases) from one threshold level to another and it takes a finite amount of time. This time is measured in terms of delay variable $d_q^{+/-}$ in hybrid modeling of BRN. We note that the automaton moves from one state to another instantaneously once the switch conditions are satisfied. So in order to measure the time delay for this transition we introduce the parametric clock variable hq which increases at a rate dhq/dt for the time the automaton remains in a particular state and its value is recorded on transition to the next state.

Consider the initial location of automaton geneQ is loc q₀ and the clock variable h_q is initially set to 0 which starts to increase with rate dhq/dt. On satisfaction of switch conditions, automaton geneQ leaves its current state loc q₀ and moves to next state loc q₁ which represents the transition of Gene q from expression level 0 to 1. At this time the value of clock variable h_q for this transition is noted in terms of delay variable as $d_{q_1}^{+}$. Similarly, the delay for moving from level 1 to 0 will be given by $d_{q_1}^{-}$.

4.5. Composing System of Automata

We compose the System of Automata by considering all individual automaton together. The transitions in each automaton are linked with the state of other automata through switch conditions and the time taken by each transition is obtained in terms of its delay variable. We recall that the switch conditions are given by the refined cooperative hits and represent the permissible dynamics in the BRN. When all the individual automaton are considered together we obtain the constraints on the range of values of respective clocks for the system to move from one state to another.

Example 5. In order to describe our methodology for obtaining the time delays, we now model the individual genes in BRN of Figure 7A with respective automaton. From ${H}$ we note the action (switch condition) for upregulation q₀ to q₁ for automaton geneQ as r₀ and effective level $LE{V}_r^{+}=\left\{0\right\}$. Therefore, the guard condition for automaton geneQ is set as r = 0 & $h_q=d_{q_1}^{+}$. Similarly for downregulation, the guard condition is r = 1 & $h_q=d_{q_1}^{-}$. In case of geneR, we have the switch condition for upregulation as pq₁₁ given by the logical function p ∧ q with effective levels as $LE{V}_p^{+}=\left\{1\right\}$ and $LE{V}_q^{+}=\left\{1\right\}$. For downregulation, we get the switch conditions and effective levels from Table 1 for the given logical function. Accordingly, the guard conditions in automaton geneR are set to model this behavior. The individual automaton for p, q and r are modeled as shown in Figure 8. The corresponding hytech file for modeling this BRN is given in Supplementary File 1. We assume that gene p is at level 1 initially and hytech file is run to get the constraints on delay parameters corresponding to each state of composed system as shown in Table 2 below.

FIGURE 8

Figure 8. The individual automaton modeling the behavior of genes p, q, and r. Each automata contains the respective guard and switch conditions.

TABLE 2

Table 2. Delay Constraints obtained for the example BRN through modeling of its corresponding system of Automata in HyTech.

4.6. Inference of Parametric Delays in Hybrid Process Hitting

The range of delay parameters alongwith the constraints on delay variables in terms of time delays is obtained as output from simulation of System of Automata which modeled the proposed Hybrid Process Hitting approach. As evident from Table 2, for each path or trajectory from one particular state to another of the State Transition Graph, we get the constraints on clock variables in terms of time delays for each gene of BRN. For instance, the delay constraint for the path from state p₀q₀r₀ to state p₁q₀r₀ requires that $d_{p_1}^{+}\le d_{q_1}^{+}$ and $d_{p_1}^{+}\le d_{r_1}^{+}$. Similarly, the delay constraints corresponding to the paths between various states are obtained from running the HyTech model.

It is highlighted that the delays $d_{q_1}^{+}$ and $d_{q_1}^{-}$ represents the accumulation and degradation time of gene q. Hence, the constraints on the delays give us the conditions in terms of time taken for accumulation or degradation of biological entities in the BRN and clearly define the dynamics of the modeled system.

We are thus able to determine the dynamical behavior of a biological regulatory network in which the constraints on time delays are deduced for each state.

5. Biological Application I: Bacteriophage Lambda

In this section we apply the proposed Hybrid Process Hitting approach to the well-known example of Bacteriophage λ (Thieffry and Thomas, 1995). This network has been widely studied because of its peculiar switch mechanism because of which this virus chooses between lysis and lysogenization. These two responses are the result of different dynamics of the BRN and thus leads to separate attractors.

5.1. Automata Network

We take into consideration the BRN model proposed in Thieffry and Thomas (1995) and shown in Figure 9. The network has 4 nodes, having expression levels 3, 4, 2, and 2, respectively, so in this case:

Σ = {cI, cro, cII, N}, L_cI = {cI₀, cI₁, cI₂},
L_cro = {cro₀, cro₁, cro₂, cro₃}, L_cII = {cII₀, cII₁} and
L_N = {N₀, N₁}.

FIGURE 9

Figure 9. BRN of Bacteriophage λ immunity response showing interactions between nodes cI, cro, cII, and N.

We determine the local transitions ${H}$ for this network from the BRN interaction graph. The associated Automata Network is shown in Figure 10.

FIGURE 10

Figure 10. The Automata Network model of Bacteriophage λ immunity response. The corresponding local transitions ${H}$ are: l₁ = l₂ = l₃ = {cI_{≤ 1}}, l₄ = {cI₂, cro₃}, l₅ = l₆ = {cI₂}l₇ = l₈ = {cII₁, cro₀}, l₉ = l₁₀ = {cII₀}, {cro _{≥ 1}}l₁₁ = {N₁, cI_{≤ 1}, cro_{≤ 2}}, l₁₂ = {N₀}, {cI₂}, {cro₃}l₁₃ = {cI _{< 1}, cro _{< 2}}, l₁₄ = {cI _{≥ 1}}, {cro _{≥ 2}}.

5.2. Obtaining the Switch Conditions

The logical function for the cooperation between cI and cro is ¬cI ∧ ¬cro as both nodes are inhibitors of node N which implies that the process cI − cro_LL produce upregulation of node N while all other processes, i.e., cI − cro_LH, cI − cro_HL, and cI − cro_HH would result in downregulation of target node N. Here we note that the nodes cI and cro are having multiple expression levels but the threshold level for interactions on node N are 1 and 2, respectively. This condition then gives us the effective levels for the process cI − cro_LL which are $LE{V}_{cI}^{-}=\left\{0\right\}$ and $LE{V}_{cro}^{-}=\left\{0,1\right\}$. Thus, the gene N would bounce from N₀ to N₁ if its regulators have the expression levels (cI₀, cro₀) or (cI₀, cro₁).

Similarly, we get the effective levels and corresponding expression levels for downregulation of node N from the other processes of cooperative sort cI − cro as follows:

$cI-cr{o}_{LH}:LE{V}_{cI}^{-}=\left\{0\right\},LE{V}_{cro}^{+}=\left\{2,3\right\}$
(cI₀, cro₂), (cI₀, cro₃)
$cI-cr{o}_{HL}:LE{V}_{cI}^{+}=\left\{1,2\right\},LE{V}_{cro}^{-}=\left\{0,1\right\}(c{I}_1,$
$cr{o}_0),\left(c{I}_1,cr{o}_1\right),\left(c{I}_2,cr{o}_0\right),\left(c{I}_2,cr{o}_1\right)$
$cI-cr{o}_{HH}:LE{V}_{cI}^{+}=\left\{1,2\right\},LE{V}_{cro}^{+}=\left\{2,3\right\}$
$\left(c{I}_1,cr{o}_2\right),\left(c{I}_1,cr{o}_3\right),\left(c{I}_2,cr{o}_2\right),\left(c{I}_2,cr{o}_3\right)$

5.3. Automaton Specification

These expression levels give us the required switch conditions for completely specifying the automaton for node N. It contains two states, namely, loc N₀ and loc N₁ since node N has two processes and are represented by state variable N = 0 and N = 1. The clock variable is h_N and the parametric delay variables are $d_{N_1}^{+}$ (resp. $d_{N_1}^{-}$) for the time taken for up (resp. down) regulation for node N. By following our proposed methodology we deduce the Invariant and guard conditions for each state of the automaton to get the Hybrid Process Hitting model for node N as shown in Figure 11.

FIGURE 11

Figure 11. Automaton for node N. The automaton is constructed to represent the dynamics in node N due to its activators (resp. inhibitors) which cooperate under the logical function ¬cI ∧ ¬cro and corresponding switch conditions. The automaton is also enriched with parametric clock variable dN which gives us the measure of time taken for the automaton to change states in terms of delays $d_{N_1}^{+}$ (resp. $d_{N_1}^{-}$).

Following the same methodology, the automaton for other nodes cI, cro and cII are constructed with their switch conditions obtained from the logical functions ¬cI ∧ ¬cro ∧ cII, ¬cI ∧ ¬cro and ¬cI ∧ ¬cro ∧ N, respectively. These logical functions are the result of refined cooperative hits obtained from Process Hitting. The cooperativity among three nodes as required for nodes cI and cII is constructed in two steps as explained earlier in section 2 above. The specifications of these automata are completed by introducing the respective parametric clock and delay variables. The automata for these other entities, cI, cII and cro are shown in Figure 12.

FIGURE 12

Figure 12. The individual automaton modeling the behavior of genes cI, cII, and cro. Each automata contains the respective guard and switch conditions.

5.4. Modeling of Clock Variables for Multi-Level Genes

In the Bacteriophage λ BRN we have two multi-level genes, i.e., cI and cro which have 3 and 4 levels, respectively. For modeling these multi-level genes in HyTech, we introduce new clock variable h_n. Consider the Automaton of gene cI which starts from initial location cI₀ and moves to next location cI₁ when the appropriate switch conditions (cII = 1, cro < 1) are met. Here if the switch conditions remain the same, the automaton would move to location cI₂, however if either of the switch condition changes then the automaton would move towards location cI₀. We note that the jump from cI₁ to cI₂ represents the accumulation of cI whereas moving from cI₁ to cI₀ shows the degradation of entity cI. Hence, we require two clock variables h_cI and h_ncI to keep track of the accumulation and degradation delays, respectively.

5.5. Composition of System of Automata

Once all the individual automaton are constructed, we now compose the system of automata by executing the HyTech file. Since we are interested in the dynamical characteristics of the transitions that take place, therefore, we use the reach command in the Analysis section of HyTech. So starting from any arbitrary initial state we are able to find the next states that the system of automata can reach alongwith the corresponding constraints on the time taken for this change of state in terms of time delays.

5.6. Results

We compose the system of automata for Bacteriophage λ and get the results on the constraints of time delays corresponding to each state that the system can reach as given in Supplementary File 2. It is highlighted that we get 42 states instead of the total 48 possible states which shows that the state transition graph has been reduced through application of Process Hitting. Since we used the reach forward command, therefore, the HyTech has computed all trajectories alongwith the delay constraints for the system of automata which it can possibly traverse with the defined switch conditions for each automaton.

We know that the phage λ BRN has two attractors, i.e., lysogenic corresponding to state 〈cI, cro, cII, N〉 = 〈2, 0, 0, 0〉 and lysis alternating between states 〈cI, cro, cII, N〉 = 〈0, 2, 0, 0〉 and 〈0, 3, 0, 0〉 (Thieffry and Thomas, 1995). From the results, we note the constraints on time delays for the paths leading to above mentioned states. For lysogenic condition 〈2, 0, 0, 0〉, gene cI reaches its highest threshold level of 2 and represses the expression of all other genes. The time delays for the path leading to this state are found to be:

$d_{c{I}_2}^{+}\le \left|d_{c{I}_1}^{-}\right|\wedge d_{c{I}_1}^{+}+d_{c{I}_2}^{+}\le d_{cr{o}_1}^{+}\wedge d_{c{I}_2}^{+}\le d_{cI{I}_1}^{+}\wedge d_{c{I}_1}^{+}\le d_{cI{I}_1}^{+}$

In lysis, the delays for the path leading to state 〈0, 2, 0, 0〉 are:

$\begin{array}{l}{d}_{cr{o}_2}^{+}\le \left|d_{cr{o}_1}^{-}\right|\wedge d_{cr{o}_2}^{+}\le d_{cI{I}_1}^{+}\wedge d_{cr{o}_1}^{+}\le d_{cI{I}_1}^{+}\wedge d_{cr{o}_1}^{+}+d_{cr{o}_2}^{+}\le \\ d_{N_1}^{+}.\end{array}$

and for the path leading to state 〈0, 3, 0, 0〉 the delays are:

$d_{cr{o}_2}^{+}\le \left|d_{cr{o}_1}^{-}\right|\wedge d_{cr{o}_3}^{+}\le \left|d_{cr{o}_2}^{-}\right|\wedge d_{cr{o}_2}^{+}\le d_{cI{I}_1}^{+}\wedge d_{cr{o}_1}^{+}\le d_{cI{I}_1}^{+}\wedge$

$d_{cr{o}_1}^{+}+d_{cr{o}_2}^{+}+d_{cr{o}_3}^{+}\le d_{N_1}^{+}\wedge d_{cr{o}_3}^{+}\le d_{cI{I}_1}^{+}\wedge$

$2d_{cr{o}_1}^{+}+2d_{cr{o}_2}^{+}+2d_{cr{o}_3}^{+}\le d_{c{I}_1}^{+}+2d_{cI{I}_1}^{+}\wedge d_{cr{o}_1}^{+}+d_{cr{o}_2}^{+}+d_{cr{o}_3}^{+}\le$
$2d_{cI{I}_1}^{+}$.

By analysing the delay constraints for each state we get the relationship between the rates of accumulation or decay of a particular entity in the BRN. These results on delay constraints are found to be slightly different from Ahmad et al. (2008) due to the difference in the way the delays are modeled. In Ahmad et al. (2008), the time delays between various levels of af a multi-level gene were modeled with the same delay variable whereas we have specified it to be different for each change in level. For instance, our modeling methodology considers $d_{cr{o}_1}^{+}$ as the time delay for cro to change from level 0 to 1 while $d_{cr{o}_2}^{+}$ gives the time delay from level 1 to 2. But in Ahmad et al. (2008), both these delays were modeled with the delay variable $d_{cro}^{+}$. Our specification of time delays is considered more suitable for modeling of dynamic behavior of genes in BRNs since the time for change amongst various levels of a biological entity are usually different from each other.

6. Biological Application II: ERBB Receptor-Regulated G1/S Transition BRN

Here we consider a comparatively large BRN of ERBB Receptor-regulated G1/S transition involved in the breast cancer and has regulations between 20 genes (Sahin et al., 2009). The BRN for this network is shown in Figure 13. The gene pRB is activated by this activation cascade network which subsequently controls the G1/S transition involved in cell divisions. The input of this network is the gene EGF which when expressed will result in potential activation of gene pRB.

FIGURE 13

Figure 13. ERBB receptor-regulated G1/S transition BRN reproduced from Sahin et al. (2009).

In multicellular organisms, growth factors play a significant role in order to maintain proper growth, development and homeostasis. In regard to receptor tyrosine kinases (RTKs), the epidermal growth factor i.e., EGF family of RTKs also known as ERBB receptors is considered as hallmark of the human cancers along with process of growth, development and physiology. The ERBB family of protein functions through homodimers and hetrodimers and thus, regulates the G1/S transition during cell cycle progression in eukaryotic organisms by modulating the activity of the Cyclin D, Cyclin E/CDK complex, the c-MYC oncogene, and the p27 kinase inhibitor.

In order to maintain a continuous simulation of EGF, the initial state of ER − α node was set to 1. Considered to be the inhibitors of CDKs, p21 and p27 were assigned the initial state of 0. This is because the expression level of both of these inhibitors was found to be high during G0/G1 arrested cells and once cells progress through S phase, their level gets decreased because of their proteasomal degradation. Similarly, the initial state of cyclin D1 was set to 0. Now, using the above mentioned initial states, all possible state transitions were evaluated until a unique stable state was obtained that had fixed levels of all the proteins as shown in Table 3.

TABLE 3

Table 3. Stable states for ERBB receptor-regulated G1/S transition obtained by PINT software (Paulevé, 2017).

We apply the proposed Hybrid Process Hitting approach to this network to obtain the Hybrid Automaton for determining the time delays involved in its dynamical evolution. First we note the switch conditions and by using them we model the Hybrid Automaton. Once this Hybrid Automaton model was run in HyTech to compute the time delays, it was observed that the software went Out of Memory. This clearly showed the limitation of HyTech as it was unable to handle 20 Hybrid Automata at one time. The problem was overcome by utilizing the powerful reduction feature of PINT software (Paulevé, 2017) which is able to provide Goal-oriented model reduction through its Cutsets feature.

Typically a goal is the activation / inhibition of some gene in the modeled Automata Networks for BRNs. From a given initial state, the goal is reachable if a sequence of steps exists which leads to a the state in which goal is present. The sequence of states leading to the goal are known as trace. So Cutsets are sets of those local states which are included in traces leading to the goal state. Thus, the goal oriented reduction would be to apply Cutsets so as to remove all traces leading to the goal. These Cutsets are implemented by the knock-out of the intermediate genes present in the network.

In the ERBB network, Normal epithelial cell lines exhibit high expression levels of ERBB1 along with ERBB2 but, have low levels of ERBB3. Unlike ERBB1, ERBB2 does not show any ligand binding activity and ERBB3 tends to have defective tyrosine kinase domain. Due to which both ERBB2 and ERBB3 are unable to activate the key signaling intermediates MEK1 and AKT1. Hence, we began our activation of network through ER − α and eliminated all the ERBB members including, ERBB1, ERBB1-2, ERBB1-3, ERBB2, ERBB2-3, ERBB3 by applying cut-sets in PINT tool using the following command:

y = y.having(ERalpha=1) .reduce_for_goal(″AKT1=1, CDK2=1, CDK4=1, CDK6=1,

CycD1=1, CycE1=1, EGF=1, ERalpha=1, ERBB1=1, ERBB1_2=1, ERBB1_3=1, ERBB2=1,

ERBB2_3=1, ERBB3=1, IGF1R=0, MEK1=1, MYC=1, p21=0, p27=0, pRB=1″)

By applying the Cutsets to knock out the genes suggested by PINT, the ERBB network was reduced to 8 entities as shown in Figure 14. It can be observed that all the entities having feedback loops have been retained whereas the ones with straight forward regulations were removed. The AN representation of the reduced network is shown in Figure 15. The reduced network is then modeled as Hybrid Automata in HyTech and time delays are computed for the transitions between various states. The HyTech code and the time delays corresponding to various states of reduced ERBB network are given in the Supplementary Files 3, 4, respectively.

FIGURE 14

Figure 14. Reduced BRN for ERBB receptor-regulated G1/S transition obtained by applying the Cutsets in PINT software (Paulevé, 2017).

FIGURE 15

Figure 15. AN representation of the reduced BRN for ERBB network.

It is noted that ERα is required for the activation of Cyclin D1. At the pivotal point of G1/S transition, cells are dedicated to enter the S phase and results in DNA replication. This process is primarily regulated by CycD1/CDK4/CDK6 complex that phosphorylates and thus, activates the retinoblastoma tumor suppressor protein pRB. For the state 〈1, 0, 1, 0, 1, 1, 0, 1〉 in which the protein pRB is expressed from the initial state 〈1, 0, 0, 0, 0, 0, 0, 0〉, the constraints on the time delays are obtained from the HyTech results as shown in Table 4.

TABLE 4

Table 4. Delay Constraints obtained for the ERBB Receptor-regulated G1/S transition BRN through modeling of its corresponding system of Automata in HyTech.

7. Conclusion and Future Work

In this paper we have proposed a new approach, Hybrid Process Hitting, through which the time delay information has been incorporated in the Process Hitting framework for analysis of large Biological Regulatory Networks. It has been accomplished by integrating the principles of hybrid modeling with those of Process Hitting. The proposed approach utilizes the advantage of Process Hitting in terms of reduction of large state graph which represents the evolution of BRN. This is achieved by considering only the permissible or possible dynamics in the state graph evolution and also by considering the combined influence in case of interactions involving multiple genes.

The feature of Cutsets is utilized for the large BRNs to knock-out some genes to achieve goal-oriented reduction through PINT software tool (Paulevé, 2017). Once the reduced BRN is obtained then the time delay variables for each gene are introduced in it which keeps track of the time taken during each step of BRN evolution. Hence, it becomes possible to obtain the constraints in terms of time delays for each gene corresponding to every evolution step of State Graph.

The time delays found through the application of Hybrid Process Hitting on the evolution of large biological networks gives us very useful information especially at the bifurcation points in the State Graph. While we can see the different paths in the State Graph through static analysis also, however, it cannot be determined which path would evolve before the others. Here, the time delay information becomes highly useful. Once we have this information, we can clearly identify the therapeutic targets to control the progress of State Graph towards desirable goals, i.e., homeostasis conditions.

It is highlighted that the proposed approach takes into consideration the concurrent evolution of all components of BRN and its modeling is realized through Hybrid Automaton. It implies that all the genes are free to evolve simultaneously as allowed by the switching conditions specified by Process Hitting. In this way, true dynamical behavior of BRN evolution from some given initial state to all the possible final (stable) states is obtained.

We have given the complete methodology for this approach with a running example on a toy BRN and it has also been applied to the simpler BRN of Bacteriophage λ and large BRN of ERBB receptor-regulated G1/S transition to show its tractability to large biological networks for computing the time delays associated with the transitional behavior of Network evolution.

Normal regulation of Insulin secretion related pathways pose an important aspect of pancreatic beta cells functionality. Several different types of substrates including metabolites, cytokines and neurotransmitters are involved in regulating insulin secretion through different pathways which also converge with each other at different points in this large interactome.

We plan to apply the proposed Hybrid Process Hitting methodology to Insulin secretion pathways considering all the known stimuli at the same time which would give us a more holistic picture of their regulation in maintaining pancreatic beta cells functional integrity. Overall, our aim is to model the collective dynamics of this insulin secretion interactome in pancreatic beta cells. This could increase our current comprehension in treating Diabetes Mellitus for which the first line of treatment is the use of oral anti-diabetic drugs which mainly target these pathways.

Author Contributions

IS and JA conceived and developed the methodology, conducted the experiments and writing of manuscript. All the authors took part in discussions, analysis and layout of results, and reviewed the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Dr Rehan Zafar Paracha, Dr. Muhammad Tariq Saeed, and Dr. Amnah Siddiqa for their valuable suggestions without which this paper would not have been possible.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.00090/full#supplementary-material

References

Ahmad, J., Bernot, G., Comet, J.-P., Lime, D., and Roux, O. (2006). Hybrid modelling and dynamical analysis of gene regulatory networks with delays. ComPlexUs 3, 231–251. doi: 10.1159/000110010

CrossRef Full Text | Google Scholar

Ahmad, J., Bourdon, J., Eveillard, D., Fromentin, J., Roux, O., and Sinoquet, C. (2009). Temporal constraints of a gene regulatory network: refining a qualitative simulation. Biosystems 98, 149–159. doi: 10.1016/j.biosystems.2009.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahmad, J., Niazi, U., Mansoor, S., Siddique, U., and Bibby, J. (2012). Formal modeling and analysis of the mal-associated biological regulatory network: insight into cerebral malaria. PLoS ONE 7:e33532. doi: 10.1371/journal.pone.0033532

PubMed Abstract | CrossRef Full Text | Google Scholar

Ahmad, J., Roux, O., Bernot, G., Comet, J.-P., and Richard, A. (2008). Analysing formal models of genetic regulatory networks with delays. Int. J. Bioinform. Res. Appl. 4, 240–262. doi: 10.1504/IJBRA.2008.019573

PubMed Abstract | CrossRef Full Text | Google Scholar

Alur, R. (1999). “Timed automata,” in Computer Aided Verification, eds N. Halbwachs and D. Peled (Berlin; Heidelberg: Springer), 8–22.

Google Scholar

Aslam, B., Ahmad, J., Ali, A., Paracha, R. Z., Tareen, S. H. K., Niazi, U., et al. (2014). On the modelling and analysis of the regulatory network of dengue virus pathogenesis and clearance. Comput. Biol. Chem. 53, 277–291. doi: 10.1016/j.compbiolchem.2014.10.003

CrossRef Full Text | Google Scholar

Ben Abdallah, E., Ribeiro, T., Magnin, M., Roux, O., and Inoue, K. (2017). Modeling delayed dynamics in biological regulatory networks from time series data. Algorithms 10:8. doi: 10.3390/a10010008

CrossRef Full Text | Google Scholar

Bernot, G., Cassez, F., Comet, J.-P., Delaplace, F., Müller, C., and Roux, O. (2007). Semantics of biological regulatory networks. Electron. Notes Theor. Comput. Sci. 180, 3–14. doi: 10.1016/j.entcs.2004.01.038

CrossRef Full Text | Google Scholar

Bernot, G., Comet, J.-P., and Khalis, Z. (2008). “Gene regulatory networks with multiplexes,” in Proc of European Simulation and Modelling Conference, ESM'2008 (Le Havre), 423–432.

Google Scholar

de Jong, H. (2002). Modeling and simulation of genetic regulatory systems: a literature review. J. Comput. Biol. 9, 67–103. doi: 10.1089/10665270252833208

PubMed Abstract | CrossRef Full Text | Google Scholar

Folschette, M., Paulevé, L., Inoue, K., Magnin, M., and Roux, O. (2012). “Concretizing the process hitting into biological regulatory networks,” in Computational Methods in Systems Biology. CMSB 2012. Lecture Notes in Computer Science, Vol. 7605, eds D. Gilbert and M. Heiner (Berlin; Heidelberg: Springer), 166–186.

Google Scholar

Folschette, M., Paulevé, L., Magnin, M., and Roux, O. (2015). Sufficient conditions for reachability in automata networks with priorities. Theor. Comput. Sci. 608, 66–83. doi: 10.1016/j.tcs.2015.08.040

CrossRef Full Text | Google Scholar

Henzinger, T. A., Ho, P.-H., and Wong-Toi, H. (1997). Hytech: a model checker for hybrid systems. Int. J. Softw. Tools Technol. Transfer 1, 110–122. doi: 10.1007/s100090050008

CrossRef Full Text | Google Scholar

Khalis, Z., Comet, J.-P., Richard, A., and Bernot, G. (2009). The smbionet method for discovering models of gene regulatory networks. Genes Genomes Genomics 3, 15–22. Available online at: http://www.i3s.unice.fr/~comet/publications/RevuesIntl/articleGGG2009-final.pdf

Google Scholar

Paulevé, L. (2017). “Pint: a static analyzer for transient dynamics of qualitative networks with IPython interface,” in CMSB 2017-15th Conference on Computational Methods for Systems Biology (Darmstadt: Springer International Publishing), 309–316. doi: 10.1007/978-3-319-67471-1_20

CrossRef Full Text

Pauleve, L. (2018). Reduction of qualitative models of biological networks for transient dynamics analysis. IEEE/ACM Trans. Comput. Biol. Bioinform. 15, 1167–1179. doi: 10.1109/TCBB.2017.2749225

PubMed Abstract | CrossRef Full Text | Google Scholar

Paulevé, L., Andrieux, G., and Koeppl, H. (2013). “Under-approximating cut sets for reachability in large scale automata networks,” in Computer Aided Verification, Lecture Notes in Computer Science, vol. 8044 eds N. Sharygina and H. Veith. (Berlin; Heidelberg: Springer), 69–84.

Paulevé, L., Magnin, M., and Roux, O. (2011). “Refining dynamics of gene regulatory networks in a stochastic π-calculus framework,” in Transactions on Computational Systems Biology XIII. Lecture Notes in Computer Science, Vol 6575, eds C. Priami, R. J. Back, I. Petre, and E. de Vink (Berlin; Heidelberg: Springer), 171–191.

Richard, A., Comet, J. P., and Bernot, G. (2006). “Formal methods for modeling biological regulatory networks,” in Modern Formal Methods and Applications, ed H. A. Gabbar (Dordrecht: Springer), 83–122.

Google Scholar

Saeed, M. T., Ahmad, J., Kanwal, S., Holowatyj, A. N., Sheikh, I. A., Zafar Paracha, R., et al. (2016). Formal modeling and analysis of the hexosamine biosynthetic pathway: role of o-linked n-acetylglucosamine transferase in oncogenesis and cancer progression. PeerJ 4:e2348. doi: 10.7717/peerj.2348

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahin, Ö., Fröhlich, H., Löbke, C., Korf, U., Burmester, S., Majety, M., et al. (2009). Modeling ERBB receptor-regulated g1/s transition to find novel targets for de novo trastuzumab resistance. BMC Syst. Biol. 3:1. doi: 10.1186/1752-0509-3-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, Y. (2012). A lie algebra approach to susceptible-infected-susceptible epidemics. Electron. J. Differ. Equat. 2012, 1–7. Available online at: https://ejde.math.txstate.edu/Volumes/2012/233/abstr.html

Google Scholar

Shang, Y. (2015a). Analytical solution for an in-host viral infection model with time-inhomogeneous rates. Acta Phys. Pol. B 46, 1567–1577. doi: 10.5506/APhysPolB.46.1567

CrossRef Full Text | Google Scholar

Shang, Y. (2015b). Deffuant model of opinion formation in one-dimensional multiplex networks. J. Phys. A Math. Theor. 48:395101. doi: 10.1088/1751-8113/48/39/395101

CrossRef Full Text | Google Scholar

Shang, Y. (2017). On the delayed scaled consensus problems. Appl. Sci. 7:713. doi: 10.3390/app7070713

CrossRef Full Text | Google Scholar

Shang, Y. (2018). Hybrid consensus for averager–copier–voter networks with non-rational agents. Chaos Solit. Fract. 110, 244–251. doi: 10.1016/j.chaos.2018.03.037

CrossRef Full Text | Google Scholar

Sheikh, I. A., Ahmad, J., and Saeed, M. T. (2016). “Modelling and simulation of biological regulatory networks by stochastic petri nets,” in Proceedings of the World Congress on Engineering and Computer Science (San Francisco, CA).

Google Scholar

Siebert, H., and Bockmayr, A. (2008). Temporal constraints in the logical analysis of regulatory networks. Theor. Comput. Sci. 391, 258–275. doi: 10.1016/j.tcs.2007.11.010

CrossRef Full Text | Google Scholar

Thieffry, D., and Thomas, R. (1995). Dynamical behaviour of biological regulatory networks–II. immunity control in bacteriophage lambda. Bull. Math. Biol. 57, 277–297.

PubMed Abstract | Google Scholar

Thomas, R. (1978). Logical analysis of systems comprising feedback loops. J. Theoret. Biol. 73, 631–656. doi: 10.1016/0022-5193(78)90127-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, R. (2013). “Kinetic logic: a Boolean approach to the analysis of complex regulatory systems,” in Proceedings of the EMBO Course “Formal Analysis of Genetic Regulation”, (Brussels: Springer Science & Business Media).

Google Scholar

Thomas, R., and d'Ari, R. (1990). Biological Feedback. Boca Raton, FL: CRC Press.

Google Scholar