Research | Open | Published: 08 January 2016
On the modeling and analysis of the biological regulatory network of NF-\({\kappa }\)B activation in HIV-1 infection
Zurah Bibi1,
Jamil Ahmad1,
Amjad Ali2,
Amnah Siddiqa1,
Shaheen Shahzad3,
Samar HK Tareen1,
Hussnain Ahmed Janjua2 &
Shah Khusro4
Complex Adaptive Systems Modelingvolume 4, Article number: 1 (2016) | Download Citation
The complex interactions between genetic machinery of HIV-1 and host immune cells mediate dynamic adaptive responses leading to Autoimmune Deficiency Syndrome. These interactions are captured as Biological Regulatory Network (BRN) which acts to maintain the viability of host cell machinery through feedback control mechanism which is a characteristic of complex adaptive systems. In this study, the BRN of immune response against HIV-1 infection is modeled to investigate the role of NF-κB and TNF-α in disease transmission using qualitative (discrete) and hybrid modeling formalisms.
Qualitative and Hybrid modeling approaches are used to model the BRN for the dynamic analysis. The qualitative model is based on the logical parameters while the hybrid model is based on the time delay parameters.
The qualitative model gives useful insights about the physiological condition observed as the homeostasis of all the entities of the BRN as well as pathophysiological behaviors representing high expression levels of NF-κB, TNF-α and HIV. Since the qualitative model is time abstracted, so a hybrid model is developed to analyze the behavior of the BRN by associating activation and inhibition time delays with each entity. HyTech tool synthesizes time delay constraints for the existence of homeostasis.
Hybrid model reveals various viability constraints that characterize the conditional existence of cyclic states (homeostasis). The resultant relations suggest larger cycle period of HIV-1 than the cycle periods of the other two entities (NF-κB and TNF-α) to maintain a homeostatic expressions of these entities.
Autoimmune Deficiency Syndrome (AIDS) is a major pandemic, caused by the human immunodeficiency virus (HIV). HIV targets CD4+ T lymphocytes and hinders the regulation and amplification of immune cells (Douek et al. 2002). This infection results in the persistent activation of immune system, and consequently, CD4+ T lymphocytes would greatly reduce in number resulting in body's ability to fight against HIV infection. There is a need to understand the mechanism of immune activation and how it gets repressed by viral overload so the major elements of immune system that plays an important role are discussed below in detail.
Tumor necrosis factor (TNF), also known as cachectin, is a cytokine primarily produced in activated immune system, and causes inflammation in acute infections (Legler et al. 2003). It is produced by activated macrophages, lymphocytes and monocytes elicits transcription as well as production of the inhibitor of kappa B (I\(\kappa\)B) kinase enzyme after linking to corresponding receptors that results in the production of transcription factor NF-κB (Vitale and Ribeiro 2007). NF-κB, after activation, induces cytokines (mainly TNF-α and IL-1\(\beta\)) and causes inflammation (Bingham 3rd 2002) to hinder the entry of viral proteins into the cellular environment and initiates different immunological processes (Bingham 3rd 2002). Studies indicate that TNF-α specifically attacks HIV-infected cells and "boosts" immune system (Herbein et al. 1996; Lane et al. 1999) clearly depicted in Fig. 1. Few in-vitro studies prove that TNF-α stimulates the production of HIV mRNA indicating that TNF-α may be dose and time dependent on viral replication (Vitale and Ribeiro 2007). Macrophages serve as a reservoir for the virus and help the virus to propagate Swami (2013) as the HIV-1 virus releases its soluble exogenous proteins (NEF, TAT, VPR, GP120) to enter the macrophages, leading to transformation of both cellular and viral transcriptional machinery. Viral glycoprotein 120 (GP120) is essential for HIV-1 viral entry by interacting with surface receptors CD4 (Cluster of differentiation 4) and CXCR4 (CXC chemokine receptor 4) or CCR5 co-receptor (chemokine receptor 5) (Fig. 1) (Herbein et al. 1994, 2010; Melikyan 2008).
HIV signaling and activation of Immune system in Macrophages (M\(\Phi\)). HIV proteins modify the signaling pathways in both infected and uninfected cells by activating the transcription factors like NF-κB and AP-1. NF-κB has binding sites for both HIV LTR (long terminal repeat) as well as for TNF-α promoter (Adapted from Herbein et al. 2010)
Another important protein that plays significant role both in innate and adaptive immune response is the transcription factor NF-κB. In the absence of foreign antigen(s) NF-κB dimers are retained in the cytoplasm through the inhibitory action of the IκB molecules (Ghosh and Karin 2002). The pro-inflammatory cytokines TNF-α and Interleukin-1\(\beta\) (IL-1\(\beta\)) induce the activation of NF-\(\kappa\)B pathway accompanied by rapid degradation of I\(\kappa\)B-\(\alpha\), releasing NF-\(\kappa\)B dimers to the nucleus (Grey 2008), where they activate HIV transcription in infected cells (Ghosh and Karin 2002).
The immune regulators, such as cytokines, interact with the viral proteins in a feedback mechanism (see Fig. 1) to maintain homeostatic behavior of the organism. Homeostasis is one of the vital creeds of self-regulation and internal organization and is defined as the synchronized physiological reactions responsible for maintaining most of the healthy states and conditions in the body. Positive and negative feedback loops induce an on/off switching mechanism (Brandman and Meyer 2008), which stops the respective activation state after the arrival of a set point or threshold (Thomas and Thieffry 1995; Thomas et al. 1998). Such feedback mechanisms represent a characteristic feature of complex adaptive systems which facilitates the emergence of complex patterns of behaviors such as homeostasis and stable states (reference: Complex Adaptive Systems Modeling: A multidisciplinary Roadmap). These positive and negative interactions among biological entities comprise a Biological Regulatory Network (BRN) providing the basis to converge into one of the many possible outcomes. The modeling and analysis of expression dynamics of BRNs can be used to increase our understanding of the underlying adaptive mechanisms in response to various stimuli but the presence of large number of unknown biological parameters e.g., kinetic rates of reactions and the expression values of molecules etc. makes it challenging. The kinetic logic modeling formalism of Rene Thomas is one of the established methods to perform qualitative analysis of such systems (Thomas and Thieffry 1995; Gagneur and Casari 2005; De Jong et al. 2004) without any prior knowledge of kinetic parameters. The parameters of the model are computed using model checking technique based on the biological observations (see "Methods"). These parameters are then used to construct a qualitative model (state graph) which captures all the possible behaviors of the BRN as cyclic trajectories or paths diverging to the stable states. Rene Thomas has proved the similarity of differential equation models and qualitative models based on kinetic logic formalism in ([55] of rehan sepsis). However, qualitative modeling is more appropriate modeling paradigm due to its ability to infer unknown parameters ([56] of rehan sepsis). This formalism has been successfully implemented to gain insights into the pathogenesis clearance of Dengue virus and identification of the putative drug targets in Cerebral Malaria and Sepsis.
Several signaling pathways have been carefully reviewed to construct a BRN based on HIV intervention of immune system (Fig. 1). In the current study, the qualitative model developed by Bibi et al. (2011) is employed to investigate the mode of immune activation in hindering HIV-1 viral proteins rapid expression. This approach facilitates in obtaining useful constraints that provides necessary conditions for the BRN to evolve and helps to analyze immune response of the body against foreign invaders. The HyTech model checker tool is used to analyze the hybrid model of the BRN (Ahmad et al. 2007) that generates delay constraints (Fromentin et al. 2010) for the existence of a particular cyclic behavior (representing homeostasis). The model findings are in line with experimental data and provides novel insights into the disease propagation mechanism which could be used for its intervention strategies.
Discrete modeling formalism
The discrete modeling formalism of René Thomas (Ahmad et al. 2012, 2007, 2009; Ahmad and Roux 2008; Thieffry and Thomas 1998; Thomas and D'Ari 1990), a type of qualitative logical modeling, has been successfully used for the construction of BRN models. This formalism uses a set of parameters that reflect the behaviors of a BRN in terms of the presence or absence of the regulatory elements (resources). Following are the formal definitions which are used to construct a discrete model of a BRN (Ahmad et al. 2012).
In a directed graph \(G = (X, A)\), where X is a set of vertices (nodes) and A is a set of edges linking the vertices, \(G^- (v)\) and \(G^+(v)\) represent the set of predeccessors and successors of a vertex \(v \in X\) respectively.
(Biological Regulatory Network) A BRN is a graph \(G = (X, A)\) where X represents the set of nodes representing biological entities and A is the set of edges representing interactions between biological entities. Each edge \(a \rightarrow b\) is labeled as (s\(_{ab}\), \(r_{ab}\)), where s\(_{ab}\) is a positive integer representing a threshold and \(r_{ab}\) \(\in \{+, -\}\) shows the type of interactions ('\(+\)' for activation and '\(-\)' for inhibition).
A graphical representation for the activation and inhibition of v \(\in\) X is given in Fig. 2, where the rate of activation or inhibition depends upon a certain threshold level (\(\theta\)). Production of v increases rapidly once its activator \(u\in X\) reaches its threshold (\(\theta\)). Likewise, the degradation of v occurs at maximum rate once its inhibitor \(w\in X\) reaches its threshold level (\(\theta\)). The threshold \(\theta\) discretizes the expression of an activator (u) or inhibitor (w) into discrete (qualitative) levels represented by positive integers.
Activation and inhibition of an entity v are represented by sigmoid curves. Step function represents its discrete approximation. The values 0 and 1 represent the discrete (qualitative) expression levels
To analyze the behavior of a BRN, it is necessary to know all the possible states and the transitions between them. There is a limit discrete expression value \(lm_{v}\) for each entity v of a BRN which is equal to its outgoing degree, such that \(\forall\) \(u \in G^+(v)\) each \(s_{uv}\) \(\in \{1,\ldots ,n_{v}\) \(\}\) where \(n_{v}\) \(\le\) \(lm_{v}\). Each entity carries its expression levels in the set \(Q_{v}= \{0,\ldots ,n_{v}\}.\)
(State) A state is a configuration of the expression levels of all entities at any instance of time and it is denoted by a tuple \(s \in \S\), where:
$$\begin{aligned} \S = \Pi _{ v \in X} Q_{v}. \end{aligned}$$
A vector is normally used to show a qualitative state \((l_{v})_{v \in X}\), where \(l_{v}\) represents the expression level of the entity v.
A set of resources represents the activators of an entity at any instant of time.
(Resources) The set of resources R for a variable \(v\in X\) at some level \(l_v\) is given as follows.
$$\begin{aligned} R_{l_v} = \{u \in G^-(v) | (l_{u} \ge s_{uv} ~and~ r_{uv} = `+\text{'}) ~or~ (l_u < s_{uv} and~ r_{uv} = `-\text{'} )\}. \end{aligned}$$
From the above definition, it can be inferred that the absence of an inhibitor below threshold is considered as a resource.
The set of logical parameters assigned to each biological entity determine the dynamics of a BRN.
(Set of Logical Parameters) The set of logical parameters K(G) of a BRN is given below.
$$\begin{aligned} \begin{array}{l} K (G) = \{K_{v,R_{l_v}} \in \{0,...n_{v}\} ~|~ l_{v} \in Q_{v} ~\forall ~v \in X\}. \end{array} \end{aligned}$$
\(K_{v, R_{l_v}}\) define the target expression level towards which the expression of entity v evolves.
Let l and \(K \in Z\,{\ge } \, 0\), the asynchronous evolution operator \(\Rsh\) (Bernot et al. 2004) is given as:
$$\begin{aligned} l\Rsh K ={\left\{ \begin{array}{ll} {l} + 1 ~ if~ {l} < K\\ {l} - 1 ~if ~ {l} > K \\ {l} ~if ~l = K \\ \end{array}\right. } \end{aligned}$$
Depending on the values of the logical parameters, the dynamics of a BRN is defined as a state graph.
(State Graph) If \(l_v\) is the level of an entity v in a state \(s \in \S\), the state graph of a BRN with transition relation \(T \subseteq \S \times \S\) is a directed graph such that \(s \rightarrow s' \in T\) iff:
there is a unique \(v \in X\) such that \(s_{{l}_{v}} \ne s'_{l_v}\) and \(s'_{l_v} = s_{l_v} \Rsh K_{v, R_{l_v}}\);
and \(s_{l_u} = s'_{l_u} \forall u \in X \setminus {v}\).
A state graph differs from its successor state by one component only, so if a state s has n elements to be evolved then it will have at most n successor states.
The discrete modeling formalism can be easily understood by a toy example. Consider the BRN shown in Fig. 3 consisting of two genes x and y (represented by x and y nodes), where y acts as an activator of x while y favors its own production (self-stimulation) and is inhibited by x. Even for arbitrary BRNs, such as this one, different sets of logical parameters can be inferred, leading to difficulties in selection of the parameters. A methodology pioneered by Bernot et al. (2004), utilises the application of model checking over the state graphs generated by the possible sets of logical parameters. Model checking is an exhaustive exploratory approach applicable to graphs to check them for certain properties expressed under Computation Tree Logic (CTL) (Clarke et al. 1999). The details of the syntax and applicability of CTL, along with the explanation of the tool SMBioNet are given in detail by Khalis et al. (2009). The SMBioNet tool was applied on this BRN, using the CTL formula: \((y=2)\rightarrow AXAF(y=2)\), (where A enforces execution of the given CTL in all paths and F imposes execution of the given CTL in at least one path Aslam et al. 2014) which states that once the protein y has achieved the level 2, it will remain at that level and induce mucus production. The application of this CTL generated a parameter set, given in state table (Table 1), which also shows the transitioning of states from a given state. Its state graph is also shown in Fig. 4, representing the behavior of the BRN by one cycle \((0, 0) \rightarrow (0, 1) \rightarrow (1, 1) \rightarrow (1, 0)\) and stable steady state (1, 2).
Table 1 The state table of genes x and y
A toy BRN consisting of genes x and y where y activates x (shown by positive sign) and x inhibits y (shown by negative sign). Self-loop of y shows that y is self-activating
The state graph of x and y genes; representing cyclic (homeostatic) and stable state (deadlock) behaviors of the system. The system evolves from silenced state (x,y) (0,0) leading towards (0,1) as in the absence of the inhibitor x, y is activated. Once the threshold value of 2 of y in the state (0, 2) is achieved (constant activation of y), it takes the system towards deadlock state (1, 2)
Hybrid modeling
In hybrid modeling, both discrete and continuous features are combined to study diverse biological properties in detail (Ahmad et al. 2009). In this modeling formalism, a sigmoid-curve is discretized into piece-wise linear curve (Ahmad et al. 2007) as illustrated in Fig. 5. The delays governing the gene evolution are characterized as time intervals, which are measured using real-valued clocks (Ahmad et al. 2009).
Evolution of clock variable \(h_v\). Here, \(dh_v\)/dt = 1 represents the rate of evolution of clock \(h_v\) associated with gene v, which evolves from low level to higher level 1 Ahmad et al. (2008). When a delay is measured then the respective clock is reset (\(h_v \leftarrow 0\))
In a BRN each gene is associated a clock to record evolution time of its expression. The clock is reset to zero after each transition whenever the expression passes from one level to another (Ahmad et al. 2009). In the current BRN, delays are taken as unvalued parameters to introduce the concept of Parametric Bio Linear hybrid Automaton (Bio-LHA) (Ahmad et al. 2007, 2009, 2012; Ahmad and Roux 2008) that could be analyzed for certain delay constraints describing the stay conditions for the behaviours of BRN.
We denote by V and \(\mathbbm {P}\), the sets of real valued variables and parameters, respectively. Let \(C =(V, \mathbbm {P})\), \(C \le (V, \mathbbm {P}))\) and \(C \ge (V, \mathbbm {P}))\) represent the set of atomic constraints using only =, ≤ and ≥, respectively.
(Parametric Bio-Linear Hybrid Automaton) A Bio-Linear Hybrid Automaton is a 6-tuple \(H=\left\langle L, l_{0}, \mathbbm {P}, V, E, Inv, Diff \right\rangle\) where:
L is a finite set of locations;
\(l_{0}\) is the initial location;
\(\mathbbm {P}\) is the finite set of parameters (delays);
V is a finite set of real-valued variables (clocks);
\(E\subseteq L \times C = (V, \mathbbm {P}) \times 2^{v} \times L\) is a finite set of edges. \(e = \{l ,g, \mathbbm {R}, l'\}\) \(\in E\) is the edge between the locations l and \(l'\) with the guard g and the reset set \(\mathbbm {R} \subseteq V\) contains the set of clocks which appears in g.
Inv assigns invariant condition as a conjunction of atomic constraints \(C \le (V, \mathbbm {P})) ~ and ~C \ge (V, \mathbbm {P}))\) to each location \(l \in L\);
Diff \(L \times V \rightarrow \{-1, 0, 1\}\) assigns an evolution rate to each continuous variable \(v \in V\) (clocks) in each location \(l \in L\).
The semantics of a parametric Bio-LHA is a timed transition system. It is defined according to the time domain T. We let T* \(= T\) \\(\{0\}\).
(Semantics of Bio-Linear Hybrid Automaton) Let \(\nu\) be a valuation for the delay parameters \(\mathbbm {P}\). The \((T,\nu )\)-semantics of a parametric Bio-LHA
$$\begin{aligned} H=\left\langle L, l_{o}, v, \mathbbm {P}, e, Inv, Diff \right\rangle \end{aligned}$$
is defined as a timed transition system \(S_\mathbbm {H} = \{S,s_i,T \rightarrow \}\) where: (1) \(S={\{(l, \nu )\} | l \in L ~and~ \nu \models Inv (L)}\); (2) \(s_i=(l_0, v_0)\) represents initial state (a configuration of initial location and initial values of real valued variables clocks); (3) the relation \(\rightarrow \subseteq S \times T \times S\) is defined for \(t \in T\) as:
Discrete transitions: \((l, v) \rightarrow (l', v')\) iff \(\exists ~\{l, g, \mathbbm {R}, l'\} \in E\) such that \(g(v)=\) true, \(v'(h)=0\) if \(h \in \mathbbm {R} ~and~ v'(h)=v(h) ~if~ h \notin \mathbbm {R}\).
Continuous transitions: For \(t \in\) T*, \((l, v) \rightarrow (l', v')\) iff \(l'= l, v'(h)= v(h) + Diff(l, h) \times t\), and for every \(t' \in [0,t], (v(h) + Diff(l, h) \times t' \models Inv(l)\).
(Temporal Zone) A zone in which the discrete level of entities remains the same, however the clocks of the entities keep ticking until the evolution (activation or inhibition) of an entity takes place (Ahmad et al. 2008, 2009).
(Temporal State Space) The temporal state of a BRN comprises of the complete set of temporal zones derived from the discrete model of the respective BRN.
In the hybrid model of a BRN, we denote \(\phi (t)\) for \(t\in \mathbb {R}_{\ge 0}\), while the sequence of points of a trajectory and the set of all points in its state pace is denoted by S (Ahmad et al. 2009). A particular trajectory is said to be viable, iff it remains within a prescribed region known as its viability domain (Ahmad and Roux 2008).
Definition 10
(Invariance kernel) An invariance kernel \(K\in S\) is the largest invariant subset of S iff for any point \(p\in K\), any trajectory starting in p is viable in K.
For the analysis of invariance kernel, the notion of Period (denoted by \(\pi )\) (Ahmad et al. 2009) is introduced.
(Period) It is defined as sum of all the delays that a gene goes through sequentially and successively, once through each of all its expression levels in a cycle (Ahmad et al. 2009). It is denoted as \(\pi\).
HyTech (Henzinger et al. 1997) is used to analyze the invariance kernel and convergence domain, the results of which are expressed as linear constraints that specify the conditions under which a particular behaviour exists. (Ahmad et al. 2007, 2009; Ahmad and Roux 2008).
Discrete and hybrid modeling of the NF-\(\kappa\)B associated BRN
The regulatory pathway in Fig. 1 is represented by a labeled directed graph as shown in Fig. 6. Qualitative model of the BRN is derived by applying discrete formalism. For this purpose only qualitative data of the network is required, i.e., sign of interactions (+ or −), threshold order and weights for each interaction. For situations where the thresholds are weakly estimated, the logical approach could be an alternative (Thieffry and Thomas 1998).
BRN of NF-\(\kappa\)B responsible for the cytokine production and regulation. HIV stimulates NF- \(\kappa\)B and TNF-\(\alpha\) production while TNF-\(\alpha\) in turn regulates HIV and NF-\(\kappa\)B that helps HIV to replicate itself. Threshold of 2 represents that HIV wins the battle to trans-activate itself even at low expression level of NF-\(\kappa\)B immediately after entering into cellular environment while TNF-\(\alpha\) will build at later points depending upon increased or greater levels of NF-\(\kappa\)B
In the feedback circuit (see Fig. 6) only the main regulators (HIV-1, TNF-\(\alpha\) and NF-\(\kappa\)B) are included and the logical formalism has been applied to observe the behaviors like multistationarity and homeostasis (Gouze 2009). The set of parameters is defined in Table 2, derived after performing several experiments (refer to Bibi et al. 2011 for detailed BRN understanding).
Table 2 State transition table: transition from one state to the next state is given, here activators act as resources while absence of an inhibitor is also taken as its resource that decides its K Parameter (e.g. absence of TNF-\(\alpha\) in Ist row act as a resource for HIV that set K Parameter for HIV as 1)
GENOTECH tool (Ahmad 2009) is used to observe the dynamical behavior of the BRN through the generation of steady states. This tool facilitates the discrete modeling of a BRN. Apart from the discrete modeling, it also aids in analysis of steady state behaviors (cyclic and dead lock states). A hybrid model is also obtained by GENOTECH that integrates clocks to compute production and degradation delays for each entity in the BRN. HyTech tool is then applied to analyze this hybrid model which generates parametric constraints or stay conditions to explain the existence of a particular state (or behavior) (Fromentin et al. 2010). The work flow is sequentially illustrated in (Fig. 7).
Work flow. Step-wise representation of methodology applied to construct BRN and subsequent state graph analysis using Graphviz tool. HyTECH model is finally developed to incorporate time-delays for the activation and inhibition of each entity
Parameter synthesis and Qualitative Model of NF-\(\kappa\)B associated BRN
The NF-\(\kappa\)B associated BRN (illustrated in Fig. 6) was first modelled in SMBioNet as an input to generate the preferred logical parameter set. Trivial parameters of each entity, the resource sets of which contained only activators or inhibitors, were restricted to '\(n_v\)' and '0' respectively. The CTL formulas used were:
$$\begin{aligned} \begin{array}{cccc} (HIV=0\,\wedge \,TNFA=0\,\wedge \,NFkB=0) &{} \rightarrow &{}\\ EF(AG(HIV=1,\wedge \,TNFA=1\,\wedge \,NFkB=2)) &{} \wedge \\ (HIV=0\,\wedge \,TNFA=0\,\wedge \,NFkB=0) &{} \rightarrow &{}\\ EX(EF(HIV=0\,\wedge \,TNFA=0\,\wedge \,NFkB=0)) &{} \wedge \\ (HIV=0\,\wedge \,TNFA=1\,\wedge \,NFkB=2) &{} \rightarrow &{}\\ AF(AG(HIV=1\,\wedge \,TNFA=1\,\wedge \,NFkB=2)) &{} \\ \end{array} \end{aligned}$$
The first formula states that from an arbitrary initial state where all three entities have the level '0', there should exist a trajectory where the system enters a deadlock state in which the HIV proteins are at level '1', TNF-\(\alpha\) is at level '1', and NF-\(\kappa\)B is at level '2'. Likewise, the second formula starts from the same arbitrary initial state and enforces that there should exist at least one trajectory from which this initial state is again reached, meaning that there should exist at least one cycle from this state. The third and final formula states that once TNF-\(\alpha\) and NF-\(\kappa\)B are at their respective \(n_v\) in the absence of HIV, then the system will always end up in the diseased steady state. These formulas generated 3 parameter sets, which differed only at the parameter \(K_{NF\mathrm {-}\kappa B}\{HIV\}\). As the trend in all the parameter sets was that of NF-\(\kappa\)B being induced in the presence of HIV, we selected the parameter set which allowed NF-\(\kappa\)B to achieve its maximum threshold level of '2' on this logical parameter assuming that full activation (optimal expression) of NF-\(\kappa\)B is lethal for infected macrophages to fight against HIV. The selected parameter set enables all the entities to interact with each other corresponding to naturally occurring phenomenon while maintaining interdependencies of participating entities upon each other (for stimulation or suppression). In addition, the logical parameters are given in Table 2.
State graph
Based on the selected logical parameters, a state graph is generated as shown in Fig. 8. There are two types of states in the graph, the states with certain incoming and outgoing arrows known as steady states while there is a state (1,1,2) with only incoming arrows (no outgoing arrow) it represents stable steady state (in red). A close analysis of the state graph reveals two cycles exiting within it (in green) that represent a probable mechanism of homoeostasis once the system encounters an HIV infection representing that how the BRN behaves after viral entry and recovery control mechanism to fight against HIV infection. In fact, these cyclic states show the periods of clinical latency, the period when virus stimulates the immune system to produce antibodies against HIV-1. The stable steady state represents a diseased state (the onset of AIDs) where all the proteins are over expressed. A total of 6 paths lead towards stable steady state (1,1,2) from the initial state (1,0,0), which indicates HIV infection. All these paths enter the homeostatic cycle via the state (1,1,0) and divert towards stable steady state in case of violation of the condition required to be in the cycle (discussed in the "Hybrid modeling results" section). The complete model of the system given in Fig. 8 shows cyclic states along with the steady state where the cyclic states are associated with green arrows while red arrows lead towards stable steady state.
State graph of the BRN shows a complete set of states. Each state represents the discrete expression levels of the entities in the order of HIV, TNF-\(\alpha\), NF-\(\kappa\)B. Green arrowsrepresent cyclic trajectories maintaining the recovery mechanism. Whereas, the paths leading into stable steady state (1,1,2) are represented by red dashed arrows
Real-time hybrid model
The qualitative modeling process generated all possible discrete states (Fig. 8) and their transitions which give further insights into the steady states behaviors. These steady states (cycles and stable states) are in agreement with the general hypothesis of wet-lab and clinical experiments. After applying the discrete modeling formalism on NF-\(\kappa\)B associated BRN, its hybrid model is obtained for the BRN with help of parametric constraints. Hybrid modeling allowed us to compute the necessary and sufficient conditions for the existence of different behaviors.
The BRN under study generated two cycles in HyTech, having some overlapping cyclic states. The presence of inequality or equality in the resulting constraints determines the stability or instability of the cycle respectfully. Of the two cycles, only one was able to converge to its invariance kernel, and is thus the only one analysed (see Table 3 for the cycle and delay constraints).
Table 3 Constraints for the invariant cycle
By varying time delays, i.e., the delay required by an entity to evolve or degrade, system dynamics can be accurately observed. This new model is then analyzed with linear hybrid model checker tool HyTech (Henzinger et al. 1997) in order to extract the constraints for cyclic states, the significant ones of which are discussed below:
Constraint 1: \(\delta ^+_{NFKB010}+|\delta ^-_{NFKB001}|+2|\delta ^-_{TNF-A011}|=\delta ^+_{HIV000} +|\delta ^-_{HIV110}|\) This equality constraint revealed the following relation between HIV and NF-\(\kappa\)B.
$$\begin{aligned} \pi (HIV) \ge \pi (NF{\text -}\kappa B) \end{aligned}$$
where \(\delta ^+_{NFKB010}\) represents production delay of NF-\(\kappa\)B in the logical state 010 while \(\delta ^-_{NFKB001}\) represents degradation delay of NF-\(\kappa\)B in the logical state 001 and \(\pi\) represents time period (sum of positive and negative time delays).
These relations capture the general trend of HIV infection and immune stimulation where the period of HIV infection is greater than the period of NF-\(\kappa\)B, this is graphically represented in Fig. 9b. This constraint shows slow propagation of viral infection that strengthens the immune system and hence lowers the risks of viral attack for surrounding cells (Bosque 2009). It is important to note that the real trend of a gene expression may require more time (as it may go through lazy stages where no increase or decrease occurs) to complete a period (Ahmad et al. 2009).
Dummy/graphical representation of HIV time period relative to TNF-\(\alpha\) a and NF-\(\kappa\)B b
Constraint 2: \(\delta ^+_{TNF-A100}\)+\(|\delta ^-_{TNF-A011}|\) \(\le \delta ^+_{HIV000}\)+\(|\delta ^-_{HIV110}|\) The relation between HIV and TNF-\(\alpha\) is also extracted as follows.
$$\begin{aligned} \pi (HIV) \ge \pi (TNF{\text{-}}\alpha ) \end{aligned}$$
This constraint is explained by the fact that cytokine production is boosted exponentially early in immune attack to combat viral attack by HIV-1. Here the delay associated with TNF-\(\alpha\) for its expression and degradation (\(\delta ^+\) production delay at state 100 and \(\delta ^-\) degradation delay at state 011) is less than the complete period of HIV at states 000 and 011. Thus, the qualitative period of TNF-\(\alpha\) is traversed faster the period of HIV, in an attampt to hamper HIV proteins entry into surrounding macrophages. It is evidenced by severe inflammation to control initial burst of viral proliferation (represented in Fig. 9a). It showed that after the release of HIV proteins into the cellular machinery, production of cytokines occurs rapidly as compared to viral replication, i.e., the period of HIV ends after the period of TNF-\(\alpha\). It leads to the condition where TNF-\(\alpha\) reaches its destination earlier as compared to HIV in order to activate or enhance the immune system.
Constraint 3: \(\delta ^+_{HIV000}\) + \(|\delta ^-_{HIV110}|\) \(\le \delta ^+_{NF-KB110}\) + \(|\delta ^-_{NF-KB001}|\)+ \(2|\delta ^-_{TNF-A011}|\) After getting infected with HIV, some patients remain healthy for some years thats shows clinical latency period. However, the virus somehow manages to progress despite the active immune system. This condition is a clear depiction of such a scenario where HIV's qualitative period progresses ahead of NF-\(\kappa\)B and TNF-\(\alpha\) causing failure of immune system. This constraint is an important candidate to represent critical points in later stages of HIV infection having inefficient T cells response leading to immune escape of HIV. It illustrates the situation when HIV gets into the immune system and starts proliferating its number via the transcription factor NF-\(\kappa\)B, immune system takes more time delay to produce inflammatory cytokines (TNF-\(\alpha\)) as compared to the viral propagation (production) (Fig. 10). This time delay is generally followed after the period of clinical latency has passed.
Fig. 10
Dummy/graphical representation of period of TNF-\(\alpha\) relative to HIV
In this study, qualitative and hybrid models of the BRN for immune activation (upon HIV infection) were constructed based on the multi-valued logical formalism of Thomas, to estimate the important dynamic properties characterized by regulation delays. On the basis of HyTech results, two important relations are drawn from the set of constraints: \(\pi (HIV)\ge \pi (NF\text{-}\kappa B)\) and \(\pi (HIV) \ge \pi (TNF\text{-}\alpha )\).
In order to maintain the homeostatic balance of human immune defense, the invasion period of HIV must span longer duration as compared to TNF-\(\alpha\) and NF-\(\kappa\)B during the initial phases of viral attack. This will tend to reduce the viral proteins ability to start its replication rapidly, giving strength to body's innate immunity to fight against HIV-1 pathogenesis. It reflects an important relationship between HIV, TNF-\(\alpha\) and NF-\(\kappa\)B. When the progression of HIV proteins is slow in comparison to immune defense mechanism, the virus normally takes much longer to infect and override the whole immune system of the infected individual. This phenomenon is followed by rapid generation of cytokines by cellular transcriptional machinery NF-\(\kappa\)B while there is a considerably slow replication and invasion of HIV in the body. It shows frequent evolution of NF-\(\kappa\)B and TNF-\(\alpha\) and relatively stable behavior of HIV in comparison to the other two. However, if the fore mentioned constraint is violated, the system would get stuck in the stable steady state (1,1,2) shown in Fig. 8, where the chances of survival diminishes. This situation is generally seen during the period of acute infection. Keeping this in view, longer survival time of patients (clinical latency periods) are observed in case of longer qualitative period of HIV than that of immune defense elements (i.e., NF-\(\kappa\)B and TNF-\(\alpha\)). Furthermore, at initial stages of HIV infection, high-levels of CD8+ T cell frequencies have been observed to correlate with the control of viral replication (Wilson et al. 2000; Allen et al. 2000; Davenport et al. 2004). From this, it can be concluded that while designing drugs against HIV-1, the delay time for HIV-1 proteins' expression must be kept greater than those of NF-\(\kappa\)B and TNF-\(\alpha\) by administrating greater dose of cytokines (TNF-\(\alpha\)).
When the time period for HIV's protein expression is lesser, it indicates that the system is attainting the stage of viral setpoint. At this stage, number of TNF-\(\alpha\) decreases in T lymphocytes and this can be seen by the shifting of the initial steady state of TNF-\(\alpha\) (cytokines) to a new equilibrium value. This new equilibrium value is lower than the preinfection value (Fauci 1996). At this stage, an exponential increase in virus load may be observed, reaching a peak and finally declining to the steady state level also called as viral setpoint (Kaufmann et al. 1998; Lindbäck et al. 2000). HIV setpoint serves as a predictor of AIDS and it gets stabilize after a period of acute HIV infection (Burg et al. 2009).
In vitro studies have shown that HIV-infected CD4+ T cells are lysed by HIV viral proteins and by mathematical modeling the primary HIV infection kinetics, it has been revealed that the viral-induced cytopathicity and susceptible T cells (also called the 'target-cell limited model') controls the infection (Somasundaran and Robinson 1987; Muller and Bonhoeffer 2003).
By observing the different behavioral patterns of HIV, TNF-\(\alpha\) and NF-\(\kappa\)B with respect to time, it can be observed that it depends upon time delays taken by these key players to incite or repress a specific response. These conditions help us to understand that how the equilibrium is maintained in the body upon arrival of any invader (or infection), e.g., if the qualitative period of HIV remains slower than TNF-\(\alpha\), no symptoms will appear until this condition is violated. In this way the immune cells (macrophages) can suppress the viral infected cells or hinder their binding to macrophages by speeding up the production of cytokines, chemokines and several other growth factors essential for immune stimulation. This model of HIV-1 infection can be further extended to explore disease behaviour to get deeper insights. The involvement of Rel Homology domain with discriminatory domains (Rel A and Rel B) having contrasting functions make NF-\(\kappa\)B a very interesting candidate for drug designing. The contrasting role of Rel B in anti-inflammation as compared to Rel A (causing inflammation) revealed interesting results in initial clinical findings. Experimental data shows that suppression of Rel A domain by introducing Rel B externally (or making it active by other means) tends to slow down the disease progression in the body and helps to inhibit inflammation caused by Rel A.
By adding more entities in the BRN understudy, exact therapeutic targets can be found by docking and simulation studies of the interacting proteins. These findings would aid in the process of drug designing. Model checking techniques could help us to design a computer aided tool for rapid immunological analysis. Real time modeling and its applications would pave the way towards predicting the preventive medicines, which will give the advancement in medicine developed for the treatment of AIDS.
Qualitative and hybrid models of the Biological Regulatory Network (BRN) of immune response against HIV-1 infection were constructed to explore the qualitative and real-time properties of the dynamics of the three entities HIV, \((NF \text{-} \kappa B)\) and \(TNF \text{-} \alpha\). The qualitative model predicted one stable state showing the over expression of all the interacting entities and two cycles representing a recovery mechanism of the BRN from HIV infection. The hybrid model was further analyzed to predict the delay constraints for the existence of the invariant cycle. Two important delay constraints were identified which show the oscillatory periods of HIV is greater than the oscillatory periods of \((NF \text{-} \kappa B)\) and \(TNF \text{-} \alpha\). These findings may help in the design of new generation drugs in the treatment of AIDS.
Douek DC, Brenchley JM, Ambrozak DR, Betts MR, Hill BJ, Okamoto Y, Casazza JP, Kuruppu J, Kunstman K, Wolinsky S, Grossman Z, Dybul M, Oxenius A, Price DA, Connors M, Koup RA (2002) HIV preferentially infects hiv-specific CD4+ t cells. Nature 417(6884):95–98
Legler DF, Micheau O, Doucey M-A, Tschopp J, Bron C (2003) Recruitment of tnf receptor 1 to lipid rafts is essential for tnf\(\alpha\)-mediated nf-\(\kappa\)b activation. Immunity 18(5):655–664
Vitale RFF, Ribeiro FDAQDEA (2007) The role of tumor necrosis factor-alpha (TNF-alpha) in bone resorption present in middle ear cholesteatoma. Braz J Otorhinolaryngol 73(1):117–121
Bingham CO 3rd (2002) The pathogenesis of rheumatoid arthritis: pivotal cytokines involved in bone degradation and inflammation. J Rheumatol 65:3–9
Herbein G, Montaner LJ, Gordon S (1996) Tumor necrosis factor alpha inhibits entry of human immunodeficiency virus type 1 into primary human macrophages: a selective role for the 75-kilodalton receptor. J Virol 70(11):7388–7397
Lane BR, Markovitz DM, Woodford NL, Rochford R, Strieter RM, Coffey MJ (1999) Tnf-alpha inhibits hiv-1 replication in peripheral blood monocytes and alveolar macrophages by inducing the production of rantes and decreasing c-c chemokine receptor 5 (ccr5) expression. J Immunol 163(7):3653–3661
Swami M (2013) HIV infections: Restricting HIV from macrophages. Nat Med 19(4):416. doi:10.1038/nm.3166
Herbein G, Gras G, Khan K, Abbas W (2010) Macrophage signaling in HIV-1 infection. Retrovirology 7(1):34. doi:10.1186/1742-4690-7-34
Melikyan GB (2008) Common principles and intermediates of viral protein-mediated fusion: the HIV-1 paradigm. Retrovirology 5(1):111
Herbein G, Keshav S, OLLIN M, Montaner L, Gordon S (1994) HIV-1 induces tumour necrosis factor and IL-1 gene expression in primary human macrophages independent of productive infection. Clin Exp Immunol 95(3):442–449
Ghosh S, Karin M (2002) Missing pieces in the nf-\(\kappa\)b puzzle. Cell 109(2):81–96
Grey ST (2008) Regulating inflammation: The ying and yang of nf-\(\kappa\)b activation. Immunol Cell Biol 86(4):299–300
Brandman O, Meyer T (2008) Feedback loops shape cellular signals in space and time. Science 322(5900):390–395
Thomas R, Thieffry D (1995) Dynamical behaviour of biological regulatory networks. I. Biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull Math Biol 57:247–276
Thomas R et al (1998) Laws for the dynamics of regulatory networks. Int J Dev Biol 42:479–485
Gagneur J, Casari G (2005) From molecular networks to qualitative cell behavior. FEBS Lett 579(8):1867–1871
De Jong H, Gouzé J-L, Hernandez C, Page M, Sari T, Geiselmann J (2004) Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bull Math Biol 66(2):301–340
Bibi Z, Ahmad J, Niazi U (2011) Dynamical modeling of the biological regulatory network of nf-kb activation in hiv-1. In: Computer Networks and Information Technology (ICCNIT), 2011 International Conference On, pp. 47—51. IEEE
Ahmad J, Bernot G, Comet J-P, Lime D, Roux O (2007) Hybrid modelling and dynamical analysis of gene regulatory networks with delays. ComPlexUs 3(4):231–251
Fromentin J, Eveillard D, Roux O (2010) Hybrid modeling of biological networks: mixing temporal and qualitative biological properties. BMC Syst Biol 4(1):79
Ahmad J, Niazi U, Mansoor S, Siddique U, Bibby J (2012) Formal modeling and analysis of the mal-associated biological regulatory network: insight into cerebral malaria. PLoS One 7(3):33532
Ahmad J, Bourdon J, Eveillard D, Fromentin J, Roux O, Sinoquet C (2009) Temporal constraints of a gene regulatory network: refining a qualitative simulation. Biosystems 98(3):149–159. doi:10.1016/j.biosystems.2009.05.002. (Evolving gene regulatory networks)
Ahmad J, Roux O (2008) Invariance kernel of biological regulatory networks. In: Proceedings of the 2008 IEEE International Conference on Bioinformatics and Biomedicine. BIBM '08, pp. 135–140. IEEE Computer Society, Washington, DC. doi:10.1109/BIBM.2008.64. http://dx.doi.org/10.1109/BIBM.2008.64
Thieffry, D, Thomas R (1998) Qualitative analysis of gene networks. Pacific Symposium on Biocomputing. pp 77–88
Thomas R, D'Ari R (1990) Biological feedback. Taylor & Francis, London. http://books.google.com.pk/books?id=WrRVk7B\_NLgC
Bernot G, Comet J-PP, Richard A, Guespin J (2004) Application of formal methods to biological regulatory networks: extending thomas' asynchronous logical approach with temporal logic. J Theor Biol 229(3):339–347. doi:10.1016/j.jtbi.2004.04.003
Clarke EM Jr, Grumberg O, Peled DA (1999) Model Checking. MIT Press, Cambridge
Khalis Z, Comet J-P, Richard A, Bernot G (2009) The SMBioNet method for discovering models of gene regulatory networks. Genes, genomes and genomics 3(1):15–22
Aslam B, Ahmad J, Ali A, Parachaa RZ, Tareen SHK, Niazi U, Saeed T (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
Ahmad J, Roux O, Bernot G, Comet J-P (2008) Analysing formal models of genetic regulatory networks with delays. Int J Bioinform Res Appl 4(3):240–262
Henzinger TA, Ho P-H, Wong-toi H (1997) Hytech: A model checker for hybrid systems. Softw Tools Technol Transfer 1:460–463
Gouze JL (1998) Positive and negative circuits in dynamical systems. J Biol Syst 6(01):11–15
Ahmad J (2009) Modélisation hybride et analyse des dynamiques des réseaux de régulations biologiques en tenant compte des délais
Bosque A, Planelles V (2009) Induction of HIV-1 latency and reactivation in primary memory CD4+ t cells. Blood 113(1):58–65
Wilson JA, Hevey M, Bakken R, Guest S, Bray M, Schmaljohn AL, Hart MK (2000) Epitopes involved in antibody-mediated protection from ebola virus. Science 287(5458):1664–6
Allen TM, O'Connor DH, Jing P, Dzuris JL, Mothé BR, Vogel TU, Dunphy E, Liebl ME, Emerson C, Wilson N et al (2000) Tat-specific cytotoxic T lymphocytes select for SIV escape variants during resolution of primary viraemia. Nature 407(6802):386–390
Davenport MP, Ribeiro RM, Perelson AS (2004) Kinetics of virus-specific CD8+ t cells and the control of human immunodeficiency virus infection. J Virol 78(18):10096–10103
Fauci AS (1996) Host factors and the pathogenesis of HIV-induced disease. Nature 384(6609):529–534
Kaufmann GR, Cunningham P, Kelleher AD, Zaunders J, Carr A, Vizzard J, Law M, Cooper DA et al (1998) Patterns of viral dynamics during primary human immunodeficiency virus type 1 infection. J Infect Dis 178(6):1812–1815
Lindbäck S, Karlsson AC, Mittler J, Blaxhult A, Carlsson M, Briheim G, Sönnerborg A, Gaines H, Group, K.I.P.H.I.S et al (2000) Viral dynamics in primary HIV-1 infection. Aids 14(15):2283–2291
Burg D, Rong L, Neumann AU, Dahari H (2009) Mathematical modeling of viral kinetics under immune control during primary HIV-1 infection. J Theor Biol 259(4):751–759
Somasundaran M, Robinson H (1987) A major mechanism of human immunodeficiency virus-induced cell killing does not involve cell fusion. J Virol 61(10):3114–3119
Muller V, Bonhoeffer S (2003) Mathematical approaches in the study of viral kinetics and drug resistance in HIV-1 infection. Curr Drug Targets Infect Disord 3(4):329–344
Conceived and designed the experiments:ZB JA. Performed the experiments: ZB. Analyzed the data: ZB JA AS SHKT. Contributed reagents/materials/analysis tools: ZB JA AS SHKT JA AA SS SK. Wrote the paper: ZB JA. All authors read and approved the final manuscript.
This research work was supported by Higher Education Commission (HEC) of Pakistan, National University of Sciences and Technology (NUST) and International Islamic University, Islamabad, Pakistan.
Research Center for Modeling and Simulation (RCMS), National University of Sciences and Technology (NUST), H13, Islamabad, 44000, Pakistan
Zurah Bibi
, Jamil Ahmad
, Amnah Siddiqa
& Samar HK Tareen
Atta-Ur-Rahman School of Applied Biosciences (ASAB), National University of Sciences and Technology (NUST), H13, Islamabad, 44000, Pakistan
Amjad Ali
& Hussnain Ahmed Janjua
Department of Bioinformatics and Biotechnology, International Islamic University, H12, Islamabad, 44000, Pakistan
Shaheen Shahzad
Department of Computer Science, University of Peshawar, Peshawar, 25120, Pakistan
Shah Khusro
Search for Zurah Bibi in:
Search for Jamil Ahmad in:
Search for Amjad Ali in:
Search for Amnah Siddiqa in:
Search for Shaheen Shahzad in:
Search for Samar HK Tareen in:
Search for Hussnain Ahmed Janjua in:
Search for Shah Khusro in:
Correspondence to Jamil Ahmad.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
TNF-\(\alpha\)
Qualitative modeling