Predictive health intelligence: Potential, limitations and sense making

1.   Introduction

Hepatocellular carcinoma (HCC), is the most frequently occurring form of primary liver cancer and the main reason for tumor-related deaths globally ^[1]. Despite advances in technology for screening, diagnosis, and treatment, there is still a significant rise in the rate of incidence and mortality worldwide. Rapidly growing HCC cells that provide a site for early vascular invasion, tend to be unaffected by currently available chemotherapeutic options ^[2]. This leads to an increased incidence of recurrence and metastasis post-surgery and results in poor prognoses ^[3]. Over the previous decade, several studies have significantly increased our understanding of the pathogenesis of HCC, but the exact molecular mechanism involved in the metastasis of HCC is still ambiguous.

There exist a plethora of signaling pathways which regulate the system to maintain a homeostatic condition. These regulations are highly significant in preventing the cells from metastasis. Once these regulations are disturbed, there is a chance that the system might move towards a pathogenic state. Thus, there is an imminent need to explore the critical molecular regulatory mechanisms in the key signaling pathways involved in HCC metastasis which will help in the identification of new therapeutic targets. TGF-β and PI3K/AKT pathways are the key signaling pathways that have been identified during the progression of HCC. The misregulation of the TGFβ signaling pathway may result in tumor development. According to the recent studies, matrix metalloproteinases (MMPs) like MMP-2, MMP-1, MMP-9, and MMP-8 along with transforming growth factor beta1 (TGF-β1) actively contribute to HCC progression. They play a significant role as the major key players in metastatic and invasive tumor cells ^[4]. MMPs or matrix metalloproteinases are a group of proteolytic enzymes that are responsible for the collective degradation of the components of proteins found in the basement membranes and extracellular matrix (ECM). Several studies indicate that the over-expression of MMP- 2, 1, 7, 9 and 14 lead to invasion as well as metastasis of HCC through the induction of epithelial to mesenchymal transition (EMT) ^[5,6,7]. It has been implied that enhanced expression of MMP-9 is highly anticipated in invasive HCC cells as compared to non-invasive cells. Moreover, during HCC, the upregulation of MMP-9 has been associated with metastasis, tumor recurrence and invasion ^[5]. Previous reports indicate that it may act as a tumor promoter or as an anti-cancer agent in highly rare situations ^[6]. This was observed from studies performed on MMP-9 knockdown mice models where the tumor incidence was significantly decreased ^[7]. In addition, overexpression of MMP-9 also correlates to tumor progression and reduced survival rate of HCC affected individuals ^[11,12]. Nevertheless, the importance of MMP-9 in HCC progression still needs to be clarified.

TGF-β1 or transforming growth factor beta1 is a widely known cytokine that mediates cellular processes like apoptosis, cell differentiation, and homeostasis. It acts as a driver of tumor progression by inducing EMT and is involved in the activation of several MMPs, including MMP-9 ^[13,14]. During breast cancer, TGF-β1 promotes metastasis by upregulating the expression of MMP-9 ^[10]. Likewise, MMP-9, 2, and 14, can induce the biological activity of TGF-β1 via cleavage of latent TGF-β-binding protein-1 ^[11]. Furthermore, MMP-9 is involved in HCC cell invasion and EMT induced by autophagy through the stimulation of TGF-β/Smad3 signaling ^[12]. A constant relationship between TGF-β1 and MMP-9 has been established through numerous in vitro studies. Moreover, it has been implied that genetically targeted TGF-β1 could alter expression intensities of MMP-9 and MMP-2 in HCC infected cells ^[13]. Despite the evidence, the molecular mechanisms of positive feedback loop between TGF-β1 and MMP-9 on HCC still remain ambiguous.

The latent TGF-β1 ligand is activated by either a surface-bound or soluble MMP-9 (Figure 1). The activated TGF-β1 ligand adheres to the cell surface type Ⅱ receptor (TβRII) which is activated by autophosphorylation and recruits the type Ⅰ receptor (TβRI). TβRI is then activated by phosphorylation. This heterotetrameric complex then stimulates the association of the intracellular signal mediators. TGF-β1 stimulates phosphatidylinositol 3-kinase (PI3K) via phosphorylating Akt, a downstream effector. The PI3K/Akt pathway exists as a non-Smad pathway that is known to contribute to EMT prompted by TGF-β ^[14]. In case of head and neck cancer, stimulation of PI3K/Akt intensifies the expression of MMP-9, leading to degradation of E-cadherin, hence promoting invasion and passage of cells ^[15]. MMP-9 production is enhanced by stimulated transcriptional activity (mediated by Akt/PKB) of nuclear factor-kappa B.

