Figures
Abstract
The discrete modeling formalism of René Thomas is a well known approach for the modeling and analysis of Biological Regulatory Networks (BRNs). This formalism uses a set of parameters which reflect the dynamics of the BRN under study. These parameters are initially unknown but may be deduced from the appropriately chosen observed dynamics of a BRN. The discrete model can be further enriched by using the model checking tool HyTech along with delay parameters. This paves the way to accurately analyse a BRN and to make predictions about critical trajectories which lead to a normal or diseased response. In this paper, we apply the formal discrete and hybrid (discrete and continuous) modeling approaches to characterize behavior of the BRN associated with MyD88-adapter-like (MAL) – a key protein involved with innate immune response to infections. In order to demonstrate the practical effectiveness of our current work, different trajectories and corresponding conditions that may lead to the development of cerebral malaria (CM) are identified. Our results suggest that the system converges towards hyperinflammation if Bruton's tyrosine kinase (BTK) remains constitutively active along with pre-existing high cytokine levels which may play an important role in CM pathogenesis.
Citation: Ahmad J, Niazi U, Mansoor S, Siddique U, Bibby J (2012) Formal Modeling and Analysis of the MAL-Associated Biological Regulatory Network: Insight into Cerebral Malaria. PLoS ONE 7(3): e33532. https://doi.org/10.1371/journal.pone.0033532
Editor: Szabolcs Semsey, Niels Bohr Institute, Denmark
Received: November 15, 2011; Accepted: February 10, 2012; Published: March 30, 2012
Copyright: © 2012 Ahmad et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Funding: This work was supported by research grants from Higher Education Commission of Pakistan (Grant No. 20-1464/R&D/09/5252) and COMSTECH-TWAS (Grant No. 09-047RG/ITC/AS_C) (Standing Committee on Scientific and Technological Cooperation (COMSTECH) and The Academy of Sciences for the Developing World (TWAS)). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Severe forms of malaria claim a huge number of lives worldwide, contributing to over a million deaths annually, mostly that of children [1]. Malaria is vectored by Anopheles mosquitoes and is a common infectious disease caused by Plasmodium parasites that readily infect blood erythrocytes [2]. In a few cases a severe pathogenesis occurs due to hyperinflammation, usually following Plasmodium falciparum infection, that may turn fatal. The blood flow through small blood vessels to the brain is severely hampered as the infected erythrocytes are sequestered by parasites causing ischaemic hypoxia and increased nitric oxide production in brain tissues, leading to coma, a condition known as diffuse encephalopathy or CM [3]–[5]. Clinical prognosis depends on factors ranging from patient's social conditions to recurrent parasitic exposure, however current evidence strongly suggests that the genetics of an individual may play an even more significant role [6]–[8].
The human innate immune system is the first line of defence against such infections and the responses include inflammation which helps to control the infection and promotes healing. Yet if left unchecked, this advantageous inflammatory response turns astray, causing effects ranging from mild allergies to severe inflammatory disorders [9]. Acute inflammatory response is initiated following an infection through the production of proinflammatory cytokines, such as the tumor necrosis factor alpha (TNF-) and the interferon gamma (INF-
) that play a prominent role in parasite destruction. The generation of inflammation is tightly regulated at multiple levels to control this production. However, in cases where pathogenesis becomes severe, chronic over production of cytokines contributes to elevated levels of a cellular messenger, induced nitric oxide synthase (iNOS). These elevated levels of iNOS plus the hypoxia caused by the parasites work in sync to create a condition of chronic hyperinflammation causing an augmentation of CM pathogenesis [3], [10].
The signal transduction pathway involved in systemic production of proinflammatory cytokines in case of malaria is initiated following the activation of the Toll like receptor 2 (TLR2) and TLR4 when they recognize glycosylphosphatidylinositols (GPIs) anchored on plasmodium membrane proteins [11]–[13]. TLRs are characteristic type I transmembrane pattern recognition receptors, used by the innate immune system to recognize conserved microbial structures or pathogen-associated molecular patterns (PAMPs). They have a conserved cytoplasmic toll-interleukin 1 receptor (TIR) domain and are included in the interleukin 1 receptor (IL-1R)/TLR super-family [14]. After stimulation by PAMPs, TLRs form dimers and begin an intricate multifaceted signalling cascade, initiated by the recruitment of adapter proteins at their cytoplasmic TIR domain. The myeloid differentiation primary response protein (MyD88) and a number of kinases (like interleukin-1 receptor-associated kinases; IRAKs) are recruited downstream. The signal transduction culminates with the activation of nuclear factor kappa-light-chain-enhancer of activated B cells (NF-B), resulting in the expression of proinflammatory cytokine genes (see Figure 1) . For details of the TLR signalling see reference [15].
The TLR2/4 signalling pathway starts with recognition (1) of PAMPs by TLRs. This activates (2) BTK which phosphorylates (3) MAL. MyD88 adapter protein and kinases are recruited (4) and activated around MAL. This eventually leads to the activation (5) of NF-B as I
B is degraded. The proinflammatory cytokine genes are expressed (6) producing INCY that are secreted (7). INCY are responsible (8) for the production of inflammation and activation of their respective receptors. This again activates NF-
B (9a) and through an alternate pathway induces the production of SOCS-1 (9b). SOCS-1 negatively regulates MAL by polyubiquitination (10a) and blocks NF-
B mediated expression (10b). Abbreviations: TLR, toll like receptors; PAMPs, pathogen associated molecular patterns; BTK, bruton's tyrosine kinase; MAL, MyD88 adapter like; MyD88, myeloid differentiation response protein; NF-
B, nuclear factor kappa-light-chain-enhancer of activated B cells; I
B, inhibitor of
B; INCY, proinflammatory cytokines; SOCS-1, suppressor of cytokine signaling-1 [15], [24], [25].
Khor et al. and others [8], [16]–[18] indicate the toll-interleukin 1 receptor domain containing adapter protein (TIRAP), also known as MyD88-adapter-like (MAL), as a key adapter protein associated with disease pathogenesis leading to hyperinflammation. They argue on the basis of single nucleotide polymorphism (SNP) analysis that wild-type MAL may lead to chronic hyperinflammation, increasing host susceptibility towards developing CM; while mutants may confer a protective effect from CM pathogenesis for its carrier host by controlling the amount of inflammation produced. Some studies have suggested MAL as an important therapeutic target for the management of CM pathogenesis [19]–[21].
Signalling regulators associated with MAL include Bruton's tyrosine kinase (BTK) which is a member of the Tec family of kinases. It acts as a positive regulator of MAL by activating it through phosphorylation [22], [23]. The inflammatory cytokines (INCY) produced not only generate inflammation but also induce SOCS-1 (suppressor of cytokine signaling-1), which degrades phosphorylated MAL through polyubiquitination. It is also known to block NF-B mediated expression, thereby functioning as a negative regulator for this signalling pathway [24], [25].
A biological regulatory network (BRN) is a set of interactions (activation and inhibition) between biological entities (e.g., proteins in a biological signalling network). The MAL associated BRN (Figure 2) was abstracted from the TLR2/4 signalling pathway (Figure 1) and includes: BTK that activates MAL through phosphorylation; SOCS-1 that degrades phosphorylated MAL by polyubiquitination; NF-B which is the key transcription factor for initiating the expression of proinflammatory genes; and the INCY which produce the inflammation.
Numerals (1 and 2) represent the threshold levels of interactions; plus (+) signs indicate activation while inactivation is indicated by a minus (−) sign. Arrows indicate the direction of activation/inactivation. The thresholds values are set according to Definition 2.
Two tools, GENOTECH [26] and HyTech [27], were used for the modeling and analysis of this BRN. In GENOTECH tool the qualitative modeling formalism of René Thomas [28]–[36] has been implemented. A BRN in GENOTECH is represented by a directed graph where nodes represent biological entities and edges represent interactions among them. The edges are labeled with integers representing concentration thresholds and signs of interactions (+for activation and −for inhibition). Each entity is assigned a set of logical parameters to generate a state graph (qualitative model) reflecting the possible steady state behaviors of the BRN. GENOTECH also generates a hybrid model by incorporating clocks that measure activation and degradation delays. This hybrid model is analyzed by HyTech which is a tool for the analysis of hybrid systems [37] to compute parametric constraints for the existence of a particular behavior.
This study uses qualitative discrete and hybrid modeling with delays to model the MAL associated BRN. The model characterizes conditions when the system may converge toward hyperinflammation. Initially, a discrete qualitative modeling of the BRN was done to predict qualitative behaviors or states which either lead to a normal (normal inflammatory) or diseased (hyperinflammatory) state. The qualitative model was then refined with clocks measuring regulation delays to obtain a hybrid model incorporating discrete and continuous behaviors. Our results suggest that the system converges towards hyperinflammation when concentration of cytokines is at high level. Furthermore, we obtained the conditions for each behavior of the MAL associated BRN in terms of delay constraints in order to select consistent models regarding delays (see HyTech results for delay constraints). These results give useful information such as, which variable is evolving with high production or degradation rate in order to follow a particular behavior or path. The delay constraints further suggest that for some given set of delay values, the behaviors (see Discussion section) would be very stable and for other set of values they would not exist.
In summary, we make use of qualitative in-silico modeling; a kind of computational method [38] used in systems biology, to characterize the behaviors of MAL associated BRN. At the same time we follow the trend in systems biology [39] to augment the descriptive and analytical power of models in systems biology by formal methods.
Results
This section presents the results obtained by using the modeling and analysis tools for the MAL associated BRN (Figure 2). The first subsection describes the results obtained by the qualitative modeling and analysis of the BRN using the GENOTECH tool [26]. The results obtained by using the HyTech [27] tool for the hybrid modeling and analysis of the BRN are given in the second subsection.
The State Graph of the MAL associated BRN
GENOTECH is a tool for the qualitative (discrete) modeling of BRNs according to Thomas' formalism (see methods). It takes a BRN and corresponding logical parameters from a user in a simple graphical user interface (GUI) and then automatically generates the state graph in which stable states, cycles and paths between any two states can be identified. GENOTECH generates the state graph (qualitative model) of the MAL associated BRN (Figure 2) representing all possible transitions from one state to the other (Figure 3), where each state indicates the concentrations of every entity at a particular time. A state is represented by BTK, MAL, NF-B, INCY and SOCS-1, in the respective order. In a stable state, the whole system converges, halts and cannot proceed to any other state. The stable states are called sinks and the network from any initial state either moves towards these sinks, or ends up in a cyclic path if there is a cycle in the network. From any state other than the stable states, a path that leads to and culminates at a sink is called a trajectory; and all trajectories culminate at some sink or remain in a cycle. It is an inherent property of any system to attain stability; therefore whenever the system is perturbed from its stable state or sink, it tends to shift back or reach another stable state. For the functional and realistic model of MAL-associated BRN, the generation of a specific stable state where every component has a zero concentration is imperative as this state represents the system before any perturbation ( also called the virgin state).
The complete state graph is obtained by using the GENOTECH tool. Definitions 2 and 4 assist in setting the values of the K-parameters. The K-parameters are set such that they result in a model coherent with the observed steady states behaviors. In the case of MAL-associated BRN these states are (0,0,0,0,0) and (0,0,1,2,1). In the state (0,0,0,0,0) the system does nothing when there is no signal of pathogen or returns to this state after proper response. In the state (0,0,1,2,1), the system produces inflammation continuously. The set of logical parameters is: ,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
and
.
In our case two stable states (also called stable steady states or sinks) are produced, indicating zero concentrations for every component in the first state; and zero concentration for BTK and MAL, 1 for NF-
B,
for INCY and
for SOCS-1
in the second state. The stable state
represents the normal behavior where infection is cleared whereas
represents the diseased behavior leading to hyperinflammation. Since BTK is the initiator for this BRN, the state
was taken as the starting state. All possible trajectories from this starting state leading to the normal and diseased states were identified (Figure 4). The trajectories leading to the sink
are representatives of a normal behavior against infection, where inflammation is produced to clear the infection and diminished afterwards.
Each circle represents a particular state (configuration) and inside the circle the values 0,1 and 2 represent qualitative levels of proteins according to the order (BTK, MAL, NF-kB, INCY, SOCS-1). The solid lines show the transition which leads to normal state (0,0,0,0,0) and the dotted lines show the transitions towards the state of hyperinflammation (0,0,1,2,1). The conditions for all the transitions in this figure are given in Table 1.
HyTech Results for Delay Constraints
All possible transitions of normal (towards (0,0,0,0,0)) and divergent (towards (0,0,1,2,1) pathways are highlighted in Figure 4 and the corresponding conditions are listed in Table 1. The state (1,0,0,0,0) is the initial state of the hybrid model. There are two possible outer transitions which are towards (0,0,0,0,0) and (1,1,0,0,0) from the initial state. The constraint shows that if the deactivation (or degradation) delay of BTK is less than or equal to the activation delay of MAL, then the next location of the system will be (0,0,0,0,0), which is a stable state representing normal behavior. The transition
(Figure 4 and Table 1 - Transition 1) can only occur if the positive delay of MAL (
) is less than or equal to the degradation delay of BTK (
), i.e. if the rate of activation of MAL is faster than the rate of degradation of BTK. The transition
can only occur if the sum of the activation delays of MAL and NF-
B will be greater than or equal to the degradation delay of BTK (see transition
in Table 1). The transition
will occur unconditionally because there is only one outer transition from state (0,1,0,0,0) and this will occur after some delay
which is needed to activate NF
B. The transitions which lead to divergence from normal pathways (shown as dotted lines in Figure 4) are important because these can provide useful information about the BRN dynamics. The state (0,1,1,1,0) has two transitions, one towards a normal and the other towards a divergent path (see Figure 4, transitions 5 and 15). The transition
is only possible if SOCS-1 becomes active before the cytokines reaches level
(INCY2). If this is not the case, then the network will be in state (0,1,1,2,0), which ultimately leads to the diseased state (0,0,1,2,1). Another very sensitive transition is
, which directly leads into a diseased state due to the high level of cytokines (level
). In order to follow a particular path all conditions (see Table 1) of individual transitions must be satisfied. The conditions for all paths leading to the stable states are given in Table 1 and 2 respectively. All the constraints of Table 1 and 2 were computed in
seconds on Intel Core i7 machine.
Discussion
MAL is imperative for TLR2/4 downstream signaling, acting both as an adapter sorter by bringing in other adapter proteins to the TLR cytoplasmic TIR domain. It also acts as a site for the assembly of several kinases important for downstream signal transduction leading to inflammation. Although production of inflammation in response to infection is not just limited to the TLR pathway, studies by Khor et al. and others have indicated a putative role of MAL in the pathogenesis of CM due to hyperinflammation [8], [16]–[18]. In this study we use discrete qualitative and hybrid modeling with delay constraints to characterize the behavior of MAL associated BRN.
Discrete Modeling
A qualitative discrete model based on the kinetic logic of Thomas is obtained using GENOTECH by providing it information of component interactions and their respective threshold concentration levels. The GENOTECH results are rendered as a state transition graph (Figure 3) referred to as the MAL state space from here on and it represents different configurations the system may reside in as an inflammatory response is initiated. The two stable states or sinks generated indicate two distinct behaviors for the system, parallel to the normal and diseased behaviors, and any configuration of the system may eventually move toward either of these two sinks through specific trajectories or paths (Figure 3 and 4). These two states are (0,0,0,0,0), indicating zero concentration for every component; and (0,0,1,2,1), indicating zero concentration for BTK and MAL, one for NF-B, two for INCY and one for SOCS-1 respectively. The sink (0,0,0,0,0) represents the normal behavior where inflammation is produced in response to infection and the system reverts to this resting state after mounting an inflammatory response. In contrast, the sink (0,0,1,2,1) represents a diseased behavior where INCY levels are very high and a condition of hyperinflammation is produced. MAL state space allows for transitions between different states forming trajectories or paths that lead toward either the (0,0,0,0,0) or the (0,0,1,2,1) sink and the significance of these transitions is that they represent the indigenous ability of the system to attain stability.
The complete state space is very large and it becomes important to isolate key trajectories that lead either to a normal or diseased sink. The shortest possible route toward the normal sink was identified and used as the reference trajectory (Figure 4). Small divergences from this trajectory that culminate at the normal state were also identified. Other critical divergences were those that lead to the diseased state and these present the perturbations the system may encounter. Divergences that eventually lead to the normal state indicate perturbations that the system can tolerate, while divergences that can only lead to the diseased state are perturbations that put too much stress on the system resulting in a collapse of the system controls.
This MAL state space suggests that the system maintains the capacity to stabilize itself after perturbations as long as the concentration of INCY for any particular state within the system remains at level 1. But once the concentration of INCY attains level 2, the system enters trajectories that lead only to the sink (0,0,1,2,1), representing the diseased behavior. It is important to note that the threshold concentration level for INCY for the diseased state is 2. This indication is supported by the experimental observations made in the case of bacterial sepsis where the pathogenesis is associated with the over production of INCY and hyperinflammation [40], [41]. The pathogenesis and tissue injury in case of sepsis is closely related to the severe malarial pathogenesis [42].
Another important indication is that the moment INCY appears, the system is at risk of entering trajectories that may only lead to the diseased state (Figure 4, transitions 23 and 24). For the normal trajectory the INCY concentration level is maintained at 1, however when concentration level of 2 for INCY is acquired, the system loses its capacity to revert toward normality and moves towards the diseased state in a progressive manner. Furthermore, since INCY are produced through several alternate pathways, the model suggests that if INCY levels are already elevated due to some other reason, like prior inflammatory condition, the chance of the system shifting towards diseased state elevates. The state (0,0,1,1,1) is indicative of this where it culminates directly at the diseased state with only one transition (Figure 4, transition 24). Moreover, the model also indicates that an important component for this BRN (apart from MAL) is BTK which remains constitutively active in the cases where the divergence from the normal trajectory leads to the diseased state (Figure 4, for example the trajectory with transitions 1, 10, 12, 14, 22, 18 and 19). BTK is important for activation of MAL through phosphorylation, and this suggests that BTK may also play an important role in pushing the system towards hyperinflammation, and warrants further investigation.
Hybrid Modeling
The rate at which a component of a BRN is regulated in relation to other components determines the direction of the evolution of the system. These rates are proportional to the threshold concentrations for each component and are modeled as unvalued parameters with delays. HyTech [27] uses these delays as unvalued parameters along with the discrete model of our BRN to calculate the conditions, in the form of constraints that are required for a transition between two states to occur. A path in the model is a sequence of states leading to a particular stable steady state and HyTech synthesizes the conditions for transitions between the states (Table 1). The conjunctions of all the constraints of the path constitute the conditions for that path to be followed (Table 2). Only when these conditions are satisfied, the system will reach a stable state representing either normal or diseased behavior.
(1,1,1,1,0) (1,1,1,2,0) & (0,1,1,1,0)
(0,1,1,2,0) (transitions 16 & 23): These transitions require that the rate of production of INCY in order to reach high level (level 2) should be faster than the induction of SOCS-1, suggesting that reducing high cytokine levels and inducing SOCS-1 can help to avoid this transition.
(1,1,1,1,1) (1,1,1,2,1) (transition 22): The conditions for this transition suggest that if BTK remains active for a long period of time, the system will lead to the diseased state. This can be due to continuous activation of MAL by BTK and adapting strategies to prevent this is suggested.
(0,1,1,1,1) (0,1,1,2,1) & (0,0,1,1,1)
(0,0,1,2,1) (transitions 26 & 24): These transitions suggests that the system may also move towards diseased state due to multiple factors being involved that include: high rate of production of INCY; slow rate of deactivation of MAL, induction of SOCS-1, and deactivation of NF-
B. To avoid these transitions, intervention at multiple levels is suggested that include: reducing high cytokine levels, NF-
B and their effects; faster rate of deactivation of MAL which in turn is affected by higher rate of SOCS-1 induction; slower rate of activation of MAL being controlled by BTK.
The conditions (Table 1) for transition between these states imply that, to avoid hyperinflammation, strategies should be designed that: i) inhibit continuous activation of MAL by BTK, ii) reduce the effect of high cytokine levels and iii) induce SOCS-1. Our results suggest that under certain conditions, MAL related inflammatory response in malaria patients, may cause hyperinflammation leading to CM. We hypothesize that constitutively active BTK keeps activating MAL to produce the diseased behavior and may be an important drug target. BTK, a tyrosine kinase enzyme, is an attractive drug target and inhibitors are already being developed against it for a number of diseases like B-cell lymphomas and some autoimmune disorders [43].
The rate at which INCY are produced within the system is crucial and will determine its fate [44]. These key divergence states have INCY at level 1 and the system is at risk of going into the diseased state if these levels are further increased. Individuals can have high INCY levels due to various reasons which may include genetic predisposition, repeated exposure to infectious agents or chronic inflammatory conditions. Such conditions may partially explain prevalence of CM in some Sub Saharan African populations where malaria is endemic all year round and cytokine levels are high due to asymptomatic parasitaemia [45]–[49]. The effects of high cytokine levels may be controlled through the use of anti-inflammatory agents. A recent study by Franklin et al., 2011 [50], emphasizes the relationship between elevated levels of inflammatory cytokines and CM. They report a promising strategy for the management of malarial severity by interfering with the release of cytokines through inactivation of TLRs using antagonists. This supports our second hypothesis that suggests an association between a preexisting high level of cytokines and the progression of malaria towards CM. We suggest that such individuals or population groups may be more susceptible to develop CM and this opens up new avenues for future investigations.
Conclusions
The normal (inflammatory) and abnormal behaviors (hyperinflammatory), presented by this qualitative model are consistent with experimental observations. They are commonly observed among population groups where in some the malarial infection is cleared after a normal inflammatory response, while in others severe malaria (CM) may develop due to hyperinflammation.
This study suggests an association between pre-existing high cytokine levels and progression of malaria towards CM, and opens up new avenues of investigation. This means that screening population groups, in areas where CM is endemic, for inflammatory cytokine levels may help identify ‘at risk’ individuals. This can be used as a predictive indicator that such an individual may progress towards CM. Multiple strategies that interfere with the TLR signal transduction, e.g. the use of BTK inhibitors and antinflammatory agents along with existing treatments should be further explored.
This study also shows that formal methods can support BRN modeling, which in turn facilitates state space exploration and trajectory computation. In future we plan to conduct the formal probabilistic analysis [51] to gain more insight into the MAL signaling network.
Materials and Methods
Qualitative modeling
Biological systems are traditionally modelled with ordinary differential equations that refer to the time derivative of each quantity (concentration, rates and temperature etc.). Due to the inherent complexity of biological systems, the exact amounts of these quantities are rarely known. In 1970, René Thomas [38] introduced a Boolean logic based method for the modeling of biological regulatory networks in which the system dynamics are modelled qualitatively. He proved the effectiveness of his modeling approach by analyzing the lambda phage genetic switch [38], [52]–[55]. According to René Thomas, the Boolean modeling suffers from some limitations because it has only two levels, or
, which are not sufficient to formalize all kinds of relevant problems. Later on, Thomas generalized the Boolean logic to ‘kinetic logic’ [28]–[36] and showed its practical effectiveness by applying it to different genetic regulatory systems. An interesting feature of ‘kinetic logic’ is that it closely approximates the differential equations based modeling [29].
The Kinetic Logic formalism of René Thomas
The kinetic logic formalism models BRNs by focussing on threshold effects. The entities (usually proteins) of biological systems can interact with each other either positively or negatively; that is the concentration level of one element can increase or decrease the rate of activation of other elements. In our setting activation means production of a particular protein, if it is absent, or its activation from an inactive state. Similarly, inactivation means inhibition of a particular protein or its degradation. Following are the two types of biological regulation which are of sigmoid (nonlinear bounded curve) nature.
Activation
If a product increases the rate of activation of product
, it is a positive regulator. In this situation the rate of activation of
increases with increasing concentration of
, which can be depicted by a sigmoidal representation as shown in Figure 5. It can be seen that there is little effect of
on the rate of activation of
for as long as it remains below a given threshold,
. Once
reaches this threshold, the rate of activation of
increases rapidly until
is saturated. In other words we can say that
was ‘inactive’ when
and ‘active’ for
, this suggests the approximation of the sigmoidal curve by a step function as shown in Figure 5(A.).
(A.) Sigmoidal (above) and Step function (below) representation of activation. (B.) Sigmoidal (above) and Step function (below) representation of degradation (inhibition).
Inhibition
If product reduces the rate of activation of product
, it is called a negative regulator or inhibitor. The effect of this negative regulation or inhibition is shown in Figure 5(B.).
Semantics of the Kinetic Logic Formalism
This section presents the semantics of Thomas's formalism with a hypothetical running example of a gene regulatory network involving three genes (x, y and z).
Definition 1 (Directed Graph). A directed graph is an ordered pair , where
is the set of all nodes and
is the set of ordered pairs called edges or arcs
An edge is considered to be directed from
to
;
is called the head and
is called the tail of the edge. In a directed graph
,
and
denote the set of predecessors and successors of a node
, respectively.
Figure 6 shows a network represented by a directed graph showing the genes in which and
are the sets of nodes and edges, respectively. In case of gene
, (The set of regulators of
)
and
=
represents the predecessors and successors of
.
A toy example of a BRN where ,
and
represent biological entities. The labels ‘+’ and ‘1’ represent the activation and the threshold concentration respectively.
Definition 2 (Biological Regulatory Network). A biological regulatory network (BRN) is a labeled directed graph , where
is a set of nodes which represents biological entities and
, is a set of all possible edges which represent the interaction between entities.
- Each edge
is labeled by a pair
, where
is positive integer (qualitative level representing a threshold) and
is either the (+) sign or the (−) sign, where (+) represents activation and (−) inhibition. For an example see Figure 6, where
,
,
and similarly
,
,
.
- Each node
has a limit
, which is equal to out-degree of
(the total number of targets of x), such that
each
where
.
- Each entity
has its abstract concentration in the set
.
It is important to consider the possible number of states and transitions between them to understand the behavior of the BRN.
Definition 3 (States). Let be a BRN. The state of a BRN is a tuple
, where
The qualitative states are represented by vector , where
denotes the level of concentration of product
. A qualitative state represents a configuration of all the elements of a BRN at any instant of time. The number of activators of a particular variable at a given level of concentration are represented by its set of resources (see the definition of Resources given below). In the following definition, we formally define the set of resources which represents all the activators of an entity at any instance of time.
Definition 4 (Resources). Let be a BRN. The set of resources
of a variable
at a level
is defined as
=
and
or
and
.
The dynamic behaviors of BRN depends on logical parameters. The set of these logical parameters is defined as .
The parameter (at a level
of
) gives the information about the evolution of
. There are three cases: 1) if
then
increases by one unit, 2) if
then
decreases by one unit, and 3), if
then
cannot evolve from its current level. In case of the running example of three genes (Figure 6), the set of resources and logical parameters are shown in Table 3.
It is convenient to describe the evolution from one level to the other by an evolution operator <$>\raster="rg1"<$> [56], which is defined as follows:
where
and
Definition 5 (State Graph). Let be a BRN and
represents the concentration level of entity
in a state
. Then the state graph of the BRN will be the directed graph
, where
is set of states and
represents a relation between states, called the transition relation, such that
if and only if:
a unique
such that
and
, and
.
From the table of resources (Table 3), interesting information regarding its dynamic behavior can be derived, c.f. Figure 7. It is clear from the state graph that there are two stable states, namely (000) and (111).
The state graph of the BRN in Figure 3. Each node represents a qualitative (discrete) state of the BRN. The values inside a state represent the concentrations of the entities ,
and
.
Discrete Modeling of MAL associated BRN using GENOTECH
GENOTECH [26] facilitates modeling of BRNs according to Thomas' formalism. It takes a BRN and the corresponding logical parameters form a user in a simple graphical user interface (GUI). It automatically produces the whole state space in which stable states, cycles and paths between any two states can be identified. The software tool GINsim (Gene Interaction Network simulation) [57] could also be used for the qualitative modeling and analysis of BRNs. GENOTECH is available for download at http://code.google.com/p/genotech/downloads/list and the following step are required for modeling a BRN in GENOTECH.
- Construction of a BRN as a labeled directed graph: The GUI of GENOTECH provides a set of commands to create, edit and save a BRN. The two drop down menus Gene and Interaction are used to construct the nodes and edges of a BRN respectively; and each biological entity is assigned a set of logical parameters (see Figure 8 and 9).
- Generation of a state graph: After the construction of the BRN, a state graph can be generated by using the command state graph in the File menu. The state graph, showing stable states in red, appears in a new window with its own set of commands to find cycles, paths between two states, and neighboring states (see Figure 10). GENOTECH also provides an option to save the state graph in DOT format [58] for visualization in Graphviz tool [59].
The BRN is constructed as a directed graph by using the Gene’ New and Interaction’New menu options. An edge between two entities shows an interaction which is labeled with threshold and the sign of interaction (+for activation and −for degradation).
Each entity of the BRN is assigned a set of logical parameters by using the properties option accessible by right clicking that entity.
The stable states are highlighted in red. The right panel shows the analysis command which include: show path which highlights paths between two selected states; show cycle which highlights existing cycles (closed path), if any, in the state graph; show neighboring states highlights all the neighboring states of the selected state. The conversion menu contains the commands to export the graph to DOT and HyTech formats.
The complexity of a state graph increases depending on the number of entities, their interactions, thresholds and logical parameters. In order to make this method more scalable, a simplification technique based on the coloration method, which extracts the desirable parts of the state graph, has been proposed in our previous work [60] – the future versions of GENOTECH will implement this technique.
Figure 1 shows the TLR2/4 molecular pathway adapted from literature [15], [24], [25]. The MAL associated BRN was abstracted from this pathway and is shown in Figure 2. The Network flow is as follows: BTK activates MAL and MAL activates NF-B after reaching threshold level of 1; INCY are activated by NF-
B at a threshold level of 1; there are two interactions of INCY which are the activations of SOCS-1 and NF-
B at level 1 and level 2 respectively; and SOCS-1 inhibits MAL and NF-
B at threshold level 1. Using the Thomas' formalism presented in the previous section, preliminary insilico experiments were performed to set the threshold orders and to determine the values of logical parameters (see Figure 3) in order to make the model consistent with naturally observed behaviors (steady states).
Hybrid modeling to enhance the discrete modeling with delays
A hybrid model of a BRN involves discrete locations associated with some continuous variables (clocks). The continuous transitions represent time that elapses in a location, whereas discrete transitions show the instantaneous change between locations [61]. The timing diagram of one normal path is shown in Figure 11. It depicts proteins evolving from one level to another (from x to x+1 or x+1 to x) in a discrete fashion (Figure 12). In reality, however, the concentrations of these proteins evolve in a nonlinear and continuous manner and this behavior cannot be represented in a discrete modeling framework. Various formalisms [62], [63] have been proposed for biological modeling to overcome this limitation of discrete modeling. Ahmad et al. [64] proposed the refinement of discrete modeling by hybrid modeling, where sigmoidal nature evolutions are modelled with piecewise linear curves. This contrasts discrete modeling, in which these evolutions are modelled using step functions (Figure 12). Considering that a delay is required for the evolution of a protein from level to
, or from
to
, it is important to introduce some additional concepts, namely time intervals and clocks.
Timing diagram showing the evolution of proteins involved in BRN. Here the concentrations are represented by two levels x and x+1 on the vertical axis. The horizontal axis represent the time of evolution. The dotted lines are the boundaries of the different configurations of the discrete concentrations. Each configuration is according to the order (BTK , MAL, NF-kB, INCY, SOCS-1).
Evolutions of an entity is shown as: sigmoidal representation of the activation (A.), discrete approximation of the activation (B.), piece-wise linear approximation of the activation (C.), sigmoidal representation of the degradation (D.), discrete approximation of the degradation (E.) and piece-wise linear approximation of the degradation (F.).
The approach of Siebert and Bockmayr [63] is similar to our approach. The authors incorporate time delays in the logical formalism of René Thomas using timed automaton [65] which is a more restrictively expressive model than hybrid automaton, since it does not deal with non-increasing variables.
Clocks are continuous variables used in timed automaton based models [65], which are a subclass of hybrid automata [66]. Each protein is associated with a clock variable that synchronously evolves with time. These clock intervals reflect the characteristics of continuous dynamics within the available discrete formalism [61], [64]. The time measured by the clock variable
between two levels is called the delay between these levels. The initial value of the clock variable is set to zero and when the value of this clock variable is equal to delay time
or
the transition between two levels takes place. The delays
or
represents time taken from
to
(positive delay or production delay) or
to
(negative delay or degradation delay), respectively (Figure 12). The speed by which the clock variables
evolve is modelled by the first order derivatives
, where
lies in the set
. This value characterizes the evolution of the associated variable, which normally represents the concentration level of a protein [61], [64].
We thus obtain a hybrid model that is suitable to represent the discrete and continuous dynamics of the systems that we consider. In this work, delays are considered as unvalued parameters, which motivate the introduction of Parametric Bio Linear Hybrid Automata (Bio-LHA) for the modeling of the MAL associated BRN, as presented in the next subsection.
Parametric Bio-LHA
Parametric Bio-LHA was originally proposed by Ahmad et al. [26] for the linear hybrid modeling of BRN. This framework refines the discrete (qualitative) model by incorporating clocks and time delays in it. The resultant model is then suitable for the computation of exact conditions in the form of delay constraints for the existence of the behaviors of BRN.
Let (resp.
) be the set of constraints using only
(resp.
) where
and
are the sets of real valued variables and parameters respectively.
Definition 6 (Parametric Bio-LHA). A parametric Bio-LHA is a tuple
where
is a finite set of locations.
is the initial location.
is a finite set of parameters (delays).
is a finite set real-valued variables (clocks).
is a finite set of edges,
represents an edge from the location
to the location
with the guard
and the reset set
; we require that the set of clocks in
is a subset of
.
assigns an invariant to any location.
maps each pair
to an evolution rate.
The semantics of a parametric Bio-LHA is a timed transition system. We define the semantics according to the time domain . We let
.
Definition 7 (Semantics of Bio-LHA). Let be a valuation for the parameters
and
represents the values of clocks in a location. The
-semantics of a parametric Bio-LHA
is defined as a timed transition system
where: (1)
; (2)
is the initial state and (3) the relation
is defined for
as:
-
discrete transitions:
iff
such that
,
if
and
if
.
-
continuous transitions: For
,
iff
,
, and for every
,
.
The partial Bio-LHA of the MAL associated BRN is shown in Figure 13. For the sake of simplicity, only three locations are considered in this model. The inequalities, such as (), represent invariants or the stay conditions for the BRN to remain in a particular configuration (location). For example, (
) requires that the clock variable
of BTK should be less than or equal to the degradation delay (time required to change level
to
) of BTK.
Partial view of Bio-LHA of the MAL associated BRN. The labels 11000, 01000 and 11100 shows the state (configuration) of BRN consisting of five entities (BTK, MAL, NF-kB, INCY, SOCS-1).
The rate dhBTK (HyTech representation of dhBTK/dt) provides an indication about the evolution from one level to another. The clock variables (i-th protein) measure the time taken by associated proteins in a particular level, and clocks are set to zero once the transition has occurred. The complete Bio-LHA involves forty eight locations and five clock variables associated with each protein involved in the BRN (Figure 2).
The next section presents the investigation of causality conditions for the normal and the divergent behaviors.
Synthesis of delay parameters using HyTech: A Linear Hybrid Model checking approach
In this study, a Bio-LHA was used for the hybrid modeling of the MAL associated BRN. As discussed above, the unvalued delays represent parameters in the Bio-LHA. In order to analyze this model, the linear hybrid model checking tool HyTech was used, which is widely applied to the modeling and verification of hybrid systems. The main motivation for using HyTech for the analysis of this class of systems is its very rich set of analysis commands as well as the parameter synthesis capabilities of this tool (Table 1 and 2).
Acknowledgments
The authors would like to thank the anonymous reviewers, Stefan Leue and Rehan Zafar for their constructive comments and suggestions to improve the quality of the paper.
Author Contributions
Conceived and designed the experiments: JA UN US. Performed the experiments: JA US. Analyzed the data: SM JA UN. Contributed reagents/materials/analysis tools: JA US. Wrote the background: SM UB JB. Wrote the methods and results sections: JA US. Wrote the discussion and conclusions: SM JA UN.
References
- 1.
Who malaria fact sheet no. 94 world health organization. URL http://www.who.int/mediacentre/factsheets/fs094/en/index.html. Accessed April 2010.
- 2. Miller LH, Baruch DI, Marsh K, Doumbo OK (2002) The pathogenic basis of malaria. Nature 415: 673–679.
- 3. Clark IA, Cowden WB (1999) Why is the pathology of falciparum worse than that of vivax malaria? Parasitol Today (Regul Ed) 15: 458–461.
- 4. Amante FH, Haque A, Stanley AC, Rivera FdeL, Randall LM, et al. (2010) Immune-mediated mechanisms of parasite tissue sequestration during experimental cerebral malaria. J Immunol 185: 3632–3642.
- 5. Hunt NH, Golenser J, Chan-Ling T, Parekh S, Rae C, et al. (2006) Immunopathogenesis of cerebral malaria. Int J Parasitol 36: 569–582.
- 6. Jain V, Armah HB, Tongren JE, Ned RM, Wilson NO, et al. (2008) Plasma IP-10, apoptotic and angiogenic factors associated with fatal cerebral malaria in India. Malar J 7: 83.
- 7. Hirayama K (2002) Genetic factors associated with development of cerebral malaria and _brotic schistosomiasis. Korean J Parasitol 40: 165–172.
- 8. Khor CC, Chapman SJ, Vannberg FO, Dunne A, Murphy C, et al. (2007) A Mal functional variant is associated with protection against invasive pneumococcal disease, bacteremia, malaria and tuberculosis. Nat Genet 39: 523–528.
- 9.
Serhan CN, GDWard PA (2010) Fundamentals of Inammation. Cambridge: Cambridge University Press.
- 10. Artavanis-Tsakonas K, Tongren JE, Riley EM (2003) The war between the malaria parasite and the immune system: immunity, immunoregulation and immunopathology. Clin Exp Immunol 133: 145–152.
- 11. Barton GM, Kagan JC (2009) A cell biological view of Toll-like receptor function: regulation through compartmentalization. Nat Rev Immunol 9: 535–542.
- 12. Krishnegowda G, Hajjar AM, Zhu J, Douglass EJ, Uematsu S, et al. (2005) Induction of proinammatory responses in macrophages by the glycosylphosphatidylinositols of Plasmodium falciparum: cell signaling receptors, glycosylphosphatidylinositol (GPI) structural requirement, and regulation of GPI activity. J Biol Chem 280: 8606–8616.
- 13. Nebl T, De Veer MJ, Schofield L (2005) Stimulation of innate immune responses by malarial glycosylphosphatidylinositol via pattern recognition receptors. Parasitology Suppl 130: 45–62.
- 14. O'Neill LA, Dinarello CA (2000) The IL-1 receptor/toll-like receptor superfamily: crucial receptors for inammation and host defense. Immunol Today 21: 206–209.
- 15. Akira S, Takeda K (2004) Toll-like receptor signalling. Nat Rev Immunol 4: 499–511.
- 16. Hamann L, Kumpf O, Schuring RP, Alpsoy E, Bedu-Addo G, et al. (2009) Low frequency of the TIRAP S180L polymorphism in Africa, and its potential role in malaria, sepsis, and leprosy. BMC Med Genet 10: 65.
- 17. Ferwerda B, Alonso S, Banahan K, McCall MB, Giamarellos-Bourboulis EJ, et al. (2009) Functional and genetic evidence that the Mal/TIRAP allele variant 180L has been selected by providing protection against septic shock. Proc Natl Acad Sci USA 106: 10272–10277.
- 18. Ramasawmy R, Cunha-Neto E, Fae KC, Borba SC, Teixeira PC, et al. (2009) Heterozygosity for the S180L variant of MAL/TIRAP, a gene expressing an adaptor protein in the Toll-like receptor pathway, is associated with lower risk of developing chronic Chagas cardiomyopathy. J Infect Dis 199: 1838–1845.
- 19. Zhu J, Mohan C (2010) Toll-like receptor signaling pathways–therapeutic opportunities. Mediators Inamm 2010: 781235.
- 20. Loiarro M, Ruggiero V, Sette C (2010) Targeting TLR/IL-1R signalling in human diseases. Mediators Inamm 2010: 674363.
- 21. Niazi UH, Bibby J, Sutcli_e MJ (2011) In-silico characterization of the effects of phosphorylated tyrosines 86 and 106 on structure and binding of MAL: insight into hyperinammatory response to infection by the human malaria parasites. J Recept Signal Transduct Res 31: 53–65.
- 22. Gray P, Dunne A, Brikos C, Jefferies CA, Doyle SL, et al. (2006) MyD88 adapter-like (Mal) is phosphorylated by Bruton's tyrosine kinase during TLR2 and TLR4 signal transduction. J Biol Chem 281: 10489–10495.
- 23. Piao W, Song C, Chen H, Wahl LM, Fitzgerald KA, et al. (2008) Tyrosine phosphorylation of MyD88 adapter-like (Mal) is critical for signal transduction and blocked in endotoxin tolerance. J Biol Chem 283: 3109–3119.
- 24. Mansell A, Smith R, Doyle SL, Gray P, Fenner JE, et al. (2006) Suppressor of cytokine signaling 1 negatively regulates Toll-like receptor signaling by mediating Mal degradation. Nat Immunol 7: 148–155.
- 25. Yoshimura A, Naka T, Kubo M (2007) SOCS proteins, cytokine signalling and immune regulation. Nat Rev Immunol 7: 454–465.
- 26.
Ahmad J (2009) Modélisation hybride et analyse des dynamiques des réseaux de régulations biologiques en tenant compte des délais. Ph.D. thesis, Ecole Centrale de Nantes, France.
- 27. Henzinger TA, Ho PH, Wong-Toi H (1997) HYTECH: A model checker for hybrid systems. International Journal on Software Tools for Technology Transfer 1: 110–122.
- 28. Thomas R (1978) Logical analysis of systems comprising feedback loops. Journal of Theoretical Biology 73: 631–656.
- 29.
Thomas R, D'Ari R (1990) Biological Feedback. CRC Press, Boca Raton, FL.
- 30. Thomas R (1991) Regulatory networks seen as asynchronous automata: A logical description. Journal of Theoretical Biology 153: 1–23.
- 31. Snoussi E, Thomas R (1993) Logical identi_cation of all steady states : the concept of feedback loop characteristic states. Bull Math Biol 55: 973–991.
- 32. Thomas R, Thieffry D, Kaufman M (1995) Dynamical behaviour of biological regulatory networks: I. biological role of feedback loops and practical use of the concept of the loop-characteristic state. Bull Math Biol 57: 247–276.
- 33. Thomas R (1998) Laws for the dynamics of regulatory circuits. Int J Dev Biol 42: 479–485.
- 34. Demongeot J, Kaufman M, Thomas R (2000) Positive feedback circuits and memory. CR Acad Sci III 323: 69–79.
- 35. Thomas R, Kaufman M (2001) Multistationarity, the basis of cell differentiation and memory. II. logical analysis of regulatory networks in terms of feedback circuits. Chaos 11: 180–195.
- 36. Thomas R, Kaufman M (2001) Multistationarity, the basis of cell differentiation and memory. I. structural conditions of multistationarity and other nontrivial behavior. Chaos 11: 65–179.
- 37.
Heemels W, Lehmann D, Lunze J, De Schutter B (2009) Introduction to hybrid systems. In: Lunze J, Lamnabhi-Lagarrigue F, editors. Handbook of Hybrid Systems Control – Theory, Tools, Applications. Cambridge, UK: Cambridge University Press, chapter 1. pp. 3–30.
- 38. Thomas R (1973) Boolean formalization of genetic control circuits. J Theor Biol 42: 563–585.
- 39. Priami C (2009) Algorithmic systems biology. Commun ACM 52: 80–88.
- 40. Tsujimoto H, Ono S, Efron PA, Scumpia PO, Moldawer LL, et al. (2008) Role of Toll-like receptors in the development of sepsis. Shock 29: 315–321.
- 41. Klucinski P, Martirosian G (2005) [Role of cytokines and pathogen associated molecular pattern receptors in sepsis]. Przegl Epidemiol 59: 695–701.
- 42. Praetorius F, Altrock G, Blees N, Schuh N, Faulde M (1999) [Imported Anopheles: in the luggage or from the airplane? A case of severe autochthonous malaria tropica near an airport]. Dtsch Med Wochenschr 124: 998–1002.
- 43. Honigberg LA, Smith AM, Sirisawad M, Verner E, Loury D, et al. (2010) The Bruton tyrosine kinase inhibitor PCI-32765 blocks B-cell activation and is efficacious in models of autoimmune disease and B-cell malignancy. Proc Natl Acad Sci USA 107: 13075–13080.
- 44. Clark IA (2007) How TNF was recognized as a key mechanism of disease. Cytokine Growth Factor Rev 18: 335–343.
- 45. Bloland PB, Boriga DA, Ruebush TK, McCormick JB, Roberts JM, et al. (1999) Longitudinal cohort study of the epidemiology of malaria infections in an area of intense malaria transmission II. Descriptive epidemiology of malaria infection and disease among children. Am J Trop Med Hyg 60: 641–648.
- 46. Bottius E, Guanzirolli A, Trape JF, Rogier C, Konate L, et al. (1996) Malaria: even more chronic in nature than previously thought; evidence for subpatent parasitaemia detectable by the polymerase chain reaction. Trans R Soc Trop Med Hyg 90: 15–19.
- 47. Njama-Meya D, Kamya MR, Dorsey G (2004) Asymptomatic parasitaemia as a risk factor for symptomatic malaria in a cohort of Ugandan children. Trop Med Int Health 9: 862–868.
- 48. Smith T, Charlwood JD, Kihonda J, Mwankusye S, Billingsley P, et al. (1993) Absence of seasonal variation in malaria parasitaemia in an area of intense seasonal transmission. Acta Trop 54: 55–72.
- 49. Crookston BT, Alder SC, Boakye I, Merrill RM, Amuasi JH, et al. (2010) Exploring the relationship between chronic undernutrition and asymptomatic malaria in Ghanaian children. Malar J 9: 39.
- 50. Franklin BS, Ishizaka ST, Lamphier M, Gusovsky F, Hansen H, et al. (2011) Therapeutical targeting of nucleic acid-sensing Toll-like receptors prevents experimental cerebral malaria. Proc Natl Acad Sci USA 108: 3689–3694.
- 51.
Calder M, Vyshemirsky V, Gilbert D, Orton R (2005) Analysis of signalling pathways using the prism model checker. pp. 179–190.
- 52. Thomas R (1979) Kinetic logic: a boolean approach to the analysis of complex regulatory systems. Lecture Notes in Biomathematics 29: 507.
- 53. Thomas R, Gathoye AM, Lambert L (1976) A complex control circuit: Regulation of immunity in temperate bacteriophages. Eur J Biochem 71: 211–227.
- 54. Thomas R (1981) On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations. Springer series in Synergetics 9: 180–193.
- 55. Thomas R, Thieffry D (1995) Dynamical behaviour of biological regulatory networks - ii. immunity control in bacteriophage lambda. Bull Math Biol 57: 277–297.
- 56. Bernot G, Cassez F, Comet JP, Delaplace F, Müller C, et al. (2007) Semantics of biological regulatory networks. Electron Notes Theor Comput Sci 180: 3–14.
- 57. Gonzalez A, Naldi A, Sanchez L, Thieffry D, Chaouiya C (2006) GINsim: A software suite for the qualitative modelling, simulation and analysis of regulatory networks. Biosystems 84: 91–100.
- 58.
Koutsofios E, North SC (1993) (1993) Drawing graphs with dot.
- 59. Ellson J, Gansner ER, Koutsofios E, North SC, Woodhull G (2001) Graphviz - Open Source Graph Drawing Tools. Graph Drawing 483–484.
- 60. Ahmad J, Bourdon J, Eveillard D, Fromentin J, Roux O, et al. (2009) Temporal constraints of a gene regulatory network: Refining a qualitative simulation. Biosystems 98: 149–159.
- 61. Ahmad J, Bernot G, Comet JP, Lime D, Roux O (2007) Hybrid modelling and dynamical analysis of gene regulatory networks with delays. ComPlexUs 3: 231–251.
- 62.
Adélïade M, Sutre G (2004) Parametric analysis and abstraction of genetic regulatory networks.
- 63.
Siebert H, Bockmayr A (2006) Incorporating time delays into the logical analysis of gene regulatory networks. pp. 169–183.
- 64. Ahmad J, Roux O, Bernot G, Comet J, Richard A (2008) Analysing formal models of genetic regulatory networks with delays. Int J Bioinformatics Res Appl 4: 240–262.
- 65. Alur R, Dill DL (1994) A theory of timed automata. Theoretical Computer Science 126: 183–235.
- 66.
Henzinger T (1996) The theory of hybrid automata. pp. 278–292. URL citeseer.ist.psu.edu/henzinger96theory.html.