Accepted Manuscript Title: On the Modelling and Analysis of the Regulatory Network ofDengue Virus Pathogenesis and Clearance Author: Babar Aslam Jamil Ahmad Amjad Ali Rehan Zafar Paracha Samar Hayat Khan Tareen Umar Niazi Tariq Saeed PII: DOI: Reference:
S1476-9271(14)00126-1 http://dx.doi.org/doi:10.1016/j.compbiolchem.2014.10.003 CBAC 6373
To appear in:
Computational Biology and Chemistry
Received date: Revised date: Accepted date:
27-5-2014 1-9-2014 6-10-2014
Please cite this article as: Babar Aslam, Jamil Ahmad, Amjad Ali, Rehan Zafar Paracha, Samar Hayat Khan Tareen, Umar Niazi, Tariq Saeed, On the Modelling and Analysis of the Regulatory Network ofDengue Virus Pathogenesis and Clearance, Computational Biology and Chemistry (2014), http://dx.doi.org/10.1016/j.compbiolchem.2014.10.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the Modelling and Analysis of the Regulatory Network of Dengue Virus Pathogenesis and Clearance Babar Aslama , Jamil Ahmadb,∗ , Amjad Alia , Rehan Zafar Parachaa , Samar Hayat Khan Tareenb , Umar Niazic , Tariq Saeedb a Atta-Ur-Rahman
cr
ip t
School of Applied Biosciences (ASAB), National University of Sciences and Technology (NUST), Islamabad, 44000, Pakistan. b Research Center for Modeling and Simulation (RCMS), National University of Sciences and Technology (NUST), Islamabad, 44000, Pakistan. c IBERS, Aberystwyth University, Edward Llwyd Building, Penglais Campus, Aberystwyth, Ceredigion Wales UK, SY23 3FG.
us
Abstract
Dengue virus can ignite both protective and pathogenic responses in human. The pathogenesis is related
an
with modified functioning of our immune system during infection. Pattern recognition receptors like Toll like receptor 3 is vital for the induction of innate immunity in case of Dengue infection. Toll like receptor 3 induces
M
TRIF mediated activation of Type 1 interferons and Fc receptor mediated induction of cytokines. Interferons have been related with clearance of Dengue virus but it has adopted modified regulatory mechanisms to counter this effect. SOCS protein is also induced due to the interferon and cytokine mediated signalling which
d
can subsequently play its part in the regulation of interferon and cytokine production. Our hypothesis in this
Ac ce pt e
study relates the pathogenesis of Dengue virus with the SOCS mediated inhibition of our innate immunity. We used the qualitative formalism of Ren´e Thomas to model the biological regulatory network of Toll like receptor 3 mediated signalling pathway in an association with pathogenesis of dengue. Logical parameters for the qualitative modelling were inferred using a model checking approach implemented in SMBioNet. A linear hybrid model, parametric linear hybrid automaton, was constructed to incorporate the activation and inhibition time delays in the qualitative model. The qualitative model captured all the possible expression dynamics of the proteins in the form of paths, some of which were observed as abstract cycles (representing homoeostasis) and diverging paths towards stable states. The analysis of the qualitative model highlighted the importance of SOCS protein in elevating propagation of dengue virus through inhibition of type 1 interferons. Detailed qualitative analysis of regulatory network endorses our hypothesis that elevated levels of cytokine subsequently induce SOCS expression which in turn results into the continuous down-regulation of Toll like receptor 3 and interferon. This may result into the Dengue pathogenesis during the stage of immunosuppression. Further analysis with HyTech (HYbrid TECHnology) tool provided us with the realtime constraints (delay constraints) of the proteins involved in the cyclic paths of the regulatory network 1
Page 1 of 41
backing the evidence provided by the qualitative analysis. The HyTech results also suggest that the role of SOCS is vital in homoeostasis Key words: Formal Modelling, TLR3, Dengue, Type 1 Interferons, SOCS-1, Ren´e Thomas, Kinetic Logic,
ip t
parametric linear hybrid automata, SMBioNet.
cr
1. Introduction
The incidences of dengue (DENV) fever have increased 30-fold in the last five decades [25]. Today
us
nearly 2.5 billion people, i.e. about 40 percent of the world population, is living in areas where there is a risk of DENV transmission. It has been reported as endemic in at least 100 countries [36]. According to a recent report of World Health Organization (WHO), around 50-100 million DENV infections occur per
an
5
annum including 500,000 cases of DENV haemorrhagic fever (DHF) and 22,000 cases of mortality worldwide [28, 62].
M
Dengue virus belongs to the genus Flavivirus, family Flaviviridae. It has four serotypes including DENV1, DENV-2, DENV-3 and DENV-4 [10, 61]. Transmission of the virus is mainly caused through the bites of 10
the female Aedes aegypti mosquito. Female Aedes aegypti acquires this virus from the blood of an infected
d
person [69, 54]. The virus infects the mid-gut of the mosquito and transfers to the salivary glands over
Ac ce pt e
a period of 8-12 days, which then is transmitted to humans on next bite [74]. The DENV virus after entering the systemic circulation replicates within the cell of mononuclear phagocyte lineage (macrophages, monocytes and B-cells). DENV infections may result in dengue fever (DF), dengue haemorrhagic fever 15
(DHF), or dengue shock syndrome (DSS) marked with high fever, skin rash, headache and severe muscular and joint pains [56, 74].
1.1. TLR3 Signalling Pathway
Recent studies implicates that the course of DENV infection has been attributed to the innate immunity, particularly generated by the pattern-recognition receptors (PRRs) such as Toll-like receptors (TLRs) [24, 20
96, 82, 93]. TLRs are the part of the innate immune system and activated by the recognition of pathogen associated molecular patterns (PAMPs) [7, 47]. The study of Iwasaki and Medzhitov on human monocytic cells U937 strongly supports the role of TLR3 for inducing immune response against DENV [46]. Figure 1 ∗ Corresponding
author Email addresses: [email protected] (Jamil Ahmad)
Preprint submitted to Computational Biology and Chemistry
October 13, 2014
Page 2 of 41
illustrates the possible signalling events generated by TLR3 during DENV infection. As key players of innate immunity, the recognition of several viruses via TLR is very important for the induction of various 25
cytokines to avert pathogenesis [46]. Several pro-inflammatory cytokines are induced as a result of the TLR3 activation and provide an interface between the innate and adaptive immune responses [43]. Modulation
ip t
of DENV infection by type 1 interferon alpha/beta (IFNα/β) has been suggested by observing DENVinfected patients and IFNα/β-receptor deficient mice [9]. In vitro analysis with several cell types and viral
30
cr
strains proved that DENV viral RNA decreased by multiple folds upon introduction of IFNα/β, however, IFN gamma (γ) showed variable responses depending on the cell type and signalling pathway [27]. At the
us
early stage of viral infection, inflammatory response through PRRs provide dominant role in priming and regulating innate immune system [55]. Yi-Ting Sai et al. in their study demonstrated the role of TLR3 in mediating strong IFNα/β induction and subsequent inhibition of DENV viral replication, which limits its
35
an
cytopathic effect [93, 27]. As Type 1 interferons upon their interaction with IFN receptors elicit antiviral, anti-proliferative and immunomodulatory roles. Shresta et al. demonstrated critical and non-overlapping
M
role of IFNα/β in resolving primary DENV infection in genetically deficient mice [78].
Ac ce pt e
d
Figure 1: Pathway of TLR3. DENV is sensed by TLR3 in endosomes [8]. TLR3 then activates TRIF which initiates the up regulation process of NFκB and IFNα/β by two distinct pathways [44]. Type 1 IFN pathway TRIF via NAP 1 and TRAF3 interacts with TBK1 which transduces a signal to Ikki which leads to phosphorylation of IRF [33]. Phosphotyrosine residues on TLR3 binds to P13K which is required for complete IRF activation [77]. Dimerized IRFs then translocated into nucleus where they bind to DNA and in conjugation with the co-activator CBP/p300 they transactivate type 1 IFN gene expression [63]. Pro inflammatory cytokine pathway. NFκB up regulation by TLR3 is mediated by TRAF6 interaction with TRIF that involves the TAK1/TAB2/TAB3 complex which results in activation of TAK1 [83]. TAK1 then phosphorylate and activate IKK. Then NFκB nuclear localization sequence is unmasked by degradation effect of phosphorylated IKB. Activation of NFκB and transcriptional induction of target genes requires signalling from TRIF via RIP1 [23]. NFκB up regulation via TRIF, TRAF 6 results in up regulation of IL-10 which leads to expression of SOCS via JAK / STAT pathway [65, 50]. SOCS can inhibit this signalling pathway at four sites including degradation of NFκB [57] , TLR3 [92] , IFNα/β and JAK / STAT. IFNα/β have been seen in regulation of DENV whereas NS4B protein of DENV inhibit the role of IFNα/β [60]. IFNα/β can also lead to induction of SOCS [67]. DENV by its ADE (Antibody dependent enhancement) mechanism through FcR signalling [16, 73, 19] up regulates IL-10 production [32, 17].
Expression of DENV non-structural (NS) proteins including NS2A, NS4A or NS4B in human cells support the elevated replication of an IFN sensitive virus [60]. Cells infested with NS4B protein or DENV do not exhibit IFN-α/β mediated induction of nuclear signal transducer and activator of transcription (STAT)1. 40
By blocking STAT1 activities, NS4B protein of DENV may interfere with the IFNα/β mediated antiviral response [60]. These findings indicate that NS4B might be involved in blocking IFNα/β signalling during DENV infections. IFNα/β activate STAT1–STAT3 heterodimers [66], followed by SOCS1 and SOCS3 mediated downregulation [80]. Qin et al. implicated in their study about the up-regulation of IFN mediated induction 3
Page 3 of 41
45
of SOCS-1, which in return down-regulate IFN signalling and reduce the expression of CD40 cells [67, 80]. They also proposed the inhibitory mechanism of JAK/STAT phosphorylation which ultimately inhibit the vital signalling cascade for the induction of cytokines [67]. Various studies have established the role of Antibody-dependent enhancement (ADE) of dengue virus
50
ip t
pathogenesis [53]. The pre-existing antibodies (Abs) are major contributors for development of DHF and DSS and fail in neutralizing the infecting serotype of DENV [38]. The Fc receptors (FcR) family offers a major
cr
example of how to promptly trigger, activation and inhibition of the signalling pathway to maintain harmony between innate and adaptive immune response [70]. In accumulation to their exact roles in prompting the
us
activation of innate effector cells, FcγRs also play a role in the regulation of B-cell activation, antigen presentation and immune-complex mediated maturation of the dendritic cells (DCs) [13]. A new hypothesis 55
termed as intrinsic ADE, suggests that DENV internalization via FcγR results in the induction of interleukin-
an
10 (IL-10) production, biased T-helper-1 (Th)1 to Th2 response and suppression of innate immunity, which leads to both high levels of antibodies (Abs) and viral load in dengue patients [94, 39]. Fujii et al. in a recent
M
study, suggested the up-regulation of IL-10 levels via FcγR [32]. As a consequence of ADE, complementary mechanism to higher viremia may be explained as FcγR mediated DENV entry due to suppressed antiviral 60
immune response. In the case of DENV infection, suppressed levels of IL-12, IFNγ, TNFα and nitric
Ac ce pt e
environment for viral replication [17].
d
oxide along with the elevated levels of anti-inflammatory cytokines tend to provide immunosuppressive
1.2. Computational Systems Biology
Computational systems biology provides a compelling platform to answer critical scientific questions 65
about signalling mechanisms based on pragmatic modelling and theoretical elucidation [51]. Multitude of biological facts have been uncovered by molecular biology including functional properties of the genes and proteins. However, for the interpretation of the dynamics of the biological system, these informations require further investigation [51]. To understand the intrinsic complexity of biological systems experimental and computational approaches in combination are speculated to solve this problem [52, 45]. Various com-
70
putational approaches assisted with the experimentally observed phenomenon in previous literature helped us to understand the behaviours of the complex biological regulatory networks involved in DENV replication and immune evasion [45]. A model that captures biological observations related to physiological and pathological conditions could also predict the underlying dynamics behind these observations: oscillations or cycles, stable states and bifurcation points. Interestingly, the bifurcation points (specific pattern of protein
75
interactions), may play a role in leading the system towards pathogenic state, can provide us an insight 4
Page 4 of 41
about the useful therapeutic target. Thus, instead of targeting proteins in the pathogenic state (abnormal physiological state) drugs may target entities involved in the interactions causing the physiological abnormalities. We predicted with the help of a qualitative model the underlying causes of the hyper-inflammatory condition in the case of cerebral malaria (CM) pathogenesis [4]. The model predicted that longer activation and interactions of two proteins BTK and MAL, and the pre-existing high levels of cytokines play key role
ip t
80
in the pathogenesis of CM. The later hypothesis was validated with a study conducted by Franklin et al.,
cr
2011 [31, 49, 40, 30, 68] which emphasizes the relationship between elevated levels of inflammatory cytokines and CM. Our recent study [64] shows a strong affinity of BTK with MAL protein. These studies suggest
85
us
that the expression dynamics of the pathway should be studied in order to have an understanding of the behaviours of the systems and to predict drug targets which would ultimately increase our confidence about the potential drugs and would accelerate the drug design process.
an
For the last few decades, formal methods have received the attention of researchers in various disciplines to formally model and verify complex systems. Model checking is an automatic formal verification technique
90
M
which exhaustively explore dynamics of a system formally expressed as automata (M) to check whether it meets specifications formally expressed as formulas (Φ) [21]. Formally, a model checking process is represented as M |= Φ. Verification could be carried out even if the system is parametric and in this case
d
the values of the parameters are synthesised for the properties to be verified [81, 42]. Thus, parameter
Ac ce pt e
synthesis is another key feature of model checking.
In the current study, TLR3 mediated signalling pathway has been studied using the generalized logical 95
approach introduced by Ren´e Thomas [87, 90]. This approach represents the qualitative analysis of dynamic behaviour of the system with a focus on the feedback mechanisms. This approach is successfully applied on some other signalling systems including the MyD88-adapter-like (MAL) associated biological regulatory network (BRN) [4], control of differentiation process in helper T cells [58], control of organ differentiation in Arabidopsis thaliana flowers [59] and segmentation during embryo-
100
genesis in Drosophila melanogaster [84, 75]. Kinetic formalism of Ren´e Thomas and its implementation in GENOTECH [4] provided us an in silico approach to test hypotheses related to DENV pathogenesis. Ren´e Thomas’ logic relies on the discretization of expression levels which can be determined by the qualitative thresholds. The evolution of these levels are driven by the logical parameters. Values of these parameters can be inferred with the help of SMBioNet [48], which implies a model checking approach by first expressing
105
the observations (properties) in Computation Tree Logic (CTL) formulas and then computes all possible sets of parameters leading to models satisfying these observations. Once the reliable set of parameters 5
Page 5 of 41
are determined, the tool GENOTECH then computes the qualitative model showing all the behaviours of the network and where important behaviours can be seen as paths (trajectories), closed paths (cycles) and stable states (see state graph in Section 3). Delays are incorporated in the qualitative model to model the 110
activation and inhibition orders of proteins; which protein activates or deactivates first now depends on their
ip t
time delays. The values of such parameters can be synthesised using the hybrid model checking tool HyTech [41].
cr
1.3. Our Contribution
In this study, we performed a formal analysis of TLR3 signalling pathway using its discrete and hybrid models. In these models, we analysed the role of SOCS protein in the context of TLR3 mediated IFNα/β
us
115
production in DENV infection. The design of this study is graphically depicted in Section 2. First, the
an
Biological Regulatory Network (BRN) is obtained (Figure 2) by keeping the prominent entities (in which we are interested) of the pathway and abstracting the rest, provided that the behaviour is preserved. For +
−
→) ‘B’ which inhibits(− →) ‘C’, and the proteins ‘A’ and ‘C’ are the example, if an entity ‘A’ activates(− −
→ C while implicitly showing the behaviour ones in which we are interested, we simply abstract it to A −
M
120
that ‘A’ inhibits ‘C’ through the activation of ‘B’. This BRN is then used to obtain the logical parame-
d
ters from SMBioNet. These parameters are used in GENOTECH which computed the qualitative cycles,
Ac ce pt e
a stable state, and transformed it into parametric linear hybrid automata in HyTech format. Cycles were showing the oscillations of all proteins while the stable state was showing the over-expression of SOCS along 125
with DENV. The delay constraints of a significant cycle, representing its invariance kernel (a set of cyclic trajectories which start and end at the same point, see Definition 12, [5]), were computed by the analysis of parametric linear hybrid automaton (hybrid model) of TLR3 signalling pathway in HyTech tool. The computation of the invariance kernel enabled us to observe the sustained cyclic nature of this cycle, as well as the possible delay values in the form of linear constraints. Our study implicates that when SOCS was
130
oscillating within the qualitative cycles then all the other proteins of the pathway including SOCS itself were showing sustained oscillations in the hybrid model. Our analysis suggests that altered oscillatory behaviours of SOCS protein play a vital role in the inhibition of TLR3 mediated IFNα/β production, which ultimately results in the down-regulated immune response with subsequent elevation in DENV titre.
6
Page 6 of 41
Figure 3: Work flow diagram of the design of the study.
2. Methods
cr
135
ip t
Figure 2: TLR3 associated BRN. Reduced form of the TLR3 induced BRN shown in Figure 1. Entities are represented by vertices whereas edges represent interaction between them. An association between two corresponding entities is represented by each edge. Interaction types are represented by sign on arc i.e. negative for inhibition (dotted arrows) and positive for activation (solid arrows). Threshold levels of entities which are necessary to activate or inhibit its successor entity (see Definition 2 in Method section and Table 2) are represented by given values on the arc.
us
Usually, biological systems are represented by conventional models having differential equations referring to the time derivative of different quantities such as concentrations, temperatures and rates. As a result of inherent complexity of biological systems, little is known exactly about these entities. In the year 1970, Ren´e
140
an
Thomas [87] presented a Boolean logic based method to model biological regulatory networks. He presented the efficacy of his modelling methodology by checking the lambda phage genetic switch [87, 59, 58]. Ren´e
M
Thomas noted the limitations in Boolean logic as it has only two levels, 0 or 1, that are not able to formalize related problems. So, Thomas generalized the Boolean logic to kinetic logic and proved the efficacy of the methodology by utilizing it in different genetic regulatory systems [85, 88, 87, 79, 90, 91, 26, 89], and found
145
Ac ce pt e
d
it to be efficient in approximating the differential equations based modelling [88]. 2.1. The Kinetic Logic formalism of Ren´e Thomas Kinetic logic formalism models biological regulatory networks (BRNs) based on the qualitative threshold levels of entities and their types of interactions [88]. Value of at least 1 as minimum threshold for any entity is required to bring a positive or negative change in the target entity [86]. Entities are representing proteins, genes or both of them that can interact with each other to bring positive or negative effects showing 150
stimulatory or inhibitory effects on the system, respectively. Positive effect shows that the expression level of an entity increases the rate of activation or production of other entity, whereas, negative effect reduces the rate of activation. Biological activation and inactivation are sigmoidal (non-linear bounded curve) in nature and are discussed below. The following definitions and concepts define the kinetic logic formalism and are mainly adopted from [4].
155
2.1.1. Activation If an entity e1 enhances the rate of synthesis or activation of another entity e2 then e1 is the positive regulator for e2 . In this case, every time the concentration of e1 enhances or reduces, it will affect the 7
Page 7 of 41
concentration of e2 in similar manner. This relationship can be shown by sigmoidal nature. The sigmoidal curve can be shown by its step function and is shown in Figure 4 A. We can see that below a threshold 160
level θ, e1 slightly changes the rate of activation of e2 . But as soon as the concentration of e1 reaches its threshold level θ, the rate of activation of e2 elevates rapidly to a level where it saturates its target entity
ip t
e2 . So, we can say that when e1 < θ it can not bring any positive change and if e1 ≥ θ then e1 can bring
cr
the positive change in e2 .
165
an
us
Figure 4: Dummy tendency graph of biological regulation. (A)Positive regulation of a target entity e2 by an activator e1 is presented in upper left corner of sigmoidal graph , whereas upper right one is its step function.(B) Negative regulation of a target entity e4 by an inhibitor e3 is presented in lower left corner sigmoidal graph, whereas right one is its step function. There is no or little effect on its target entity when the regulator remains below threshold level (represented by θ on X-axis). The rate of activation or inactivation of target entity leads it to its saturated level when regulator reaches to its threshold concentration level.
2.1.2. Inhibition
M
If an entity e3 decreases the rate of activation of another entity e4 , then e3 can be called as negative regulator or inhibitor of e4 and the process is called as negative regulation or inhibition. The process of
d
inhibition also depends on the threshold levels of e3 showing that whenever the concentration level of e3
170
Ac ce pt e
reaches its threshold level then e3 will be efficient in decreasing the activation levels of e4 , up to a level where e4 can be completely inactivated or its affect in the system is abolished. The affect of inhibition is also of sigmoidal nature and is shown in Figure 4 B along with its step function. 2.2. Semantics of the Kinetic Logic Formalism
This part shows the semantics of Thomas formalism with a toy example of a BRN comprising of three entities (e1 , e2 and e3 ) shown in Figure 5. 175
Figure 5: Dummy BRN. An example of BRN in which interaction between genes e1 , e2 and e3 are labelled with + sign representing activation and integer ”1” representing the threshold levels.
Definition 1 (Directed Graph). A directed graph G(N, E) is a tuple, where: set of all nodes is N and set of ordered pair of edges or arcs is denoted by E ⊆ N × N 8
Page 8 of 41
(e1 , e2 ) is an edge which is considered to be directed from e1 to e2 , where e1 is called the head and e2 is 180
called the tail. In G(N, E), G− (e) and G+ (e) denote sets of predecessors and successors, respectively, of a node e ∈ N . N = {e1 , e2 , e3 } and E = {(e1 , e2 ), (e2 , e3 ), (e3 , e1 )} are sets of nodes and edges, respectively,
and G+ (e1 ) = {e2 } represents the predecessors and successors of e1 , respectively.
Definition 2 (Biological Regulatory Network). A biological regulatory network (BRN) is a labelled directed graph G(N, E), where the set of nodes N represents biological entities and interactions between entities are represented by the set of edges E ⊆ N × N , in which each edge ei → ej is labelled by a couple (jei ej , ηei ej ), jei ej being a positive integer representing qualitative threshold and ηei ej being either (+) or (-) sign representing activation or inactivation, respectively. Moreover, we define:
190
us
cr
185
ip t
which are shown as network of genes in Figure 5 . In case of gene e1 , (set of regulators of e1 ) G− (e1 ) = {e3 }
The outdegree of ei is the total number of targets of ei , it is denoted by pei and it is such that ∀ek ∈ G+ (ei ) jei ,ek pei .
an
Set Zei = {0, 1, ...., rei } of an entity ei represents its qualitative levels (concentration levels) where rei = max{jei ,ek |ek ∈ G + (ei )} To capture the behaviours of a BRN, it is important to define the set of all possible number of qualitative states.
M
195
d
Definition 3 (Qualitative States). A qualitative state s ∈ S of a BRN G(N, E) is a tuple which denotes a configuration of the qualitative expression levels of all entities. The set of states is :
200
Ac ce pt e
S = Πei ∈N Zei
Vector (sxei )∀ei ∈N represents the qualitative states, where level of expression (concentration) of a product ei is denoted by xei .
One associates with each qualitative state, the set of resources of each entity which corresponds to the set of predecessors which effectively control the entity (either a present activator or an absent inhibitor).
205
Definition 4 (Resources). At level xei , the set of resources Qxei of an entity ei ∈ N is defined as Qxei = {ej ∈ G− (ei )|(xej ≥ jej ei and ηej ei = +) or (xej < jej ei and ηej ei = −)}. In the above definition, it is to be noted that the inhibitor is considered as an activator in its absence.
210
Definition 5 (Logical Parameters). Logical parameters influence the dynamic behaviour of a BRN. The set of logical parameters are defined as K(G) = {Kei (Qxei ) ∈ Zei ∀ei ∈ N }. Logical parameters and set of resources are shown in Table 1 for the running example of three entities (Figure 5). To define the transition between two qualitative states, it is convenient to define the evolution operator [11]: 9
Page 9 of 41
Definition 6 (Evolution Operator). The evolution operator is defined as, ⎧ ⎨ xei + 1 if xei < Kei (Qxei ); xe − 1 if xei > Kei (Qxei ); xei Kei (Qxei ) = ⎩ i if xei = Kei (Qxei ). xei
cr
where xei and Kei (Qxei ) ∈ {0} ∪ N
ip t
Figure 6: State graph of the dummy BRN. State graph of BRN presented in Figure 5. Qualitative (discrete) states are presented as nodes. Qualitative levels of entities e1 , e2 and e3 , respectively, are shown as values inside a node. The state (0,0,0), and (1,1,1) represent deadlocks, whereas as transition from the state (1,0,1) to (1,0,0) creates a cyclic trajectory.
us
Table 1: Table of states, resources and logical parameters
xe2
xe3
Qxe1
Qxe2
Qxe3
Ke1 (Qxe1 )
Ke2 (xxe2 )
Ke3 (Qxe3 )
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1
{} {e3 } {} {e3 } {} {e3 } {} {e3 }
{} {} {} {} {e1 } {e1 } {e1 } {e1 }
{} {} {e2 } {e2 } {} {} {e2 } {e2 }
0 1 0 1 0 1 0 1
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
d
M
an
xe1
215
220
Ac ce pt e
All the states (xe ), set of resources (Qxe ) and logical parameters Ke (Qxe ) of the dummy BRN (Figure 5)are listed.
Definition 7 (State Graph). Let G(N, E) be a BRN, let s = (sxei )i ∈ S be a state and let consider the family of logical parameters K(G). The state graph of the model is the directed graph noted G = (S, T ) where S is the set of states and T ⊆ S × S is the relation between states called transition relation such that s → s ∈ T if and only if: 1. ∃ a unique ei N such that sxei = sxe and sxe = sxei Kei (Qxei ), and i i 2. ∀ ej N \ {ei } sxe = sxej . j
Table 1 represents the state table of running example which shows all the possible states, set of resources and logical parameters. This table helps to construct the state graph which captures all the behaviours of the running example. Two stable states namely (000) and (111) can be observed clearly in the state state graph shown in Figure 6. 225
2.3. Discrete Modelling of DENV associated BRN using GENOTECH GENOTECH [4] constructs a qualitative model of BRNs according to Thomas’ formalism and produce the state graph, which shows the dynamical behaviour of a BRN. Modelling is facilitated by GUIs, which 10
Page 10 of 41
allow to construct a BRN as directed labelled graph. The tool can perform the basic analyses of the state graph like the identification of stable states, cycles and path finding between two states. The state graph 230
can be saved in DOT format [34] for further study in the Graphviz tool [29].
ip t
2.4. Model Checking Although Ren´e Thomas framework captures the qualitative dynamics of biological regulatory networks, its semantics are derived from the values of logical parameters. Once the state graph is generated by assuming
235
cr
a certain value of these parameters, it is still required to manually verify the model against some known observations. This manual process is quite time expensive even for small BRNs containing only two entities.
us
Model Checking is one of the formal techniques in computer science that verifies models against specification formally expressed as temporal logic formulas. Bernot et al.[12] applied model checking to reverse engineer
an
the values of logical parameters from biological observations specified as CTL (Computational Tree Logic) formulas. This technique is based on exhaustive enumeration of parameters values, each time a parameter 240
combination is used to generate the state graph. Finally, by invoking a model checker, the verification
M
of CTL formulas is carried out and all parameters that are coherent with CTL formulas are selected as set of acceptable parameters. Thus model checking solves this problem by not only expressing biological
245
Ac ce pt e
2.4.1. Quantifiers
d
behaviours in a formal temporal language but also reverse engineering the values of logical parameters.
A path is a sequence of states in a state graph. In Computational Tree Logic (CTL), properties are checked in the state graph by applying quantifiers to the paths originating from a given state, as well as the states occurring on each of the originating paths. Description of each of the quantifiers is given below: A : This is a path quantifier which enforces that a given property should hold in all paths originating from the given state. The quantifier itself is read as “For all paths”.
250
E : Known as “Existential Quantifier”, this is also a path quantifier which enforces that a given property must hold in at least one path originating from the given state. The quantifier is read as “There exists a path”. G : This quantifier is known as the “Global Quantifier” and is a state quantifier which enforces that a property holds in all states of a path originating from the given state, inclusive of the given state as
255
well. It is read as “globally”. 11
Page 11 of 41
F : The “Future Quantifier” is the second state quantifier and enforces that a given property must hold in one of the future states in the path originating from the given state. It is read as “in future” or “eventually”. The future quantifiers also covers the current/given state as well when checking the property. X : The “Next Quantifier” is the third state quantifier and enforces that a given property must hold
ip t
260
in the immediate successor state. It is read as “next”.
cr
The usual format for a CTL formula is φ = {A/E}{G/F/X}{ψ/φ} where ψ is a predicate representing
us
a property and φ another CTL formula. Then, this format allows nesting of quantifiers. 2.4.2. SMBioNet
SMBioNet [12, 71, 48, 72] is a tool for the parameter enumeration of Biological Regulatory Networks
an
265
(BRNs), based on the Discrete Framework of Ren´e Thomas [87]. Given a model of BRN in the form of Thomas network and behavioural properties (observations), expressed as CTL formulas, SMBioNet exhaus-
M
tively enumerates all compatible parametrizations by generating a state graph for each parameter combination and by verifying the formulas on each state graph. The verification of CTL property is performed by invoking model checker NuSMV [20]. The parameter combinations are reduced by applying Snoussi and
d
270
observability constraints [15]. Finally, all the models that satisfy CTL properties are short-listed. SM-
Ac ce pt e
BioNet has been applied in studies like tail resorption in tadpole metamorphosis [48] and immunity control in bacteriophage lambda [71].
2.5. Hybrid Modelling to Enhance Qualitative Modelling with Delays 275
Hybrid modelling merges the discrete modelling with the continuous domain of time via special variables called ’clocks’. Likewise, the transitions are characterised as discrete and continuous where, Discrete transitions represent instantaneous switching of locations or states, and Continuous transitions represent elapsing of time within the same location or state. The natural environment of the cell is difficult to accurately model through discrete modelling frame-
280
works. Take the concentration of a protein for instance – the discrete model represent their expression levels as positive integers where a protein is at ‘0’ level at one instance, and immediately jumps to ‘1’ level the very next instance. This does not accurately replicate the continuous process of increasing expressions that proteins go through in the natural cellular environment. Ahmad et al. [6] proposed the use of piecewise 12
Page 12 of 41
linear curves for the modelling of the sigmoidal nature of protein expression (Figure 7), as opposed to the 285
step functions used in discrete models.
ip t
Figure 7: Activation and Inhibition delays. The clock ha 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+/− .
Clock variables are used to measure the ‘delays’ (the time duration) that needs to pass between two
cr
consecutive expression levels. Thus a clock variable h is associated with each protein in the BRN. The initial values of each h are set to zero, which then approach either d+ or d− . d+ signifies production delay,
290
us
that is, the delay required to increase the associated proteins concentration level by 1. Similarly, d− signifies the degradation delay, that is, the delay to decrease the said concentration by a single level. The rate of evolution of each h is given by the first order derivative dh/dt = r where r ∈ {0, 1, −1} [6, 2].
an
In most cases, the exact values of the delays associated with the proteins are not known which is why unvalued parametric delays are instead used [4]. Thus, the hybrid model was constructed using Parametric
295
M
Bio Linear Hybrid Automaton [1] defined below and given as Supplementary File 3. Let C = (X, P ), C ≤ (X, P ), and C ≥ (X, P ) be the set of constraints using only =,≤, and ≥, respectively. Here, X and P are the sets of real valued variables and parameters, respectively.
Ac ce pt e
d
Definition 8 (Parametric Bio Linear Hybrid Automaton (Bio-LHA). A parametric Bio linear hybrid automaton B is a tuple (L, l0 , X, P, E, Inv, Dif ) where, L is a finite set of locations, 300
l0 ∈ L is the initial location,
P is a finite set of parameters (delays),
X is a finite set of real-valued variable (clocks), E ⊆ L×C = (X, P )×2X ×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 and the reset set R ⊆ X, 305
Inv : L → C ≤ (X, P ) ∪ C ≥ (X, P ) 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 T, where T∗ = T \ {0}.
310
Figure 8 represents the conversion of the discrete state graph of Figure 6 into a Bio-LHA. Definition 9 (Semantics of Bio-LHA). Let γ be a valuation for the parameters P and ν represents the values of clocks in a location. The (T, γ)-semantics of a parametric Bio-LHA B = (L, 0 , X, P, E, Inv , Dif ) is defined as a timed transition system B = (S, s0 , T, →) where: (1) S = {( , ν) | ∈ L and ν |= Inv ( )}; (2) s0 is the initial state and (3) the relation →⊆ S × T × S is defined for t ∈ T as: 13
Page 13 of 41
315
ip t
Figure 8: The Parametric Bio Linear Hybrid automaton of the running example of the dummy BRN based on its state graph represented in Figure 2. Each location is marked with its configuration (state) name (for example (1,0,0)), the invariants (conjuncted constraints), and the entity rates. Each transition is marked with a guard and clock reset of the respective changing entity. Location (1,0,0) is the initial location, represented as a small arrow on its upper-left side. Locations (0,0,0) and (1,1,1) are representing deadlock states as: (i) there are no outgoing transitions from these locations and (ii) all rates are equal to zero. The rates of the entities in each location are based on the evolution in the first and second successor locations. 0
discrete transitions: ( , ν) −−→ ( , ν ) iff ∃( , g, R, ) ∈ E such that g(ν) = true, ν (h) = 0 if h ∈ R and ν (h) = ν(h) if h ∈ / R.
cr
t
us
continuous transitions: For t ∈ T∗ , ( , ν) −−→ ( , ν ) iff = , ν (h) = ν(h) + Dif ( , h) × t, and for every t ∈ [0, t], (ν(h) + Dif ( , h) × t ) |= Inv ( ), where |= represents satisfaction operator. Using the semantics of the Bio-LHA, Ahmad et al. [2] were then able to define the temporal state space and the invariance kernel set which have been adapted below.
an
320
325
M
Definition 10 (Temporal Zone). Temporal zone is defined as a region where time elapses until a discrete transition between states takes place. Definition 11 (Temporal State Space). The temporal state of a BRN is composed of the complete set of temporal zones derived from the discrete model of the said BRN.
d
In the hybrid model of the BRN, the sequence of points in a trajectory is denoted by φ(t) for t ∈ R ≥ 0 [3].
Ac ce pt e
A particular trajectory is said to be viable according to viability domain S, iff it remains within a prescribed region S [3]. 330
Definition 12 (Invariance kernel). A trajectory φ(t) is said to be viable in S if φ(t) ∈ S ∀ t ≥ 0. A subset K of S is an invariant if for any point p ∈ K, a trajectory starting in p is viable in K. An invariance kernel K is the largest invariant subset of S.
3. Results
3.1. Logical Parameters
The BRN given in Figure 2 comprises four entities; namely DENV, IFNα/β, TLR3 and SOCS. DENV 335
activates TLR3 and SOCS, which then formulates a positive feedback loop with IFNα/β due to the even number of negative interactions involved in this loop. The presence of SOCS inhibits the synthesis of IFNα/β and TLR3, while its own production is positively regulated by IFNα/β. Thresholds fixed by variables TLR3, IFNα/β, SOCS and DENV lead to a total of 16 possible states in the state graph. The state of the system at any given time is represented by a vector (TLR3, IFNα/β, SOCS, DENV). Before the infestations of DENV,
340
all the proteins that participate in this network are absent. Hence, (TLR3=0 & IFNα/β=0 & SOCS=0 14
Page 14 of 41
& DENV=0) constitutes a starting state of the system. The presence of two responses within the same network, derived from the literature, implies that there exist two possible paths starting from the initial state i.e. from all zeros state, can proceed towards DENV clearance and a pathogenic state [95]. The DENV clearance or normal state is known to be characterized by low titres of DENV, whereas the pathogenic state is characterized by its high titre and elevated levels of SOCS [97]. These two observations are encoded in
ip t
345
cr
CTL formulas as described in following equations.
P athogenic = (SOCS = 1 & DEN V = 1) Φα = Init ⇒ EF (AG(pathogenic))
an
Φβ = Init ⇒ AX(EF (Init))
us
Init = (T LR3 = 0 & IF N α/β = 0 & SOCS = 0 & DEN V = 0)
(2) (3) (4) (5)
M
Φ = Φα ∧ Φβ
(1)
The desirable qualitative model of DENV should contain closed or cyclic paths, which originate from an initial state, represented by all zeros state, and then return back to its original position from where
diseased state can be achieved if homoeostasis could not be maintained. Along with the presence of closed
Ac ce pt e
350
d
it started. Closed or cyclic paths actually represent the homoeostasis of the immune system, however, a
paths or cycles in the dynamics of DENV infection, pathogenesis can be encountered if dynamics proceed towards the stable steady state (0,0,1,1). The first formula, denoted as Φα models the behaviour of the BRN that reflects that pathogenic state is reachable from an initial state and remains as its future state. If the future state persists for a long time then we can call this as stable state. 355
The second property, given as Φβ expresses that starting from an initial state, the immune system is able to reset its configurations by going back to its original position and maintain the immune homoeostasis. The formula Φ combines (conjuncts) these two properties and is used for the selection of models by SMBioNet. SMBioNet selects 43 sets of parameters (each set represents a specific qualitative model)(Supplementary file 1) that satisfy the above mentioned CTL formula. These qualitative models were then further analysed
360
in GENOTECH by generating state graphs, which allowed us to analyse the expression patterns of the proteins, cyclic paths, stable states and the dynamics of the BRN. The selected set of parameters given in Table 2 were used to model the BRN Figure 2, which gives a state graph shown in Figure 9. This model was selected since it incorporates the above mentioned observations in the form of a stable state (0,0,1,1) 15
Page 15 of 41
and homoeostasis, represented by the qualitative cycles visible in Figure 9. In particular, we are interested 365
in the cycle which follows well ordered regulation pattern evident in the TLR3 pathway (Figure 1) where the occurrence of DENV activates TLR3 which in turn activates IFN before the activation of SOCS. The activation of IFN suppresses DENV resulting in the down regulation of TLR3 and in the up regulation of
ip t
SOCS. Finally, IFN is down regulated by the presence of SOCS, which is then degraded in the absence of its activators.
This biologically plausible pattern is found in 9 models out of 43, and is represented as (0,0,0,0) →
cr
370
(0,0,0,1) → (1,0,0,1) → (1,1,0,1) → (1,1,0,0) → (0,1,0,0) → (0,1,1,0) → (0,0,1,0) → (0,0,0,0). The detailed
us
method for the selection of a single parameter set (given in Table 2) from the 9 sets is provided in Supplementary file 5. It is important to note that the change in the parameters (that is, if any parametrisation from the other enumerated sets is selected) will change the qualitative dynamics of the BRN. This in turn will change the hybrid model (see section 3.3), affecting the delay constraints themselves, as well as their
an
375
M
biological implications and significance.
Ac ce pt e
d
Figure 9: State graph of physiological signalling of TLR3. The states, represented in oval shapes are named by the qualitative expression levels of each entity in the order TLR3, IFNα/β, SOCS, and DENV. The normal stable state (0,0,0,0) highlighted in blue is taken as initial state. The recommended diseased state after pathogenesis of DENV is (0,0,1,1) which is presented in red colour. Transitions from one state to the other are shown by different types of edges. Bold edges (condition 1) represent the transitions which may eventually lead to the clearance of DENV infection. Short dashed edges (condition 2) represent the transitions which make diverge the system towards the SOCS and DENV mediated pathogenic state. Dotted edges (condition 3) reflect that some states can converge towards the cyclic behaviour of the system.
3.2. Qualitative Model
This section presents the results obtained by the modelling and analysis in GENOTECH, which simulate the model according to Ren´e Thomas’ formalism. This tool is used to model BRN along with the usage of 380
set of logical parameters obtained from SMBioNet. As output, this gives a state graph (Figure 9) showing the paths between any two states (transitions), stable states and possible cycles among different states. The state graph represents all possible transitions from one state to another, where each state shows the expression levels of every entity at a particular time. From a stable state, which is like a sink, the system does not evolve to other states so has been termed as stable (pathogenic) state.
385
The state (0,0,0,0) represents that all entities are at zero level while (0,0,0,1) shows the presence of DENV. These two states are included in all the cycles of the BRN as shown by solid arrows (Figure 9). The tendencies of divergence from cycles towards stable state (0,0,1,1) are also shown by dashed edges. The cycles in the state graph represent the homoeostasis of the BRN to overcome or control the level of the viral 16
Page 16 of 41
Table 2: Possible parameter values for selected models. The third column lists the possible range of the parameter values while the fourth column lists the selected values of the parameters by SMBioNet tool.
cr
ip t
Selected set of values 0 0 0 1 0 0 0 1 0 1 1 1 0 1 1 1 0 1
us
Range of Values 0 0,1 0,1 0,1 0 0,1 0,1 0,1 0,1 0,1 0,1 0,1 0 0,1 0,1 0,1 0 1
an
Resource {} {DEN V } {SOCS} {DEN V, SOCS} {} {DEN V } {SOCS} {T LR3} {DEN V, SOCS} {SOCS, T LR3} {DEN V, T LR3} {DEN V, SOCS, T LR3} {} {DEN V } {IF N α/β} {DEN V, IF N α/β} {} {IF N α/β}
M
Parameter KT LR3 KT LR3 KT LR3 KT LR3 KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KSOCS KSOCS KSOCS KSOCS KDEN V KDEN V
390
the pathogenic behaviour.
d
load and clear infection. The stable state shows the elevated levels of both DENV and SOCS representing
Ac ce pt e
As these cycles are time-abstract cycles therefore, their exact nature (sustained oscillation or damped oscillation) cannot be predicted in qualitative models. Thus, with time-delays, their nature can be predicted in the analysis of the hybrid model, where a convergence of the invariance kernel algorithm shows sustained oscillations, while non-convergence shows damped oscillations.
395
3.3. HyTech Results
In this section, we present the HyTech results regarding the invariance kernel of the BRN. Invariance kernels are linear time constraints represented as conjuncted relations between different clocks of the automaton. These delays provide us the timing (of activation and inhibition) for the involved proteins which govern a particular set of trajectories resulting in a cycle. If a conjuncted constraint relation of an invariance 400
kernel becomes false, the trajectories will deviate from cycles and might even result in a stable (pathogenic) state. Table 3 shows the complete trajectory of the cycle discussed in section 3.1, along with its invariance kernel. As the other cycles in the state graph (see Figure 6) themselves are small perturbations of this 17
Page 17 of 41
Table 3: Invariance Kernel of the significant cycle. The invariance kernel dictates the delay constraints that are satisfied in this cycle.
Invariance Kernel
Conjunction of constraints I-VIII:
ip t
(0, 0, 0, 0) → (0, 0, 0, 1) → (1, 0, 0, 1) → (1, 1, 0, 1) → (1, 1, 0, 0) → (0, 1, 0, 0) → (0, 1, 1, 0) → (0, 0, 1, 0) → (0, 0, 0, 0)
− |d− DEN V | ≤ |dT LR3 |
II.
+ − + − − d+ IF N + dDEN V + |dT LR3 | ≤ dSOCS + |dSOCS | + |dDEN V |
III.
− − + − d+ T LR3 + |dIF N | + |dDEN V | ≤ dSOCS + |dSOCS |
IV.
+ − − d+ SOCS ≤ dT LR3 + |dIF N | + |dDEN V |
V.
− + + − d+ SOCS + |dSOCS | ≤ dIF N + dDEN V + |dT LR3 |
VI.
− + d+ T LR3 + |dT LR3 | ≤ dSOCS
VII.
− + − + d+ SOCS + |dSOCS | ≤ dT LR3 + |dT LR3 | + dDEN V
VIII.
− + − d+ DEN V + |dT LR3 | ≤ dSOCS + |dSOCS |
us
cr
I.
M
an
Qualitative Cycle
one, the invariance kernel was computed only for this cycle, detailing the delay constraints governing the trajectory. The table also allows us to list the relations (≤, , =) between the individual delays in
d
405
Ac ce pt e
the form of a matrix (given in Table 4) computed by using a linear constraint solver over finite ranges for each delay (the computed sets are provided in Supplementary File 2).
4. Discussion
Dengue infection is the form of critical disease that shows the symptoms even before the activation 410
of adaptive immune response [76]. This infection also presents viral particles to various organs, such as liver, spleen, lymph nodes, brain and spinal cord by inhibiting the IFNα/β mediated signalling pathways [78]. Previous studies showed that DENV uses its NS proteins to down regulated the IFN induction, which is very necessary to control the DENV infection in early phases of its progression [37]. More than 300 IFN stimulated genes are induced as a result of the activation of Janus Kinase (JAK), Tyrosine Kinase
415
(Tyk) signal transducers and activation of transcription (STAT)1/2 and interferon regulatory factor (IRF) when IFNs bind to cell surface receptors [14]. For the clearance of DENV particularly in early stages of infection, TLR3 mediated IFN production is very critical in determining course of pathogenesis. However, uncontrolled TLR signalling results in chronic inflammation, infectious disease propagation and autoimmune 18
Page 18 of 41
Table 4: Relation Matrix of the significant cycle which depicts binary relations between the states. Each entry of the matrix represents whether a delay ’a’ would be greater than, equal to, or less than a delay ’b’. The matrix is read as drow ∼ dcolumn where ∼ ∈ {≤, , ≥, =}.
Relation Matrix |d− T LR3 |
d+ IF N
|d− IF N |
d+ SOCS
|d− SOCS |
d+ DEN V
d+ T LR3
=
≥, ≤
≥, ≤
≥, ≤
≥
≥, ≤
=
|d− SOCS |
≥, ≤
≥, ≤
≥, ≤
≥, ≤
≥, ≤
d+ DEN V
≥, ≤
≥, ≤
≥, ≤
≥, ≤
≥, ≤
|d− DEN V |
≥, ≤
≤
≥, ≤
≥, ≤
cr
us
≥, ≤
>
=
≤
≥, ≤
≥
=
≥, ≤
≥, ≤
≥, ≤
=
an
≥, ≤
M
S1476-9271(14)00126-1 http://dx.doi.org/doi:10.1016/j.compbiolchem.2014.10.003 CBAC 6373
To appear in:
Computational Biology and Chemistry
Received date: Revised date: Accepted date:
27-5-2014 1-9-2014 6-10-2014
Please cite this article as: Babar Aslam, Jamil Ahmad, Amjad Ali, Rehan Zafar Paracha, Samar Hayat Khan Tareen, Umar Niazi, Tariq Saeed, On the Modelling and Analysis of the Regulatory Network ofDengue Virus Pathogenesis and Clearance, Computational Biology and Chemistry (2014), http://dx.doi.org/10.1016/j.compbiolchem.2014.10.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
On the Modelling and Analysis of the Regulatory Network of Dengue Virus Pathogenesis and Clearance Babar Aslama , Jamil Ahmadb,∗ , Amjad Alia , Rehan Zafar Parachaa , Samar Hayat Khan Tareenb , Umar Niazic , Tariq Saeedb a Atta-Ur-Rahman
cr
ip t
School of Applied Biosciences (ASAB), National University of Sciences and Technology (NUST), Islamabad, 44000, Pakistan. b Research Center for Modeling and Simulation (RCMS), National University of Sciences and Technology (NUST), Islamabad, 44000, Pakistan. c IBERS, Aberystwyth University, Edward Llwyd Building, Penglais Campus, Aberystwyth, Ceredigion Wales UK, SY23 3FG.
us
Abstract
Dengue virus can ignite both protective and pathogenic responses in human. The pathogenesis is related
an
with modified functioning of our immune system during infection. Pattern recognition receptors like Toll like receptor 3 is vital for the induction of innate immunity in case of Dengue infection. Toll like receptor 3 induces
M
TRIF mediated activation of Type 1 interferons and Fc receptor mediated induction of cytokines. Interferons have been related with clearance of Dengue virus but it has adopted modified regulatory mechanisms to counter this effect. SOCS protein is also induced due to the interferon and cytokine mediated signalling which
d
can subsequently play its part in the regulation of interferon and cytokine production. Our hypothesis in this
Ac ce pt e
study relates the pathogenesis of Dengue virus with the SOCS mediated inhibition of our innate immunity. We used the qualitative formalism of Ren´e Thomas to model the biological regulatory network of Toll like receptor 3 mediated signalling pathway in an association with pathogenesis of dengue. Logical parameters for the qualitative modelling were inferred using a model checking approach implemented in SMBioNet. A linear hybrid model, parametric linear hybrid automaton, was constructed to incorporate the activation and inhibition time delays in the qualitative model. The qualitative model captured all the possible expression dynamics of the proteins in the form of paths, some of which were observed as abstract cycles (representing homoeostasis) and diverging paths towards stable states. The analysis of the qualitative model highlighted the importance of SOCS protein in elevating propagation of dengue virus through inhibition of type 1 interferons. Detailed qualitative analysis of regulatory network endorses our hypothesis that elevated levels of cytokine subsequently induce SOCS expression which in turn results into the continuous down-regulation of Toll like receptor 3 and interferon. This may result into the Dengue pathogenesis during the stage of immunosuppression. Further analysis with HyTech (HYbrid TECHnology) tool provided us with the realtime constraints (delay constraints) of the proteins involved in the cyclic paths of the regulatory network 1
Page 1 of 41
backing the evidence provided by the qualitative analysis. The HyTech results also suggest that the role of SOCS is vital in homoeostasis Key words: Formal Modelling, TLR3, Dengue, Type 1 Interferons, SOCS-1, Ren´e Thomas, Kinetic Logic,
ip t
parametric linear hybrid automata, SMBioNet.
cr
1. Introduction
The incidences of dengue (DENV) fever have increased 30-fold in the last five decades [25]. Today
us
nearly 2.5 billion people, i.e. about 40 percent of the world population, is living in areas where there is a risk of DENV transmission. It has been reported as endemic in at least 100 countries [36]. According to a recent report of World Health Organization (WHO), around 50-100 million DENV infections occur per
an
5
annum including 500,000 cases of DENV haemorrhagic fever (DHF) and 22,000 cases of mortality worldwide [28, 62].
M
Dengue virus belongs to the genus Flavivirus, family Flaviviridae. It has four serotypes including DENV1, DENV-2, DENV-3 and DENV-4 [10, 61]. Transmission of the virus is mainly caused through the bites of 10
the female Aedes aegypti mosquito. Female Aedes aegypti acquires this virus from the blood of an infected
d
person [69, 54]. The virus infects the mid-gut of the mosquito and transfers to the salivary glands over
Ac ce pt e
a period of 8-12 days, which then is transmitted to humans on next bite [74]. The DENV virus after entering the systemic circulation replicates within the cell of mononuclear phagocyte lineage (macrophages, monocytes and B-cells). DENV infections may result in dengue fever (DF), dengue haemorrhagic fever 15
(DHF), or dengue shock syndrome (DSS) marked with high fever, skin rash, headache and severe muscular and joint pains [56, 74].
1.1. TLR3 Signalling Pathway
Recent studies implicates that the course of DENV infection has been attributed to the innate immunity, particularly generated by the pattern-recognition receptors (PRRs) such as Toll-like receptors (TLRs) [24, 20
96, 82, 93]. TLRs are the part of the innate immune system and activated by the recognition of pathogen associated molecular patterns (PAMPs) [7, 47]. The study of Iwasaki and Medzhitov on human monocytic cells U937 strongly supports the role of TLR3 for inducing immune response against DENV [46]. Figure 1 ∗ Corresponding
author Email addresses: [email protected] (Jamil Ahmad)
Preprint submitted to Computational Biology and Chemistry
October 13, 2014
Page 2 of 41
illustrates the possible signalling events generated by TLR3 during DENV infection. As key players of innate immunity, the recognition of several viruses via TLR is very important for the induction of various 25
cytokines to avert pathogenesis [46]. Several pro-inflammatory cytokines are induced as a result of the TLR3 activation and provide an interface between the innate and adaptive immune responses [43]. Modulation
ip t
of DENV infection by type 1 interferon alpha/beta (IFNα/β) has been suggested by observing DENVinfected patients and IFNα/β-receptor deficient mice [9]. In vitro analysis with several cell types and viral
30
cr
strains proved that DENV viral RNA decreased by multiple folds upon introduction of IFNα/β, however, IFN gamma (γ) showed variable responses depending on the cell type and signalling pathway [27]. At the
us
early stage of viral infection, inflammatory response through PRRs provide dominant role in priming and regulating innate immune system [55]. Yi-Ting Sai et al. in their study demonstrated the role of TLR3 in mediating strong IFNα/β induction and subsequent inhibition of DENV viral replication, which limits its
35
an
cytopathic effect [93, 27]. As Type 1 interferons upon their interaction with IFN receptors elicit antiviral, anti-proliferative and immunomodulatory roles. Shresta et al. demonstrated critical and non-overlapping
M
role of IFNα/β in resolving primary DENV infection in genetically deficient mice [78].
Ac ce pt e
d
Figure 1: Pathway of TLR3. DENV is sensed by TLR3 in endosomes [8]. TLR3 then activates TRIF which initiates the up regulation process of NFκB and IFNα/β by two distinct pathways [44]. Type 1 IFN pathway TRIF via NAP 1 and TRAF3 interacts with TBK1 which transduces a signal to Ikki which leads to phosphorylation of IRF [33]. Phosphotyrosine residues on TLR3 binds to P13K which is required for complete IRF activation [77]. Dimerized IRFs then translocated into nucleus where they bind to DNA and in conjugation with the co-activator CBP/p300 they transactivate type 1 IFN gene expression [63]. Pro inflammatory cytokine pathway. NFκB up regulation by TLR3 is mediated by TRAF6 interaction with TRIF that involves the TAK1/TAB2/TAB3 complex which results in activation of TAK1 [83]. TAK1 then phosphorylate and activate IKK. Then NFκB nuclear localization sequence is unmasked by degradation effect of phosphorylated IKB. Activation of NFκB and transcriptional induction of target genes requires signalling from TRIF via RIP1 [23]. NFκB up regulation via TRIF, TRAF 6 results in up regulation of IL-10 which leads to expression of SOCS via JAK / STAT pathway [65, 50]. SOCS can inhibit this signalling pathway at four sites including degradation of NFκB [57] , TLR3 [92] , IFNα/β and JAK / STAT. IFNα/β have been seen in regulation of DENV whereas NS4B protein of DENV inhibit the role of IFNα/β [60]. IFNα/β can also lead to induction of SOCS [67]. DENV by its ADE (Antibody dependent enhancement) mechanism through FcR signalling [16, 73, 19] up regulates IL-10 production [32, 17].
Expression of DENV non-structural (NS) proteins including NS2A, NS4A or NS4B in human cells support the elevated replication of an IFN sensitive virus [60]. Cells infested with NS4B protein or DENV do not exhibit IFN-α/β mediated induction of nuclear signal transducer and activator of transcription (STAT)1. 40
By blocking STAT1 activities, NS4B protein of DENV may interfere with the IFNα/β mediated antiviral response [60]. These findings indicate that NS4B might be involved in blocking IFNα/β signalling during DENV infections. IFNα/β activate STAT1–STAT3 heterodimers [66], followed by SOCS1 and SOCS3 mediated downregulation [80]. Qin et al. implicated in their study about the up-regulation of IFN mediated induction 3
Page 3 of 41
45
of SOCS-1, which in return down-regulate IFN signalling and reduce the expression of CD40 cells [67, 80]. They also proposed the inhibitory mechanism of JAK/STAT phosphorylation which ultimately inhibit the vital signalling cascade for the induction of cytokines [67]. Various studies have established the role of Antibody-dependent enhancement (ADE) of dengue virus
50
ip t
pathogenesis [53]. The pre-existing antibodies (Abs) are major contributors for development of DHF and DSS and fail in neutralizing the infecting serotype of DENV [38]. The Fc receptors (FcR) family offers a major
cr
example of how to promptly trigger, activation and inhibition of the signalling pathway to maintain harmony between innate and adaptive immune response [70]. In accumulation to their exact roles in prompting the
us
activation of innate effector cells, FcγRs also play a role in the regulation of B-cell activation, antigen presentation and immune-complex mediated maturation of the dendritic cells (DCs) [13]. A new hypothesis 55
termed as intrinsic ADE, suggests that DENV internalization via FcγR results in the induction of interleukin-
an
10 (IL-10) production, biased T-helper-1 (Th)1 to Th2 response and suppression of innate immunity, which leads to both high levels of antibodies (Abs) and viral load in dengue patients [94, 39]. Fujii et al. in a recent
M
study, suggested the up-regulation of IL-10 levels via FcγR [32]. As a consequence of ADE, complementary mechanism to higher viremia may be explained as FcγR mediated DENV entry due to suppressed antiviral 60
immune response. In the case of DENV infection, suppressed levels of IL-12, IFNγ, TNFα and nitric
Ac ce pt e
environment for viral replication [17].
d
oxide along with the elevated levels of anti-inflammatory cytokines tend to provide immunosuppressive
1.2. Computational Systems Biology
Computational systems biology provides a compelling platform to answer critical scientific questions 65
about signalling mechanisms based on pragmatic modelling and theoretical elucidation [51]. Multitude of biological facts have been uncovered by molecular biology including functional properties of the genes and proteins. However, for the interpretation of the dynamics of the biological system, these informations require further investigation [51]. To understand the intrinsic complexity of biological systems experimental and computational approaches in combination are speculated to solve this problem [52, 45]. Various com-
70
putational approaches assisted with the experimentally observed phenomenon in previous literature helped us to understand the behaviours of the complex biological regulatory networks involved in DENV replication and immune evasion [45]. A model that captures biological observations related to physiological and pathological conditions could also predict the underlying dynamics behind these observations: oscillations or cycles, stable states and bifurcation points. Interestingly, the bifurcation points (specific pattern of protein
75
interactions), may play a role in leading the system towards pathogenic state, can provide us an insight 4
Page 4 of 41
about the useful therapeutic target. Thus, instead of targeting proteins in the pathogenic state (abnormal physiological state) drugs may target entities involved in the interactions causing the physiological abnormalities. We predicted with the help of a qualitative model the underlying causes of the hyper-inflammatory condition in the case of cerebral malaria (CM) pathogenesis [4]. The model predicted that longer activation and interactions of two proteins BTK and MAL, and the pre-existing high levels of cytokines play key role
ip t
80
in the pathogenesis of CM. The later hypothesis was validated with a study conducted by Franklin et al.,
cr
2011 [31, 49, 40, 30, 68] which emphasizes the relationship between elevated levels of inflammatory cytokines and CM. Our recent study [64] shows a strong affinity of BTK with MAL protein. These studies suggest
85
us
that the expression dynamics of the pathway should be studied in order to have an understanding of the behaviours of the systems and to predict drug targets which would ultimately increase our confidence about the potential drugs and would accelerate the drug design process.
an
For the last few decades, formal methods have received the attention of researchers in various disciplines to formally model and verify complex systems. Model checking is an automatic formal verification technique
90
M
which exhaustively explore dynamics of a system formally expressed as automata (M) to check whether it meets specifications formally expressed as formulas (Φ) [21]. Formally, a model checking process is represented as M |= Φ. Verification could be carried out even if the system is parametric and in this case
d
the values of the parameters are synthesised for the properties to be verified [81, 42]. Thus, parameter
Ac ce pt e
synthesis is another key feature of model checking.
In the current study, TLR3 mediated signalling pathway has been studied using the generalized logical 95
approach introduced by Ren´e Thomas [87, 90]. This approach represents the qualitative analysis of dynamic behaviour of the system with a focus on the feedback mechanisms. This approach is successfully applied on some other signalling systems including the MyD88-adapter-like (MAL) associated biological regulatory network (BRN) [4], control of differentiation process in helper T cells [58], control of organ differentiation in Arabidopsis thaliana flowers [59] and segmentation during embryo-
100
genesis in Drosophila melanogaster [84, 75]. Kinetic formalism of Ren´e Thomas and its implementation in GENOTECH [4] provided us an in silico approach to test hypotheses related to DENV pathogenesis. Ren´e Thomas’ logic relies on the discretization of expression levels which can be determined by the qualitative thresholds. The evolution of these levels are driven by the logical parameters. Values of these parameters can be inferred with the help of SMBioNet [48], which implies a model checking approach by first expressing
105
the observations (properties) in Computation Tree Logic (CTL) formulas and then computes all possible sets of parameters leading to models satisfying these observations. Once the reliable set of parameters 5
Page 5 of 41
are determined, the tool GENOTECH then computes the qualitative model showing all the behaviours of the network and where important behaviours can be seen as paths (trajectories), closed paths (cycles) and stable states (see state graph in Section 3). Delays are incorporated in the qualitative model to model the 110
activation and inhibition orders of proteins; which protein activates or deactivates first now depends on their
ip t
time delays. The values of such parameters can be synthesised using the hybrid model checking tool HyTech [41].
cr
1.3. Our Contribution
In this study, we performed a formal analysis of TLR3 signalling pathway using its discrete and hybrid models. In these models, we analysed the role of SOCS protein in the context of TLR3 mediated IFNα/β
us
115
production in DENV infection. The design of this study is graphically depicted in Section 2. First, the
an
Biological Regulatory Network (BRN) is obtained (Figure 2) by keeping the prominent entities (in which we are interested) of the pathway and abstracting the rest, provided that the behaviour is preserved. For +
−
→) ‘B’ which inhibits(− →) ‘C’, and the proteins ‘A’ and ‘C’ are the example, if an entity ‘A’ activates(− −
→ C while implicitly showing the behaviour ones in which we are interested, we simply abstract it to A −
M
120
that ‘A’ inhibits ‘C’ through the activation of ‘B’. This BRN is then used to obtain the logical parame-
d
ters from SMBioNet. These parameters are used in GENOTECH which computed the qualitative cycles,
Ac ce pt e
a stable state, and transformed it into parametric linear hybrid automata in HyTech format. Cycles were showing the oscillations of all proteins while the stable state was showing the over-expression of SOCS along 125
with DENV. The delay constraints of a significant cycle, representing its invariance kernel (a set of cyclic trajectories which start and end at the same point, see Definition 12, [5]), were computed by the analysis of parametric linear hybrid automaton (hybrid model) of TLR3 signalling pathway in HyTech tool. The computation of the invariance kernel enabled us to observe the sustained cyclic nature of this cycle, as well as the possible delay values in the form of linear constraints. Our study implicates that when SOCS was
130
oscillating within the qualitative cycles then all the other proteins of the pathway including SOCS itself were showing sustained oscillations in the hybrid model. Our analysis suggests that altered oscillatory behaviours of SOCS protein play a vital role in the inhibition of TLR3 mediated IFNα/β production, which ultimately results in the down-regulated immune response with subsequent elevation in DENV titre.
6
Page 6 of 41
Figure 3: Work flow diagram of the design of the study.
2. Methods
cr
135
ip t
Figure 2: TLR3 associated BRN. Reduced form of the TLR3 induced BRN shown in Figure 1. Entities are represented by vertices whereas edges represent interaction between them. An association between two corresponding entities is represented by each edge. Interaction types are represented by sign on arc i.e. negative for inhibition (dotted arrows) and positive for activation (solid arrows). Threshold levels of entities which are necessary to activate or inhibit its successor entity (see Definition 2 in Method section and Table 2) are represented by given values on the arc.
us
Usually, biological systems are represented by conventional models having differential equations referring to the time derivative of different quantities such as concentrations, temperatures and rates. As a result of inherent complexity of biological systems, little is known exactly about these entities. In the year 1970, Ren´e
140
an
Thomas [87] presented a Boolean logic based method to model biological regulatory networks. He presented the efficacy of his modelling methodology by checking the lambda phage genetic switch [87, 59, 58]. Ren´e
M
Thomas noted the limitations in Boolean logic as it has only two levels, 0 or 1, that are not able to formalize related problems. So, Thomas generalized the Boolean logic to kinetic logic and proved the efficacy of the methodology by utilizing it in different genetic regulatory systems [85, 88, 87, 79, 90, 91, 26, 89], and found
145
Ac ce pt e
d
it to be efficient in approximating the differential equations based modelling [88]. 2.1. The Kinetic Logic formalism of Ren´e Thomas Kinetic logic formalism models biological regulatory networks (BRNs) based on the qualitative threshold levels of entities and their types of interactions [88]. Value of at least 1 as minimum threshold for any entity is required to bring a positive or negative change in the target entity [86]. Entities are representing proteins, genes or both of them that can interact with each other to bring positive or negative effects showing 150
stimulatory or inhibitory effects on the system, respectively. Positive effect shows that the expression level of an entity increases the rate of activation or production of other entity, whereas, negative effect reduces the rate of activation. Biological activation and inactivation are sigmoidal (non-linear bounded curve) in nature and are discussed below. The following definitions and concepts define the kinetic logic formalism and are mainly adopted from [4].
155
2.1.1. Activation If an entity e1 enhances the rate of synthesis or activation of another entity e2 then e1 is the positive regulator for e2 . In this case, every time the concentration of e1 enhances or reduces, it will affect the 7
Page 7 of 41
concentration of e2 in similar manner. This relationship can be shown by sigmoidal nature. The sigmoidal curve can be shown by its step function and is shown in Figure 4 A. We can see that below a threshold 160
level θ, e1 slightly changes the rate of activation of e2 . But as soon as the concentration of e1 reaches its threshold level θ, the rate of activation of e2 elevates rapidly to a level where it saturates its target entity
ip t
e2 . So, we can say that when e1 < θ it can not bring any positive change and if e1 ≥ θ then e1 can bring
cr
the positive change in e2 .
165
an
us
Figure 4: Dummy tendency graph of biological regulation. (A)Positive regulation of a target entity e2 by an activator e1 is presented in upper left corner of sigmoidal graph , whereas upper right one is its step function.(B) Negative regulation of a target entity e4 by an inhibitor e3 is presented in lower left corner sigmoidal graph, whereas right one is its step function. There is no or little effect on its target entity when the regulator remains below threshold level (represented by θ on X-axis). The rate of activation or inactivation of target entity leads it to its saturated level when regulator reaches to its threshold concentration level.
2.1.2. Inhibition
M
If an entity e3 decreases the rate of activation of another entity e4 , then e3 can be called as negative regulator or inhibitor of e4 and the process is called as negative regulation or inhibition. The process of
d
inhibition also depends on the threshold levels of e3 showing that whenever the concentration level of e3
170
Ac ce pt e
reaches its threshold level then e3 will be efficient in decreasing the activation levels of e4 , up to a level where e4 can be completely inactivated or its affect in the system is abolished. The affect of inhibition is also of sigmoidal nature and is shown in Figure 4 B along with its step function. 2.2. Semantics of the Kinetic Logic Formalism
This part shows the semantics of Thomas formalism with a toy example of a BRN comprising of three entities (e1 , e2 and e3 ) shown in Figure 5. 175
Figure 5: Dummy BRN. An example of BRN in which interaction between genes e1 , e2 and e3 are labelled with + sign representing activation and integer ”1” representing the threshold levels.
Definition 1 (Directed Graph). A directed graph G(N, E) is a tuple, where: set of all nodes is N and set of ordered pair of edges or arcs is denoted by E ⊆ N × N 8
Page 8 of 41
(e1 , e2 ) is an edge which is considered to be directed from e1 to e2 , where e1 is called the head and e2 is 180
called the tail. In G(N, E), G− (e) and G+ (e) denote sets of predecessors and successors, respectively, of a node e ∈ N . N = {e1 , e2 , e3 } and E = {(e1 , e2 ), (e2 , e3 ), (e3 , e1 )} are sets of nodes and edges, respectively,
and G+ (e1 ) = {e2 } represents the predecessors and successors of e1 , respectively.
Definition 2 (Biological Regulatory Network). A biological regulatory network (BRN) is a labelled directed graph G(N, E), where the set of nodes N represents biological entities and interactions between entities are represented by the set of edges E ⊆ N × N , in which each edge ei → ej is labelled by a couple (jei ej , ηei ej ), jei ej being a positive integer representing qualitative threshold and ηei ej being either (+) or (-) sign representing activation or inactivation, respectively. Moreover, we define:
190
us
cr
185
ip t
which are shown as network of genes in Figure 5 . In case of gene e1 , (set of regulators of e1 ) G− (e1 ) = {e3 }
The outdegree of ei is the total number of targets of ei , it is denoted by pei and it is such that ∀ek ∈ G+ (ei ) jei ,ek pei .
an
Set Zei = {0, 1, ...., rei } of an entity ei represents its qualitative levels (concentration levels) where rei = max{jei ,ek |ek ∈ G + (ei )} To capture the behaviours of a BRN, it is important to define the set of all possible number of qualitative states.
M
195
d
Definition 3 (Qualitative States). A qualitative state s ∈ S of a BRN G(N, E) is a tuple which denotes a configuration of the qualitative expression levels of all entities. The set of states is :
200
Ac ce pt e
S = Πei ∈N Zei
Vector (sxei )∀ei ∈N represents the qualitative states, where level of expression (concentration) of a product ei is denoted by xei .
One associates with each qualitative state, the set of resources of each entity which corresponds to the set of predecessors which effectively control the entity (either a present activator or an absent inhibitor).
205
Definition 4 (Resources). At level xei , the set of resources Qxei of an entity ei ∈ N is defined as Qxei = {ej ∈ G− (ei )|(xej ≥ jej ei and ηej ei = +) or (xej < jej ei and ηej ei = −)}. In the above definition, it is to be noted that the inhibitor is considered as an activator in its absence.
210
Definition 5 (Logical Parameters). Logical parameters influence the dynamic behaviour of a BRN. The set of logical parameters are defined as K(G) = {Kei (Qxei ) ∈ Zei ∀ei ∈ N }. Logical parameters and set of resources are shown in Table 1 for the running example of three entities (Figure 5). To define the transition between two qualitative states, it is convenient to define the evolution operator [11]: 9
Page 9 of 41
Definition 6 (Evolution Operator). The evolution operator is defined as, ⎧ ⎨ xei + 1 if xei < Kei (Qxei ); xe − 1 if xei > Kei (Qxei ); xei Kei (Qxei ) = ⎩ i if xei = Kei (Qxei ). xei
cr
where xei and Kei (Qxei ) ∈ {0} ∪ N
ip t
Figure 6: State graph of the dummy BRN. State graph of BRN presented in Figure 5. Qualitative (discrete) states are presented as nodes. Qualitative levels of entities e1 , e2 and e3 , respectively, are shown as values inside a node. The state (0,0,0), and (1,1,1) represent deadlocks, whereas as transition from the state (1,0,1) to (1,0,0) creates a cyclic trajectory.
us
Table 1: Table of states, resources and logical parameters
xe2
xe3
Qxe1
Qxe2
Qxe3
Ke1 (Qxe1 )
Ke2 (xxe2 )
Ke3 (Qxe3 )
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
0 1 0 1 0 1 0 1
{} {e3 } {} {e3 } {} {e3 } {} {e3 }
{} {} {} {} {e1 } {e1 } {e1 } {e1 }
{} {} {e2 } {e2 } {} {} {e2 } {e2 }
0 1 0 1 0 1 0 1
0 0 0 0 1 1 1 1
0 0 1 1 0 0 1 1
d
M
an
xe1
215
220
Ac ce pt e
All the states (xe ), set of resources (Qxe ) and logical parameters Ke (Qxe ) of the dummy BRN (Figure 5)are listed.
Definition 7 (State Graph). Let G(N, E) be a BRN, let s = (sxei )i ∈ S be a state and let consider the family of logical parameters K(G). The state graph of the model is the directed graph noted G = (S, T ) where S is the set of states and T ⊆ S × S is the relation between states called transition relation such that s → s ∈ T if and only if: 1. ∃ a unique ei N such that sxei = sxe and sxe = sxei Kei (Qxei ), and i i 2. ∀ ej N \ {ei } sxe = sxej . j
Table 1 represents the state table of running example which shows all the possible states, set of resources and logical parameters. This table helps to construct the state graph which captures all the behaviours of the running example. Two stable states namely (000) and (111) can be observed clearly in the state state graph shown in Figure 6. 225
2.3. Discrete Modelling of DENV associated BRN using GENOTECH GENOTECH [4] constructs a qualitative model of BRNs according to Thomas’ formalism and produce the state graph, which shows the dynamical behaviour of a BRN. Modelling is facilitated by GUIs, which 10
Page 10 of 41
allow to construct a BRN as directed labelled graph. The tool can perform the basic analyses of the state graph like the identification of stable states, cycles and path finding between two states. The state graph 230
can be saved in DOT format [34] for further study in the Graphviz tool [29].
ip t
2.4. Model Checking Although Ren´e Thomas framework captures the qualitative dynamics of biological regulatory networks, its semantics are derived from the values of logical parameters. Once the state graph is generated by assuming
235
cr
a certain value of these parameters, it is still required to manually verify the model against some known observations. This manual process is quite time expensive even for small BRNs containing only two entities.
us
Model Checking is one of the formal techniques in computer science that verifies models against specification formally expressed as temporal logic formulas. Bernot et al.[12] applied model checking to reverse engineer
an
the values of logical parameters from biological observations specified as CTL (Computational Tree Logic) formulas. This technique is based on exhaustive enumeration of parameters values, each time a parameter 240
combination is used to generate the state graph. Finally, by invoking a model checker, the verification
M
of CTL formulas is carried out and all parameters that are coherent with CTL formulas are selected as set of acceptable parameters. Thus model checking solves this problem by not only expressing biological
245
Ac ce pt e
2.4.1. Quantifiers
d
behaviours in a formal temporal language but also reverse engineering the values of logical parameters.
A path is a sequence of states in a state graph. In Computational Tree Logic (CTL), properties are checked in the state graph by applying quantifiers to the paths originating from a given state, as well as the states occurring on each of the originating paths. Description of each of the quantifiers is given below: A : This is a path quantifier which enforces that a given property should hold in all paths originating from the given state. The quantifier itself is read as “For all paths”.
250
E : Known as “Existential Quantifier”, this is also a path quantifier which enforces that a given property must hold in at least one path originating from the given state. The quantifier is read as “There exists a path”. G : This quantifier is known as the “Global Quantifier” and is a state quantifier which enforces that a property holds in all states of a path originating from the given state, inclusive of the given state as
255
well. It is read as “globally”. 11
Page 11 of 41
F : The “Future Quantifier” is the second state quantifier and enforces that a given property must hold in one of the future states in the path originating from the given state. It is read as “in future” or “eventually”. The future quantifiers also covers the current/given state as well when checking the property. X : The “Next Quantifier” is the third state quantifier and enforces that a given property must hold
ip t
260
in the immediate successor state. It is read as “next”.
cr
The usual format for a CTL formula is φ = {A/E}{G/F/X}{ψ/φ} where ψ is a predicate representing
us
a property and φ another CTL formula. Then, this format allows nesting of quantifiers. 2.4.2. SMBioNet
SMBioNet [12, 71, 48, 72] is a tool for the parameter enumeration of Biological Regulatory Networks
an
265
(BRNs), based on the Discrete Framework of Ren´e Thomas [87]. Given a model of BRN in the form of Thomas network and behavioural properties (observations), expressed as CTL formulas, SMBioNet exhaus-
M
tively enumerates all compatible parametrizations by generating a state graph for each parameter combination and by verifying the formulas on each state graph. The verification of CTL property is performed by invoking model checker NuSMV [20]. The parameter combinations are reduced by applying Snoussi and
d
270
observability constraints [15]. Finally, all the models that satisfy CTL properties are short-listed. SM-
Ac ce pt e
BioNet has been applied in studies like tail resorption in tadpole metamorphosis [48] and immunity control in bacteriophage lambda [71].
2.5. Hybrid Modelling to Enhance Qualitative Modelling with Delays 275
Hybrid modelling merges the discrete modelling with the continuous domain of time via special variables called ’clocks’. Likewise, the transitions are characterised as discrete and continuous where, Discrete transitions represent instantaneous switching of locations or states, and Continuous transitions represent elapsing of time within the same location or state. The natural environment of the cell is difficult to accurately model through discrete modelling frame-
280
works. Take the concentration of a protein for instance – the discrete model represent their expression levels as positive integers where a protein is at ‘0’ level at one instance, and immediately jumps to ‘1’ level the very next instance. This does not accurately replicate the continuous process of increasing expressions that proteins go through in the natural cellular environment. Ahmad et al. [6] proposed the use of piecewise 12
Page 12 of 41
linear curves for the modelling of the sigmoidal nature of protein expression (Figure 7), as opposed to the 285
step functions used in discrete models.
ip t
Figure 7: Activation and Inhibition delays. The clock ha 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+/− .
Clock variables are used to measure the ‘delays’ (the time duration) that needs to pass between two
cr
consecutive expression levels. Thus a clock variable h is associated with each protein in the BRN. The initial values of each h are set to zero, which then approach either d+ or d− . d+ signifies production delay,
290
us
that is, the delay required to increase the associated proteins concentration level by 1. Similarly, d− signifies the degradation delay, that is, the delay to decrease the said concentration by a single level. The rate of evolution of each h is given by the first order derivative dh/dt = r where r ∈ {0, 1, −1} [6, 2].
an
In most cases, the exact values of the delays associated with the proteins are not known which is why unvalued parametric delays are instead used [4]. Thus, the hybrid model was constructed using Parametric
295
M
Bio Linear Hybrid Automaton [1] defined below and given as Supplementary File 3. Let C = (X, P ), C ≤ (X, P ), and C ≥ (X, P ) be the set of constraints using only =,≤, and ≥, respectively. Here, X and P are the sets of real valued variables and parameters, respectively.
Ac ce pt e
d
Definition 8 (Parametric Bio Linear Hybrid Automaton (Bio-LHA). A parametric Bio linear hybrid automaton B is a tuple (L, l0 , X, P, E, Inv, Dif ) where, L is a finite set of locations, 300
l0 ∈ L is the initial location,
P is a finite set of parameters (delays),
X is a finite set of real-valued variable (clocks), E ⊆ L×C = (X, P )×2X ×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 and the reset set R ⊆ X, 305
Inv : L → C ≤ (X, P ) ∪ C ≥ (X, P ) 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 T, where T∗ = T \ {0}.
310
Figure 8 represents the conversion of the discrete state graph of Figure 6 into a Bio-LHA. Definition 9 (Semantics of Bio-LHA). Let γ be a valuation for the parameters P and ν represents the values of clocks in a location. The (T, γ)-semantics of a parametric Bio-LHA B = (L, 0 , X, P, E, Inv , Dif ) is defined as a timed transition system B = (S, s0 , T, →) where: (1) S = {( , ν) | ∈ L and ν |= Inv ( )}; (2) s0 is the initial state and (3) the relation →⊆ S × T × S is defined for t ∈ T as: 13
Page 13 of 41
315
ip t
Figure 8: The Parametric Bio Linear Hybrid automaton of the running example of the dummy BRN based on its state graph represented in Figure 2. Each location is marked with its configuration (state) name (for example (1,0,0)), the invariants (conjuncted constraints), and the entity rates. Each transition is marked with a guard and clock reset of the respective changing entity. Location (1,0,0) is the initial location, represented as a small arrow on its upper-left side. Locations (0,0,0) and (1,1,1) are representing deadlock states as: (i) there are no outgoing transitions from these locations and (ii) all rates are equal to zero. The rates of the entities in each location are based on the evolution in the first and second successor locations. 0
discrete transitions: ( , ν) −−→ ( , ν ) iff ∃( , g, R, ) ∈ E such that g(ν) = true, ν (h) = 0 if h ∈ R and ν (h) = ν(h) if h ∈ / R.
cr
t
us
continuous transitions: For t ∈ T∗ , ( , ν) −−→ ( , ν ) iff = , ν (h) = ν(h) + Dif ( , h) × t, and for every t ∈ [0, t], (ν(h) + Dif ( , h) × t ) |= Inv ( ), where |= represents satisfaction operator. Using the semantics of the Bio-LHA, Ahmad et al. [2] were then able to define the temporal state space and the invariance kernel set which have been adapted below.
an
320
325
M
Definition 10 (Temporal Zone). Temporal zone is defined as a region where time elapses until a discrete transition between states takes place. Definition 11 (Temporal State Space). The temporal state of a BRN is composed of the complete set of temporal zones derived from the discrete model of the said BRN.
d
In the hybrid model of the BRN, the sequence of points in a trajectory is denoted by φ(t) for t ∈ R ≥ 0 [3].
Ac ce pt e
A particular trajectory is said to be viable according to viability domain S, iff it remains within a prescribed region S [3]. 330
Definition 12 (Invariance kernel). A trajectory φ(t) is said to be viable in S if φ(t) ∈ S ∀ t ≥ 0. A subset K of S is an invariant if for any point p ∈ K, a trajectory starting in p is viable in K. An invariance kernel K is the largest invariant subset of S.
3. Results
3.1. Logical Parameters
The BRN given in Figure 2 comprises four entities; namely DENV, IFNα/β, TLR3 and SOCS. DENV 335
activates TLR3 and SOCS, which then formulates a positive feedback loop with IFNα/β due to the even number of negative interactions involved in this loop. The presence of SOCS inhibits the synthesis of IFNα/β and TLR3, while its own production is positively regulated by IFNα/β. Thresholds fixed by variables TLR3, IFNα/β, SOCS and DENV lead to a total of 16 possible states in the state graph. The state of the system at any given time is represented by a vector (TLR3, IFNα/β, SOCS, DENV). Before the infestations of DENV,
340
all the proteins that participate in this network are absent. Hence, (TLR3=0 & IFNα/β=0 & SOCS=0 14
Page 14 of 41
& DENV=0) constitutes a starting state of the system. The presence of two responses within the same network, derived from the literature, implies that there exist two possible paths starting from the initial state i.e. from all zeros state, can proceed towards DENV clearance and a pathogenic state [95]. The DENV clearance or normal state is known to be characterized by low titres of DENV, whereas the pathogenic state is characterized by its high titre and elevated levels of SOCS [97]. These two observations are encoded in
ip t
345
cr
CTL formulas as described in following equations.
P athogenic = (SOCS = 1 & DEN V = 1) Φα = Init ⇒ EF (AG(pathogenic))
an
Φβ = Init ⇒ AX(EF (Init))
us
Init = (T LR3 = 0 & IF N α/β = 0 & SOCS = 0 & DEN V = 0)
(2) (3) (4) (5)
M
Φ = Φα ∧ Φβ
(1)
The desirable qualitative model of DENV should contain closed or cyclic paths, which originate from an initial state, represented by all zeros state, and then return back to its original position from where
diseased state can be achieved if homoeostasis could not be maintained. Along with the presence of closed
Ac ce pt e
350
d
it started. Closed or cyclic paths actually represent the homoeostasis of the immune system, however, a
paths or cycles in the dynamics of DENV infection, pathogenesis can be encountered if dynamics proceed towards the stable steady state (0,0,1,1). The first formula, denoted as Φα models the behaviour of the BRN that reflects that pathogenic state is reachable from an initial state and remains as its future state. If the future state persists for a long time then we can call this as stable state. 355
The second property, given as Φβ expresses that starting from an initial state, the immune system is able to reset its configurations by going back to its original position and maintain the immune homoeostasis. The formula Φ combines (conjuncts) these two properties and is used for the selection of models by SMBioNet. SMBioNet selects 43 sets of parameters (each set represents a specific qualitative model)(Supplementary file 1) that satisfy the above mentioned CTL formula. These qualitative models were then further analysed
360
in GENOTECH by generating state graphs, which allowed us to analyse the expression patterns of the proteins, cyclic paths, stable states and the dynamics of the BRN. The selected set of parameters given in Table 2 were used to model the BRN Figure 2, which gives a state graph shown in Figure 9. This model was selected since it incorporates the above mentioned observations in the form of a stable state (0,0,1,1) 15
Page 15 of 41
and homoeostasis, represented by the qualitative cycles visible in Figure 9. In particular, we are interested 365
in the cycle which follows well ordered regulation pattern evident in the TLR3 pathway (Figure 1) where the occurrence of DENV activates TLR3 which in turn activates IFN before the activation of SOCS. The activation of IFN suppresses DENV resulting in the down regulation of TLR3 and in the up regulation of
ip t
SOCS. Finally, IFN is down regulated by the presence of SOCS, which is then degraded in the absence of its activators.
This biologically plausible pattern is found in 9 models out of 43, and is represented as (0,0,0,0) →
cr
370
(0,0,0,1) → (1,0,0,1) → (1,1,0,1) → (1,1,0,0) → (0,1,0,0) → (0,1,1,0) → (0,0,1,0) → (0,0,0,0). The detailed
us
method for the selection of a single parameter set (given in Table 2) from the 9 sets is provided in Supplementary file 5. It is important to note that the change in the parameters (that is, if any parametrisation from the other enumerated sets is selected) will change the qualitative dynamics of the BRN. This in turn will change the hybrid model (see section 3.3), affecting the delay constraints themselves, as well as their
an
375
M
biological implications and significance.
Ac ce pt e
d
Figure 9: State graph of physiological signalling of TLR3. The states, represented in oval shapes are named by the qualitative expression levels of each entity in the order TLR3, IFNα/β, SOCS, and DENV. The normal stable state (0,0,0,0) highlighted in blue is taken as initial state. The recommended diseased state after pathogenesis of DENV is (0,0,1,1) which is presented in red colour. Transitions from one state to the other are shown by different types of edges. Bold edges (condition 1) represent the transitions which may eventually lead to the clearance of DENV infection. Short dashed edges (condition 2) represent the transitions which make diverge the system towards the SOCS and DENV mediated pathogenic state. Dotted edges (condition 3) reflect that some states can converge towards the cyclic behaviour of the system.
3.2. Qualitative Model
This section presents the results obtained by the modelling and analysis in GENOTECH, which simulate the model according to Ren´e Thomas’ formalism. This tool is used to model BRN along with the usage of 380
set of logical parameters obtained from SMBioNet. As output, this gives a state graph (Figure 9) showing the paths between any two states (transitions), stable states and possible cycles among different states. The state graph represents all possible transitions from one state to another, where each state shows the expression levels of every entity at a particular time. From a stable state, which is like a sink, the system does not evolve to other states so has been termed as stable (pathogenic) state.
385
The state (0,0,0,0) represents that all entities are at zero level while (0,0,0,1) shows the presence of DENV. These two states are included in all the cycles of the BRN as shown by solid arrows (Figure 9). The tendencies of divergence from cycles towards stable state (0,0,1,1) are also shown by dashed edges. The cycles in the state graph represent the homoeostasis of the BRN to overcome or control the level of the viral 16
Page 16 of 41
Table 2: Possible parameter values for selected models. The third column lists the possible range of the parameter values while the fourth column lists the selected values of the parameters by SMBioNet tool.
cr
ip t
Selected set of values 0 0 0 1 0 0 0 1 0 1 1 1 0 1 1 1 0 1
us
Range of Values 0 0,1 0,1 0,1 0 0,1 0,1 0,1 0,1 0,1 0,1 0,1 0 0,1 0,1 0,1 0 1
an
Resource {} {DEN V } {SOCS} {DEN V, SOCS} {} {DEN V } {SOCS} {T LR3} {DEN V, SOCS} {SOCS, T LR3} {DEN V, T LR3} {DEN V, SOCS, T LR3} {} {DEN V } {IF N α/β} {DEN V, IF N α/β} {} {IF N α/β}
M
Parameter KT LR3 KT LR3 KT LR3 KT LR3 KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KIF N α/β KSOCS KSOCS KSOCS KSOCS KDEN V KDEN V
390
the pathogenic behaviour.
d
load and clear infection. The stable state shows the elevated levels of both DENV and SOCS representing
Ac ce pt e
As these cycles are time-abstract cycles therefore, their exact nature (sustained oscillation or damped oscillation) cannot be predicted in qualitative models. Thus, with time-delays, their nature can be predicted in the analysis of the hybrid model, where a convergence of the invariance kernel algorithm shows sustained oscillations, while non-convergence shows damped oscillations.
395
3.3. HyTech Results
In this section, we present the HyTech results regarding the invariance kernel of the BRN. Invariance kernels are linear time constraints represented as conjuncted relations between different clocks of the automaton. These delays provide us the timing (of activation and inhibition) for the involved proteins which govern a particular set of trajectories resulting in a cycle. If a conjuncted constraint relation of an invariance 400
kernel becomes false, the trajectories will deviate from cycles and might even result in a stable (pathogenic) state. Table 3 shows the complete trajectory of the cycle discussed in section 3.1, along with its invariance kernel. As the other cycles in the state graph (see Figure 6) themselves are small perturbations of this 17
Page 17 of 41
Table 3: Invariance Kernel of the significant cycle. The invariance kernel dictates the delay constraints that are satisfied in this cycle.
Invariance Kernel
Conjunction of constraints I-VIII:
ip t
(0, 0, 0, 0) → (0, 0, 0, 1) → (1, 0, 0, 1) → (1, 1, 0, 1) → (1, 1, 0, 0) → (0, 1, 0, 0) → (0, 1, 1, 0) → (0, 0, 1, 0) → (0, 0, 0, 0)
− |d− DEN V | ≤ |dT LR3 |
II.
+ − + − − d+ IF N + dDEN V + |dT LR3 | ≤ dSOCS + |dSOCS | + |dDEN V |
III.
− − + − d+ T LR3 + |dIF N | + |dDEN V | ≤ dSOCS + |dSOCS |
IV.
+ − − d+ SOCS ≤ dT LR3 + |dIF N | + |dDEN V |
V.
− + + − d+ SOCS + |dSOCS | ≤ dIF N + dDEN V + |dT LR3 |
VI.
− + d+ T LR3 + |dT LR3 | ≤ dSOCS
VII.
− + − + d+ SOCS + |dSOCS | ≤ dT LR3 + |dT LR3 | + dDEN V
VIII.
− + − d+ DEN V + |dT LR3 | ≤ dSOCS + |dSOCS |
us
cr
I.
M
an
Qualitative Cycle
one, the invariance kernel was computed only for this cycle, detailing the delay constraints governing the trajectory. The table also allows us to list the relations (≤, , =) between the individual delays in
d
405
Ac ce pt e
the form of a matrix (given in Table 4) computed by using a linear constraint solver over finite ranges for each delay (the computed sets are provided in Supplementary File 2).
4. Discussion
Dengue infection is the form of critical disease that shows the symptoms even before the activation 410
of adaptive immune response [76]. This infection also presents viral particles to various organs, such as liver, spleen, lymph nodes, brain and spinal cord by inhibiting the IFNα/β mediated signalling pathways [78]. Previous studies showed that DENV uses its NS proteins to down regulated the IFN induction, which is very necessary to control the DENV infection in early phases of its progression [37]. More than 300 IFN stimulated genes are induced as a result of the activation of Janus Kinase (JAK), Tyrosine Kinase
415
(Tyk) signal transducers and activation of transcription (STAT)1/2 and interferon regulatory factor (IRF) when IFNs bind to cell surface receptors [14]. For the clearance of DENV particularly in early stages of infection, TLR3 mediated IFN production is very critical in determining course of pathogenesis. However, uncontrolled TLR signalling results in chronic inflammation, infectious disease propagation and autoimmune 18
Page 18 of 41
Table 4: Relation Matrix of the significant cycle which depicts binary relations between the states. Each entry of the matrix represents whether a delay ’a’ would be greater than, equal to, or less than a delay ’b’. The matrix is read as drow ∼ dcolumn where ∼ ∈ {≤, , ≥, =}.
Relation Matrix |d− T LR3 |
d+ IF N
|d− IF N |
d+ SOCS
|d− SOCS |
d+ DEN V
d+ T LR3
=
≥, ≤
≥, ≤
≥, ≤
≥
≥, ≤
=
|d− SOCS |
≥, ≤
≥, ≤
≥, ≤
≥, ≤
≥, ≤
d+ DEN V
≥, ≤
≥, ≤
≥, ≤
≥, ≤
≥, ≤
|d− DEN V |
≥, ≤
≤
≥, ≤
≥, ≤
cr
us
≥, ≤
>
=
≤
≥, ≤
≥
=
≥, ≤
≥, ≤
≥, ≤
=
an
≥, ≤
M