It is well established that HCC progression is caused by the progression of liver fibrosis i.e., deposition of ECM proteins. Henceforth, hepatic stellate cells (HSCs) along with EMT reportedly play a part in HCC progression and metastasis ^[16]. Also, TGF-β1 stimulates fibrogenesis through induction of the HSCs, thus stimulating HCC progression ^[17]. Also, MMPs are responsible for regulating the remodeling of ECM, hence they are critical in HCC progression. It has been previously testified that TGF-β1, MMP-8, 2 and 9 can trigger EMT in HCC cells and lead to pathogenesis ^[16]. These studies reveal that MMP-9 and TGF-β1 can induce progression of HCC by triggering the activation of HSCs and/or EMT.

1.1.   Computational modeling

Modulation of gene expression determines the overall dynamics of each cell and this is a complex process ^[18]. Computational and mathematical approaches in systems biology expedite to observe dynamics of a biological system where entities like RNA, DNA, proteins, enzymes and other biological molecules are involved ^[19]. The qualitative modeling framework is amongst the most recognized methods to study the dynamics of gene expression ^[25,26] by constructing BRNs (biological regulatory networks). Systems in biology are usually modeled using partial or ordinary differential equations, which depend on time derivatives of expression levels, along with kinetic rates of entities. The system is represented in a detailed manner by using quantitative models, which requires data which is either specific or too difficult to find. Hence, qualitative models (BRN, petri net) are preferred to understand the complex dynamics of the biological systems. Construction of the BRN requires qualitative thresholds and logical parameters, which can be easily adjusted according to the system dynamics and biological observations.

1.2.   Our contribution

In this study, we examined the prospective interplay between TGF-β1 and MMP-9 in HCC progression and homeostasis. A BRN was constructed comprising significant entities encompassing TGF-β1 signaling to determine the inhibition and activation relationships. Different parameters of the model are generated from prior experimental knowledge by using the model checking technique. The parameters are then rendered into a qualitative model to analyze the significant system behaviors, such as, cycles, stable states, and paths that lead to recovery and pathogenesis. Therefore, extensive analysis of the resulting model is performed to find these paths. This was followed by a prediction of cycles that symbolize the condition of homeostasis in the overall system, along with deadlock states that reveal the over-expression of the entities, from where the system can't move towards recovery. The results of our study highlight the reciprocal and persistent activation of TGF-β1 and MMP-9 as salient contributors in pathogenesis. Finally, we convert the model into a continuous petri net (PN) to study the evolution of proteins and their relative expression levels based on time. Together, our findings suggest that the TGF-β1 and the MMP-9 feedback loop is critical in tumor progression and targeting these mediators may represent a novel therapeutic strategy to make progress in the treatment of cancer.

2.   Materials and method

A biological system is comprised of numerous elements that are tightly regulated through signaling networks in which the expression level of an entity may either trigger or inhibit the synthesis rate of additional entities or itself. Generally, conventional models (such as ODE and PDE models) are used to represent the entities in a biological system through differential equations that implicate time derivatives for the estimation of various quantities like levels of concentrations, the rate of reactions and temperatures. These models are continuous in nature and rely on specific kinetic values, which are typically not available for all the entities in the network and have been assumed. Hence, the modeling approaches based on graph theory make use of linear methods to observe the behaviors exhibited by biological systems. The methods employed in the current study are divided into consequent sub-sections presented in Figure 2.

2.1.   Qualitative modeling framework

In the 70s, a method of qualitative modeling framework was put forward by Rene` Thomas based on Boolean logic to model BRNs ^[26,27]. The framework comprises various qualitative parameters and thresholds with discrete variables, having only two levels of expression (0 or 1) that can be used to analyze the dynamics of a biological regulatory network. This approach makes use of interaction graphs to imply the basic topology of a BRN.

2.2.   Definitions:

Definition 1 (directed graph): "A directed graph is an ordered pair D (V, E), where:

Ⅰ. V is the set of all nodes and

Ⅱ.E⊆V × V is the set of ordered pairs called edges or arcs"

An edge e.g. (a, b) is directed from an entity or node "a" to "b", where "a" is the tail and "b" is the head of that respective edge. In a directed graph, D^-(x) and D⁺(x) denote the set of predecessors and successors of a specific node x$\epsilon$V, respectively".

Definition 2 (biological regulatory network): "A BRN is a labeled directed graph D (V, E), where V is a set of nodes which represents biological entities and E ⊆ V × V is a set of all possible edges, which represent the interaction between entities".

Ⅰ. "Each edge can be labeled with a pair of variables (σ, ψ), where σ represents the qualitative threshold levels and is a positive integer and ψ is "+" or "–" representing the type of interaction, which can either be "activation" or "inhibition", respectively".

Ⅱ. "Each node e.g. "a" has a limit (l_a), in its threshold level, which is equal to its out-degree (the total number of outgoing edges from "a"). This relation can be presented by ∀b$\epsilon$D⁺(a) and σ_ab$\epsilon${1, 2, 3, ….., r_a} where r_a ≤ l_a which means that the threshold levels of entity "a" can be set within a range "1" to "total number of outgoing edges" and because it has only one outgoing edge towards predecessor "b" so the threshold level which can be set for it can only be "1" ".

Ⅲ. "Each entity, e.g. "a", has its abstract expression in the set Z_a = {0, 1, 2, …., r_a}".

Definition 3 (states): "The state of a BRN is a tuple s$\epsilon$M, where M in terms of entity "a" is:

The qualitative states are represented by vector (M_v)_a$\epsilon$V, where v denotes the level of expression of an entity like "a". According to this definition, M is the Cartesian product of the sets of abstract expressions of all entities. A qualitative state represents a configuration of all the elements of a BRN at any instant of time. The number of activators of a variable at a given level of expression is represented by its set of resources (see the definition of resources given below)".

Definition 4 (resources). "The set of resources Rv_a of a variable a$\epsilon$V at a level v is defined as $R{v_a} = \{ b \epsilon {D^ - }\left(a \right)|{v_b} \ge {\sigma _{ba}}\; {\rm{and}}\; {\psi _{ba}} = +)\; {\rm{or}}\; ({v_{b < }}{\sigma _{ba}}\; {\rm{and}}\; {\psi _{ba}} = -)\}$. The dynamic behaviors of BRN depend on logical parameters. The set of these logical parameters is defined as $K\left(D \right) = \{ {K_a}\left({{R_{va}}} \right) \epsilon {Z_a}\forall a \epsilon V\}$.

The parameter K_a(Rv_a) (at a level v of a) gives the information about the evolution of a. There are three cases: 1) if v_a > K_a(Rv_a) then v_a increases by one unit 2) if v_a > K_a(Rv_a) then v_a decreases by one unit and 3) if c = K_a(Rv_a) then v_a cannot evolve from its current level. It is convenient to describe the evolution from one level to another by an evolution operator "$\Rsh$" [68], which is defined in terms of entity "a" as follows:

Where v_a and K_a(R_va)∈$\mathbb{Z}$_≥0

Definition 5 (state graph): "Let D be the BRN and v_a represents the expression level of an entity e.g. "a" in a state s$\epsilon$M. Then the state graph of the BRN will be the directed graph G = (S, T), where S is set of states and T⊆S×S represents a relation between states, called the transition relation, such that s→s'$\epsilon$T if and only if:

Ⅰ Ǝ a unique a$\epsilon$V such that s_va≠s_va' and s_va' = s_va $\Rsh$ K_a(R_va) and

Ⅱ ∀b$\epsilon$V\{a}s_vb' = s_vb"

2.3.   Inference of parameters through model checking

Computation tree logic (CTL) formulae are formed using various state and path quantifiers. State quantifiers represent the state which holds a property, wherein the path quantifier defines a path from the present state, where a property holds ^[22].

Following are the descriptions of the path and state quantifiers:

∀: This quantifier is read as "For all paths". This state quantifier enforces that all paths initiating from the present state must be holding the given property.

$\exists$: This is read as "There exists a path". It specifies that a minimum of one path starting from the present state must be holding the given characteristic.

$\square$: This state quantifier is read as "globally", which specifies that in a path initiating from a current state, all the states, must hold the given property, including the current state.

O: This quantifier is read as "Next". It enforces that the state which is the immediate successor of the present state must be holding the given property.

$\underline{\Diamond}$: This is a state quantifier which specifies that in a path initiating from a current state, that state or at least one of the future states must hold the given property. This is known as the "Future Quantifier".

SMBioNet ^[23] is software for estimating parameters by employing the approach of Bernot et al. ^[24]. The BRN based on Rene Thomas formalism, where CTL formulae and the set of interactions in the form of activation and inhibition are provided to this software. Furthermore, it computes the probable sets of parameters which are further confirmed by a model checking software, NuSMV ^[25]. NuSMV employs the Snoussi and observability constraints to select the parameter combinations according to the CTL formula. The combinations are further validated to select the models which best suits the behaviors exhibited by our BRN.

Definition 6 (betweenness centrality): "For a state graph R = (S, T) of an interaction graph G = (V; E), let x, y and z be the distinct qualitative states in $\mathcal{R}$, and let σ_{x, y} be the total number of trajectories from state x to state y, and let σ_{y, x} be the total number of trajectories from qualitative state y to x, passing through a state z. Let Ox represents the set of all ordered pairs, (y, x) such that x, y, and z are all distinct. Then, the betweenness centrality of the qualitative state z can be computed from the following equation:"

2.4.   Construction and analysis of BRN:

A BRN is constructed based on the Boolean logic formalism established by Rene Thomas to model the interplay between TGF-β1 and MMP-9 in hepatocellular carcinoma using the GINSIM version 2.4 ^[26] (available at http://ginsim.org/downloads) and GENOTECH tools ^[27] (https://github.com/DrJamilAhmad/GENOTECH/blob/master/GenoTechE.jar). The BRN comprises a set of nodes that represent the biological entities, connected by edges that are directed from one node/vertex to another ^[24]. The edges represent positive and negative interactions amongst the entities. Positive interactions (illustrated by + 1) depict activation, while negative interactions (illustrated by -1) depict inhibition. Experimental knowledge is encrypted in a model checker (SMbioNET) to assign logical parameters to each entity. Furthermore, the BRN is simulated as a state graph to identify the important trajectories, cycles (homeostasis) and steady states (disease progression). The network analysis is performed in Cytoscape, to identify various paths leading to a pathogenic state or the state of homeostasis.

2.5.   Conversion of BRN to petri net (PN) model

A PN is a directed bipartite graph where places (illustrated by circles) and transitions (illustrated by squares) denote entities of a pathway and the processes amongst them, respectively. Moreover, the transition and places are connected by directed arcs to allow the flow of tokens in the pathway. By firing transitions, the source can affect the number of tokens consigned to the target, called the token-count. Therefore, this enables the signals to circulate through the directed protein interactions in a cellular pathway.

For the analysis of signaling networks, PN has emerged as a reliable tool. This approach allows the user to vary inputs, to create a flow of signal through the network based on connectivity of the network, eliminating the need for kinetic parameters. So we employed this approach by converting our BRN to a continuous PN to study the system behavior in a continuous and timed manner.

Formally, continuous PNs formalism is defined as follows:

Definition 7 (continuous petri nets): "A continuous petri net is a quintuple CPN = (P, T, f, v, m0), where:

P, T is finite, non-empty, disjoint sets. P is the set of continuous places. T is the set of continuous transitions."

"f: ((P × T) ∪ (T × P)) → R0 + defines the set of directed arcs, weighted by non-negative integer values"

"v: T→H is a function, which assigns a firing rate function ht to each transition t, whereby H: = ∪t ∈T {ht | ht : R|•t|→R +} is the set of all firing rate functions, v(t) = ht for all transitions t ∈ T, and |•t| represents the cardinality of the preplaces of transition t (reaction's precursor)"

"m0: P→R0 + gives the initial marking."

We use the tool GINsim ^[26] to export a standard PN file which can be imported in SNOOPY ^[28] for further analysis.

3.   Results

3.1.   Pathway abstraction and BRN model construction based on prior knowledge

A detailed literature review was carried out to identify pathways and entities that are critical in inducing the feedback loop between TGF-β1 and MMP-9. The extended pathway shown in Figure 1 was further reduced to a BRN (Figure 3) such that the impact of interactions amongst all the entities was significantly highlighted. This particular method of abstraction is explained in detail by Naldi et al 2009 and Saadatpour et al ^[35,36]. The reduced BRN comprises many well-characterized regulatory motifs that give rise to a particular functionality of the system. The BRN is constructed with five nodes representing the entities: TGF-β1, MMP-9, NFκB, TIMP-1 and HSC. The integers, + 1 and -1 are labeled with directed edges to represent activation (illustrated with a straight line) and inhibition (illustrated with a red dashed line) induced by the cellular and non-cellular components. The BRN highlights a positive feedback loop between TGF-β1 and MMP-9 through reciprocal activation. This loop exerts certain changes that direct the system to move away from its state of equilibrium, causing instability and eventually leading towards the state of disease i.e., HCC. A positive feedback loop consists of an even set of positive elements ^[31]. Whereas, a negative feedback loop is composed of an odd set of negative interactions, for instance, the interaction found between timp-1 and HSC tend to move towards either homeostasis or an interrupted behavior. Analysis of these motifs may provide significant insights into the probable dynamics of a system. The presence of TIMP-1 negatively regulates the expression of MMP-9, while its production is positively regulated by HSCs. Nonetheless, the dynamic behavior of a complex system, that comprises both negative and positive feedback loops, can be rendered with an appropriate set of parameters.

3.2.   Parameter inference and selection via model checking

The constructed BRN has 5 entities: TGF-β1, MMP-9, NFκB, TIMP-1 and HSC (Figure 3). Each of these five entities comprises several discrete parameters that show the level of each attribute implicated in the model.

SMBioNet (selection of models of biological networks) ^[23] has been applied to compute parameters according to known qualitative observations in the form of computational tree logic (CTL) formulae (Table 1).

It takes BRN as inputs and selects the models which verify these formulae. The formula was thoughtfully derived from biological observations indicating an increase in expression from an initial state with the help of state and path quantifiers. The CTL formula is based on the fact that the enhanced expression of TGF-β1, MMP-9, and TIMP-1 is often associated with invasiveness and metastasis of ^[13,32,35]. TGF-β1 and MMP-9 are known to promote the progression of HCC via stimulating HSCs and NF-κB is a key transcription factor in the regulation of MMP-9 expression ^[35]. Hence, the CTL formula encodes these observations that during HCC pathogenesis all the entities are upregulated due to dysregulation of the tightly regulated elements in the network (i.e., TGF-β1 = 1, MMP-9 = 1, TIMP-1 = 1). The analysis from SMBioNet resulted in four sets of parameters in the form of models which satisfied the CTL formula's properties (Table 2).

Amongst these, the fourth set was selected. The fourth model was selected based on prior knowledge from the literature that in the presence of its activators, MMP-9 is upregulated as the HSCs frequently secrete the protein along with TGF-β and NFκB that too are responsible for the activation of MMP-9. In qualitative modeling, these selected parameters aid in observing the system dynamics in the form of cycles, stable states (SS) and paths. The selected parameters are presented in Table 3.

3.3.   Analysis of the state graph: Insights into pathogenic and recovery states

A state graph of the BRN was simulated from the selected parameters (Table 3) using the GINsim software ^[26]. The logical set of parameters of the four models, which satisfied the CTL properties, were analyzed using the state graph (as shown in Figure 4, starting from M1 to M4). In this representation, a logical variable associated with each pertinent element of the system where "0" represents the normal to down-regulation of an entity, whereas "1" represents up-regulation of an entity. Following model selection, the simulated state graph consists of 32 nodes and 80 edges which were rendered using Cytoscape ^[36] (available at http://www.cytoscape.org/) (Figure 4) and the states were sorted on the basis of betweenness centrality. The state graph represents all possible transitions in the form of states that denote relative expressions of all entities at a specific time. The order of the qualitative states is TGF-β1, NFκB, MMP-9, HSC, and TIMP-1. The deadlock state (1, 1, 1, 1, 1) shows up-regulation of the all the entities. The model also generates several cycles, computed by GENOTECH software, which reveals normal system behavior categorized by down-regulation of TGF-β1 and MMP-9, with a normal expression of other entities. The state of the system, in a state graph, is denoted by a vector that includes the level of expression of all the entities. The analysis revealed that the system can either lead to a disease state or a recovery state. The recovery state is denoted as (TGF-β1 = 0, NFκB = 0, MMP9 = 0, HSC = 0, TIMP1 = 0), where all the entities exhibit low expression. On the contrary, the pathogenic state is indicated by high expressions of TGF-β1 and MMP-9 that induce the overexpression of the other entities leading to a deadlock. Under normal circumstances, biological systems maintain homeostasis where the system remains in a cycle. Hence, a desired qualitative model should demonstrate a closed path showing homeostasis along with trajectories that lead to pathogenesis. It is also revealed through our analysis that there is one significant state which results in either a deadlock state or can lead to recovery. This state (TGF-β1 = 1, NFκB = 0, MMP9 = 1, HSC = 0, TIMP1 = 1) as shown in Figure 4 highlights that TGF-β1 and MMP-9 exhibit high expression levels, but they seem to be regulated by the expression TIMP-1, which is an inhibitor of MMP-9. As soon as TIMP-1 is down-regulated, both the entities overexpress each other resulting in the overexpression of NF-kappaB. The expression of MMP-9 is further upregulated by the high expression level of NF-kappaB. Also, NF-κB is known to augment the invasive ability of tumor cells ^[37]. This indicates that the system moves toward HCC pathogenesis. On the other hand, if TGF-β1 or MMP-9's expression is down-regulated the system may move towards the state of recovery, indicating the fact that inactivation of MMP-9 is crucial to recovering from the disease. Since MMPs play a key role in stimulating tumor metastasis and angiogenesis ^[38], it is confirmed by our results that agents required for down-regulating the expression levels may prove useful in handling HCC treatment outcome.

3.4.   Selection and analysis of the most probable homeostatic cycles

It is significant to identify the most apparent biological cycle since homeostasis is represented by cycles. Our model generated 30 cycles of varying lengths and one of them was selected on the basis of betweenness centrality rendered in the Cytoscape tool ^[36]. Betweenness centrality is defined as "the number of shortest paths that go through a node among all shortest paths between all possible pairs of nodes" ^[39]. All states were sorted based on their measure of betweenness centrality as presented in Figure 5, where nodes with a larger diameter and darker color signify states with high betweenness centrality. Amongst all the cycles, the cycle with highest betweenness centrality was as follows: (0, 0, 1, 0, 1) → (1, 0, 1, 0, 1) → (1, 0, 0, 0, 1) → (1, 0, 0, 0, 0) → (1, 0, 0, 1, 0) → (1, 1, 0, 1, 0) → (1, 1, 0, 1, 1) → (0, 1, 0, 1, 1) → (0, 1, 1, 1, 1) → (0, 0, 1, 1, 1) → (0, 0, 1, 0, 1).

The cycle initiates at the point where the expression of MMP-9 is tightly regulated by its inhibitor (TIMP-1). In the following state, the expression profile of TGF-β1 shifts from down-regulation to up-regulation. As soon as this state is met, MMP-9 is downregulated in the next state of the cycle, preventing activation of the loop and signaling the system to remain in a regulated cycle. The cycle also reveals that down-regulation of either TGF-β1 or MMP-9 is essential for homeostasis. Any deviation from this cycle would head towards a pathogenic stable state i.e., (1, 1, 1, 1, 1).

3.5.   Conversion of BRN to continuous petri nets

It is very challenging to find defined kinetic parameters for each reaction; hence we depend upon the PN model so that we are not restricted to certain values and expressions. The signaling pathways in biology have evolved to an extent where the network connectivity has an alleviating effect on the systems' entities and their associated interactions. Hence, the PN analysis results in an easy interpretation of the simulations run to predict the biological systems behavior under variable stimuli, be it internal or external.

The BRN modeled in GINSIM was exported into a standard petri nets (PN) format where the places and transitions were discrete in nature. Since the biological entities tend to be continuous in nature, the standard petri net was then converted into a continuous petri net using SNOOPY tool (version 2) (Figure 6). Modelling with PN involves the use of simulations to study the changes in relative expression levels of biological entities i.e., genes/proteins over time. There are two places for each network entity in the converted PN, i.e., active and inactive state. The activated states are represented by a prime place, and the inactivated state is shown by a complementary state. Furthermore, the related transitions are labeled as "p" and "n" to represent the activating and inactivating signals, respectively. The transitions are created by GINsim to model the parameters generated in SMBioNet. This aids in analyzing the BRN dynamics over time. Over the course of disease or homeostasis, the biological entities are activated in a variable manner. Hence, the need to include inactivated places for proteins is significant so to understand the inhibition signals that are received intracellularly. The resulting PN in our study is simulated on the basis of mass action kinetics for continuous transitions, where the presence of tokens is represented by the markings in the places.

3.6.   Analysis of the continuous evolution of key entities

Simulation makes obvious the time-based changes of the concentration of molecules and thus enables to graphically represent their paths and check the refined properties of the system ^[40]. Figure 7 shows the simulation results in the case of recovery where MMP-9 and TGF-β decline towards the state where the system maintains homeostasis leading to the downregulation of all the other entities (i.e., HSCs, NFκB, and TIMP-1). It reveals that the downregulation of MMP-9 results in lowering the expression of TGF-β1. Furthermore, NFκB also declines towards homeostasis mediating the downregulation of TIMP-1 to recover from the pathogenic state.

In contrast, Figure 8 depicts a clear picture of how all the entities are increasing towards the disease state and then becomes constant as it reaches a pathogenic stable state. It is observed that the overexpression of NFκB results in an increase in MMP-9 that may enhance invasion and progression. Moreover, high expression of MMP-9 and TGF-β1 can induce progression of HCC via induction of HSCs. Also, the relative expression level of MMP-9 is higher than that of TIMP-1 in the simulation result. It can be clearly seen in Figure 8 how the overexpression of MMP-9 leads to the activation and upregulation of TGF-β confirming our hypothesis where the reciprocal activation of both entities is considered as the potential molecular mechanism in accelerating the progression of HCC.

The proposed premise has been suggested in a recent study where the idea has been put forward for further analysis ^[13]. The role of TGF-β1 has been highlighted in HCC by means of two mechanisms: First, as a paracrine or autocrine growth factor through an intrinsic activity and, second, by stimulating microenvironment changes by an extrinsic activity ^[41]. In various experimental models, inhibitors of TGF-β1 have revealed to hinder the progression and growth of HCC by modulating EMT. TGF-β inhibitors arise as a novel therapeutic approach, which targets the cancer growth and its microenvironment. As TGF-β is present at almost every stage of tumorigenesis, blocking it may have a more powerful effect than agents that block specific processes during the pathogenesis of HCC, such as neo-angiogenesis ^[42]. TGF-β inhibition is directed against the most central pathological processes which lead to chronic liver disease and consequently to HCC. Over-expression of MMP-9 in tumor progression via cell invasion has been correlated with invasion and metastasis. Inhibitors for suppressing the expression of MMP-9 have been revealed through recent studies ^[43,44,45]. Taken together, therapies targeting both mediators may lead to medical breakthroughs in the treatment of cancer.

4.   Discussion

A consistent relationship between TGF-β1 and MMP-9 has been well-established in studies carried out on various fibroblasts and cancer cells ^[13]. The high expression of MMP-9 has been associated to invasion and metastasis in multiple tumors because of its capacity to remodel the tissues via the ECM, along with its ability to cause angiogenesis and degradation of the basement membrane ^[46,47]. In HCC, increased mRNA and protein levels of MMP-9 have been associated with highly invasive cancer cells ^[5]. Moreover, TGF-β1 has proven to be an independent biomarker to indicate survival in HCC ^[48]. By directly activating various signaling pathways and transcription factors, TGF-β1 can upregulate the expression of MMP-9, leading to an increased proliferation in various tumor cells ^[16]. Similarly, the link between the two key players was demonstrated in some studies, that MMP-9 increases the TGF-β1 bioactivity by proteolytically cleaving the latent peptide of TGF-β1. Finally, this relationship between MMP-9 and TGF-β1 forms a potential regulatory feedback loop correlated with fibrogenesis and HCC via the activation of hepatic stellate cells (HSCs) ^[16]. Therefore, it is significant to identify the mechanism of the positive interplay in the TGFβ1-associated BRN that plays a key role in triggering the signaling pathway that leads to pathogenesis.

Formal modeling approaches are extensively applicable because of their computational capacity of testing exhaustively. These approaches are being used successfully over a few years for the purpose of modeling and validation of intricate systems in biology ^[49]. Boolean logic formalism is a well-established method for modeling a qualitative BRN that decrypts the dynamics of the qualitative model in the form of a directed graph. A qualitative state is represented by a node, whereas the edges demonstrate the interactions between the nodes ^[21,50,51]. The qualitative model of the TGFβ1-associated BRN was simulated as a state graph, which resulted in cycles and a deadlock stable state. The most apparent biological cycle was selected on the basis of betweenness centrality that shows the expression profile TGF-β1 and MMP-9 is tightly regulated by TIMP-1, which controls the up-regulation and down-regulation of both the entities throughout the cycle, so only one of the key players exhibits expression in a given state. On the other hand, the down-regulation of TIMP-1 results in the overexpression of the TGF-β1 and MMP-9 via reciprocal activation. This tends to lead the system towards a deadlock stable state (1, 1, 1, 1, 1).

To maintain homeostasis, the expression profiles of genes go through numerous levels (such as low or high) during the regulatory mechanisms. The regulation and activation of the TGF-β1 and MMPs is found to be perturbed in tumor tissues and has been associated with increased progression of the tumor ^[52,53]. It has been suggested that when the MMPs activate latent TGF-β1, they can mediate metastasis and invasion in the stroma by exhibiting intense tumor-promoting effects ^[54]. The deadlock state generated by our model exhibits this behavior, where MMP-9 and TGF-β1 are found to be overexpressed. In recent studies, activated TGF-β1 has been demonstrated to mediate the critical balance of the remodeling of extracellular matrix by modulating the expression level of TIMPs and MMPs which is demonstrated in the qualitative cycle ^[54]. The simulation from the continuous petri net (PN) model was also consistent with the finding that the reciprocal activation of both the major entities (i.e., TGF-β1 and MMP-9) is considered as the potential molecular mechanism in accelerating the progression of HCC. It also suggests that TGF-β inhibitors may be used as a novel therapeutic approach, which targets the cancer growth and its microenvironment. As TGF-β1 is present at almost every stage of tumorigenesis, blocking it may have a more powerful effect than agents that block specific processes during the pathogenesis of HCC, such as neo-angiogenesis ^[42]. Inhibitors for suppressing the expression of MMP-9 have been revealed through recent studies ^[43,44,45]. Taken together, further wet-lab based experiments are required to explore the roles of TGF-β1, MMP-9, and TIMP-1 to target the key players in the tumor microenvironment.

5.   Conclusion

The impact of reciprocal activation of TGF-β1 and MMP-9 in HCC has been elucidated in prior experiments. However, the precise molecular mechanism behind this interplay has not been established yet. In this study, we employed a robust computational technique to analyze the loop between TGF-β and MMP-9, which triggers invasion and metastasis in HCC. We found that both entities must be expressed in a regulated manner and that the inactivation of MMP-9 may aid in restoring the system in a recovery state and preserve homeostasis. Moreover, when MMP-9 forms a positive regulatory loop with TGF-β1, the system approaches a point from where recovery is beyond the bounds of possibility i.e., the deadlock stable state. These findings propose that this loop serves as a significant mediator of EMT resulting in HCC growth and development. Taken together, these findings may aid in the development of treatment strategies that are designed to target both TGF-β and MMP-9. Our constructed model could aid in further understanding of the dynamics of the disease by adding and exploring other entities (genes/proteins/cells) involved in HCC. Furthermore, the model can be extended with necessary perturbations and interventions to observe the effect of the entities on expression patterns.

Acknowledgements

We thank our colleagues from Atta-ur-Rahman School of Applied Biosciences (ASAB) and Research Center for Modeling and Simulation (RCMS) who provided thoughtful insight and great expertise that assisted the research.

Conflict of interest

All authors declare no conflicts of interest in this paper.