Integrated Computational Model of Intracellular Signaling and microRNA Regulation Predicts the Network Balances and Timing Constraints Critical to the Hepatic Stellate Cell Activation Process

Activation and deactivation of hepatic stellate cells (HSCs) is an important mechanism contributing to both healthy liver function and development of liver diseases, which relies on the interplay between numerous signaling pathways. There is accumulating evidence for the regulatory role of microRNAs that are downstream from these pathways in HSC activation. However, the relative contribution of these pathways and interacting microRNA regulators to the activation process is unknown. We pursued a computational modeling approach to explore the timing and regulatory balances that are critical to HSC activation and quiescence. We developed an integrated model incorporating three signaling pathways with crosstalk (NF-κB, STAT3 and TGF-β) and two microRNAs (miR-146a, miR-21) that are differentially regulated by these pathways. Simulations demonstrated that TGF-β-mediated regulation of microRNAs is critical to drive the HSC phenotypic switch from quiescence (miR-146ahigh miR-21low) to an activated state (miR-146alow miR-21high). We found that the relative timing between peak NF-κB and STAT3 activation plays a key role driving the initial dynamics of miR-146a. We observed re-quiescence from the activated OPEN ACCESS Processes 2014, 2 774 HSC state upon termination of cytokine stimuli. Our integrated model of signaling and microRNA regulation provides a new computational platform for investigating the mechanisms driving HSC molecular state phenotypes in normal and pathological liver physiology.


Introduction
Hepatic stellate cells (HSC), although important for many aspects of liver function in healthy tissue, are primary drivers of liver diseases such as fibrosis and cirrhosis [1,2].In healthy livers, the main role of HSCs appears to be storage and transport of retinoids as well as modulating the innate immune response [3][4][5].Activation of HSCs occurs in response to numerous stimuli, including induction of activating factors or loss of repressive signals due to genetic or post-transcriptional regulation.During liver regeneration induced by chemical injury (e.g., CCl4 exposure) or partial hepatectomy, HSCs activate and secrete several pro-regenerative proteins, including hepatocyte growth factor (HGF) and epidermal growth factor (EGF) [6][7][8].During diseases associated with chronic inflammation, HSCs alter their gene and protein expression profiles, change their morphology, and deposit fibrous extracellular matrix (ECM), causing scarring and leading to fibrosis and cirrhosis [1].
The process of HSC activation during chronic inflammation is governed by autocrine and paracrine signaling factors from other non-parenchymal and parenchymal cells and is likely highly dependent on the local microenvironment [2,9].It is known that the initial changes in HSCs in response to injury could be a result of paracrine stimulation from other non-parenchymal cells.Several transcription factors play a role in HSC activation.To initiate HSC activation, inflammatory molecules such as IL-1, IL-6, and TGF-β bind to receptors on the stellate cell membrane [10].The presence of constitutive NF-κB activating inflammatory pathways has been reported in HSCs [11].A recent study found IL-1 driven NF-κB inducing multiple MMPs provoking HSC activation [12].IL-6 stimulation of STAT3 increases collagen production in HSCs leading to acceleration of fibrosis.STAT3 directly up-regulates TGF-β expression in HSCs [10].TGF-β initiates a signaling cascade that results in the phosphorylation of SMAD2 (pSMAD2) and further upregulates SMAD7 expression.SMAD7 acts as an anti-fibrogenic factor through an auto-inhibitory feedback loop by inhibiting TGF-β mediated phosphorylation of SMAD2.This feedback loop leads to a restriction of the intensity and duration of TGF-β signaling [13].pSMAD2 forms a complex with SMAD4 and the complex translocates to the nucleus, after which they interact with other factors to regulate gene expression of collagen and other fibrogenic factors [10,14].Interaction between SMAD, NF-κB and STAT3 pathways providing positive and negative feedback modulates HSC activation [4,[10][11][12]14,15].Increased TGF-β expression, lead to increased ECM production by stellate cells and altered ECM composition, with higher fraction of fibrous collagens and lower fractions of basement membrane collagens Activation of HSCs also results in increased levels of growth factors (PDGF, EGF, FGF-2, and others), which induce HSC proliferation, and increased levels of chemotactic proteins (PDGF, MCP-1, CXCR3), which attract HSCs to an activation site [16].
Increased proliferation, chemotaxis, and production of activating signals by active HSCs act as positive feedbacks to amplify HSC activation once it has begun to occur in a subset of cells [2,8].
MicroRNAs are small, non-coding RNAs approximately 21 nucleotides in length which are involved in the regulation of multiple cellular processes [17].Recently, the role of microRNAs governing the activation and function of HSCs has begun to be appreciated.Knockdown of miR-27a and miR-27b in activated HSCs was shown to allow reversion to a quiescent phenotype [18].miR-15b and miR-16 are known to be downregulated during the activation of HSCs [19].Recent studies have revealed the role of miR-29b in HSC activation [20].Another microRNA, miR-146a, has emerged as a potential modulator of HSC activation.miR-146a expression is shown to be downregulated during rat HSC activation [21].A recent in vitro study reported miR-146a expression to be decreased by almost 17-fold in activated HSCs [22].In contrast, when miR-146a was overexpressed in activated HSCs in vitro, the cells appear to take on an anti-fibrotic phenotype, expressing low mRNA levels of IL-6, high mRNA levels of basement collagen (Col1a1), and lower levels of NF-κB binding activity (a mediator of inflammatory response).miR-21 is one of the well-studied microRNAs that modulates cell cycle progression and proliferation during liver regeneration.The role of miR-21 in hepatic stellate cell activation has also been studied [23].Higher expression of miR-21 and the existence of an auto-feedback loop connecting miR-21 to TGF-β pathway activation in activated HSCs have been reported [24].In contrast to miR-146a, miR-21 has been associated with increased cardiac fibrosis in vivo and increased ECM deposition in vitro [25].In HSCs activated in vitro, miR-21 expression increased by over five-fold indicating that miR-21 may play a role modulating fibrosis in HSCs [22].Several other papers have reported such a differential activation of microRNAs miR-21 and miR-146a in the context of liver and other diseases [26][27][28].Alterations in the relative balance between this pair of microRNAs, here termed "microRNA switch", can therefore act a marker for HSC activation state.Activated HSCs have a relatively high expression of miR-21 to miR-146a, whereas quiescent HSCs have a low miR-21 to miR-146a ratio of expression.In summary, the opposing expression and regulation of these microRNAs and their interaction with the cytokine induced signaling pathways are associated with HSC activation.Therefore, we used these two microRNAs as representative markers of HSC phenotype in our model and formulated our analysis based on their differential regulation.
Although the microRNAs miR-146a and miR-21 appear to be useful as biomarkers for HSC activation state, their role in HSC activation remains unclear.The dynamics of miR-21 and miR-146a during HSC activation have not been fully characterized.Several feed forward and feedback interactions between these microRNAs and signaling pathways have been demonstrated in multiple studies [23,25,29,30].The known interactions between signaling pathways and these microRNAs lead to multiple possibilities for how these microRNA biomarkers likely vary during the early stages of HSC activation.We pursued a computational model based approach to explore these possibilities.Based on an array of in vitro and in vivo literature, we developed a novel computational model of HSC activation through IL-1β, IL-6, and TGF-β signaling with feedback from miR-21 and miR-146a.We utilized this model to identify the network balances and timing constraints that may be critical to altering the dynamic switch between miR-21 and miR-146a during HSC activation.
Activation of HSCs is an inherently dynamic multi-scale process due to cell-cell interactions and the tissue microenvironment regulating paracrine signaling and activating intracellular pathways.For the present purpose, we approximated the multi-scale nature of the regulation by considering microRNA based molecular markers as indicators of the HSC cellular state at physiological scales.This approximation constrains our model formulation to focus on key pathway components regulating the switch between microRNAs.By creating an integrated model based on observed molecular and phenotypic states, we are able to study how changes at the molecular level likely translate to changes at the cellular level.

Methods
The integrated signaling pathways and microRNA regulation model were developed as a Petri net system [31].Petri nets are weighted directed graphs that make use of graph theoretical methods, which can efficiently simulate biological signal transduction, metabolic and gene networks.In Petri nets, molecular components are represented as place nodes and biochemical reactions are considered transition nodes.To better represent the complexity of our biological system, we modeled the network as a stochastic Petri net, where stochasticity was introduced as an exponentially distributed firing rate for each of the components (i.e., the likelihood of the waiting times of the nodes to activate had an exponential distribution).The Petri net model was implemented using the software Snoopy (version 1.13, Brandenburg University of Technology, Cottbus, Germany, 2014), which provided a graphical platform for model construction, computational modeling and simulation of the network [32].Transition nodes were modeled to operate by either mass action or Michaelis-Menten kinetics.The kinetic parameters were manually tuned through an iterative process.Initial value ranges were selected from several sources-TGF-β-receptor binding values were estimated from [33], activation of STAT3 and inhibition by SOCS3 from [34], and downstream SMAD/SMURF signaling in the TGF-β pathway from [35].A total of 250 Monte Carlo simulations were run using the Gillespie algorithm for each analysis scenario.In the case of transient signaling, cytokine inputs were modeled as pulse functions using pre-scheduled transitions.For persistent signaling, the initial step input signals were held throughout the simulation.

microRNAs as Markers of Hepatic Stellate Cell Molecular Phenotypes
Activated stellate cells play a major role in the development of liver fibrosis by contributing to the pro-fibrogenic pathways and other signaling pathways leading to ECM production.It has been shown that autoregulatory feedback loops between these pathways and microRNAs influence the fate of the activated stellate cells.Altered microRNA expression levels have been previously reported in activated HSCs [22,23].In particular, miR-21 has been shown to be upregulated in both activated HSCs and regenerating livers, and is known to be over expressed in certain cancer types [22,36].In contrast, miR-146a was shown to be significantly down regulated during HSC activation (~17-fold decrease in expression levels) [22].We analyzed the intensity values obtained from the microarray data for quiescent and activated hepatic stellate cells to generate relative microRNA levels [22].The expression level of these microRNAs varied in opposite fashion during HSC activation (Figure 1).Based on these findings, we postulated that anti-correlation between these two microRNAs and their relative expression levels may be playing a role in regulating the central signaling pathways in activated HSCs and may be used as representative markers of the HSC activation state.Our approach is the first step towards addressing the testable hypothesis that the microRNA behavior based kinetics may contribute to the cellular phenotypic changes leading to HSC activation.

Figure 1.
Changes in microRNA expression levels between quiescent and activated hepatic stellate cells.Data are derived from published microarray data on microRNA expression in primary hepatic stellate cells collected from Wistar rats (Gene Expression Omnibus accession GSE19462) [22].

Integrated Model of Signaling and microRNA Regulation
An overview of the modeled signaling pathways and components is shown in Figure 2. Activation or inhibition arrows ending on species nodes represent influences on production or degradation, respectively.The arrows ending on other arrows represent modification of the corresponding reaction/interaction.For example, PEL1/TRAF6 levels positively influence the transformation of NF-κB inactive into NF-κB active form.NF-κB active species lead to higher production of miR-146a and miR-21.Similarly, miR-146a levels negatively influence the activation of PEL1/TRAF6 complex, and also reduce levels of SMAD4.In order to explore the dynamics of these microRNAs leading to an activated HSC phenotype, we developed a stochastic Petri net model incorporating several well-characterized pathways that have been shown to play a role in the activation of HSCs and stimulation of liver fibrogenesis.
A Petri net implementation of the integrated network is shown in Figure 3.Our model includes IL-1/NF-κB, IL-6/STAT3, and TGF-β/SMAD signaling pathways interacting with miR-146a and miR-21.We explored the contributions of individual pathways to the differential regulation of microRNAs to identify necessary and sufficient regulatory relationships to account for the microRNA switch during HSC activation.We started with the IL-1/NF-κB pathway alone and continued to sequentially expand our model to include IL-6/STAT3 and TGF-β/SMAD pathways.

Model of the IL-1/NF-κB Pathway
Constitutive activity of the transcription factor nuclear factor-κB (NF-κB) has been strongly implicated in the promotion of liver fibrogenesis and HSC activation [37,38].In our model, the proinflammatory cytokine IL-1 activates the NF-κB pathway (Figure A1).The production of IL-1 and IL-1R leads to the formation of an IL-1/IL-1R complex in a concentration-dependent manner.This complex can either dissociate into its individual components or proceed to the activation of the next component in the pathway via signal transduction [39].During persistent signal transduction IL-1/IL-1R activates the PELI1/TRAF6 (Pellino 1/tumor necrosis factor receptor-associated factor 6) complex.NF-κB is represented in the model as distributed across two likely states: an inactive species that is activated by PELI1/TRAF6 complex, and an active species that is capable of upregulating its target genes.The balance between the production and degradation maintains the steady state levels of the inactive species.In the presence of active PELI1/TRAF6 complex, the transition from inactive to active form of NF-κB is modeled using mass action kinetics.The activation of NF-κB is modeled as a reversible reaction approximately accounting for shuttling between cytoplasm and nucleus that is not explicitly represented in the model.In biological systems, the PELI1/TRAF6 complex is the result of interactions between several other proteins, such as IL-1 receptor-associated kinases (IRAK1/IRAK4), and causes activation of NF-κB by mediating its release from inhibitory protein IκB [39].Our model approximates these intermediary reactions as a single reaction with apparent kinetic parameter values yielding an NF-κB activation profile similar to experimental observations.Once activated, NF-κB increases expression of miR-21 and miR-146a [29,30].In the model, microRNA production by activated NF-κB is governed by mass action kinetics.Previous experimental studies have shown that miR-21 targets Peli1, thereby creating a potential negative feedback targeting NF-κB signaling [29].miR-146a has been shown to modulate the IL-1 signaling pathway by targeting IRAK1 and TRAF6 to induce their down regulation [29,30].Based on these observations, miR-21 and miR-146a was modeled as inhibiting the PELI1/TRAF6 complex.This results in a negative feedback loop between the microRNAs and NF-κB activation.
We tested whether activation of NF-κB pathway is sufficient to lead to the experimentally observed switch between miR-21 and miR-146a levels.Simulations showed that both miR-21 and miR-146a levels closely followed the profile of active NF-κB (Figure 4A).Both microRNAs exhibited an initial increase in their levels, followed by a drop-off to steady state levels higher than baseline.Despite the presence of inhibitory feedback, there was no phenotypic switch from a miR-21 low miR-146a high state to a miR-21 high miR-146a low state.Based on the network structure in which NF-κB serves as an activator of both microRNAs, the switch from a miR-21 low miR-146a high state to a miR-21 high miR-146a low state can occur only with differential rates of activation of the two microRNAs.The model was expanded to include additional signaling pathways in order to accomplish such differential microRNA regulation.
It is known that NF-κB acts primarily as a mediator of immune and inflammatory responses in liver cells [37].However, NF-κB in has been shown to function as an anti-apoptotic factor and promote the expression of numerous growth factors and cytokines, including IL-6, in HSCs [2,40].IL-6 plays a key role in activating the transcription factor STAT3, which has recently been found to be a necessary factor in liver fibrosis [10,40].In addition, STAT3 has been shown to induce miR-21 expression in a cancer context to inhibit tumor suppressor pathways [41].We tested whether the addition of the IL-6/STAT3 pathway to our model was sufficient to yield a microRNA switch.[42].IL-6 binds to its receptor, IL-6R, to form the IL-6/IL-6R complex (Figure A2).Binding of ligand-receptor complex IL-6/IL-6R and its dissociation are modeled using mass action kinetics.The IL-6/IL-6R complex leads to activation of STAT3.Our model approximates an intermediate step whereby IL-6/IL-6R activates the protein Janus kinase (JAK), causing STAT3 phosphorylation [42].Once activated, phosphorylated STAT3 (pSTAT3) is able to upregulate the expression of its target genes.Suppressor of cytokine signaling 3 (SOCS3), activated by pSTAT3, targets the signaling intermediary gp130 protein to inhibit STAT3 activation [42].The kinetics of this negative feedback were modeled using Michaelis-Menten formalism with SOCS3 acting as a competitive inhibitor that reduces STAT3 activation.STAT3 activation in HSCs increases miR-21 expression significantly [43].Our model approximates the upregulation of miR-21 by pSTAT3 using mass action kinetics.
Simultaneous activation of NF-κB and STAT3 pathways led to higher levels of miR-21 and a modest increase in steady state miR-146a levels (Figure 4B).However, miR-146a levels remained high, contrary to experimental observations [8], indicating that the activation of these pathways is insufficient to induce the microRNA switch.These results indicated the requirement for negative regulation of one or both of the microRNAs to account for the switch during HSC activation.
Previous experimental studies have shown that STAT3 promotes production of TGF-β, a key factor in driving HSC-mediated fibrogenesis [10].Based on these results, the STAT3 pathway was modified as leading to production of TGF-β.Additionally, crosstalk between TGF-β pathway and microRNAs in liver has been well studied [44,45].Based on these findings we expanded the model to include a TGF-β-mediated pathway leading to activation of genes driving fibrogenesis and likely regulate miR-21 and miR-146a expression.

Model of the TGF-β/SMAD Pathway
Active TGF-β binds to TGF-β receptor I (TβRI), forming TGF-β:TβRI complex (Figure A3).The ligand receptor complex is then able to dissociate, or promote the activation of SMAD2.This activation step accounts for the ligand-receptor complex-mediated phosphorylation of SMAD2 [33,46].Inactive SMAD2 is present at a steady state level balanced by its constitutive production and degradation.The active, phosphorylated SMAD2 species (pSMAD2) is able to bind to another protein, SMAD4, with a 2:1 stoichiometry of SMAD2 vs. SMAD4.This complex, termed pSMAD2/SMAD4 in the model, leads to increase in expression of SMURF2 and SMAD7, which act as negative regulators of TGF-β signaling.The translocation of pSMAD2/SMAD4 into the nucleus has been approximated in the activation rates specified in the model.SMAD7 exerts its influence by blocking SMAD2 phosphorylation and is also able to form a complex with SMURF2 to promote degradation of TβRI [47,48].miR-146a is involved in the inhibition of TGF-β-induced HSC proliferation and activation by decreasing SMAD4 protein expression [45].miR-21 is known to directly target and inhibit SMAD7, thereby creating as a positive feedback to the pro-fibrogenic TGF-β pathway [44,45].
Simulations of the expanded model did not impact the dynamic profiles of microRNAs due to the lack of regulation by the TGF-β pathway.This led us to hypothesize that the shift in balance of the microRNAs is driven by the pro-fibrotic target genes of the TGF-β pathway.Previous studies identified the activation of fibrotic genes by the TGF-β pathway [10].Experimental studies in tissues other than liver have demonstrated an opposing influence of the TGF-β pathway on miR-21 and miR-146a expression [45,49].However, downstream effectors of the TGF-β/SMAD pathway on these microRNAs in HSCs are presently unknown.Based on this phenomenological relationship, we introduced an empirical pro-fibrotic gene (FG) as a downstream target of the TGF-β/SMAD pathway.In our model, activation of FG modulates microRNA expression by upregulating miR-21 and downregulating miR-146a via mass action kinetics.

Crosstalk between Signaling Pathways and Effect on microRNA Levels
Simulations of the integrated network with crosstalk between STAT3 and TGF-β pathway as well as the microRNAs were able to reproduce the expected phenotypic switch of microRNA levels (Figure 5).Our results were consistent with the experimentally observed dynamics of microRNA levels.miR-21 levels showed robust increase upon cytokine stimuli.Following the peak expression, miR-21 levels approached a steady state and were maintained well above baseline expression.However, the initial upregulation of miR-146a due to NF-κB activation was followed by a decline due to the negative regulation mediated by increase in FG levels.Over the course of the simulation, miR-146a expression increased steadily to below baseline levels, closely mimicking the in vitro experimental data on HSC activation [22].These results emphasized the role of crosstalk between the TGF-β pathway and microRNAs to lead to a pro-fibrogenic activated HSC phenotype.
Our results highlight the relative dynamic contribution of exogenous signaling molecules such as cytokines IL-1 and IL-6 produced by other liver cells, as well as endogenous molecules such as TGF-β produced by HSCs, to the microRNA switch underlying HSC activation.

Effects of Timing Differences in NF-κB and STAT3 Activation
We examined the effect of differential timing in activation of NF-κB and STAT3 pathways on the dynamic regulation of microRNAs.Previous experimental results indicated that NF-κB activation occurs earlier, with activity spiking rapidly before decreasing to a lower steady state, whereas peak STAT3 activation occurs later followed by a steady decline [50].However, we considered an alternative scenario in which peak NF-κB activation occurred after that of STAT3, which is likely to happen in the case of altered cell-cell signaling interactions, to evaluate the effects on downstream microRNA expression.Our simulation results revealed that earlier NF-κB peak activity than that of STAT3 resulted in an initial transient increase of both miR-146a and miR-21 expression along with a delayed microRNA switch (Figure 6A).Upon STAT3 activation, miR-21 expression levels further increased and miR-146a levels decreased before both reached a steady state.The alternative case of later NF-κB peak activity led to an initial monotonic decrease in miR-146a levels to levels below the baseline expression, and a persistent increase in miR-21 levels (Figure 6B).Both microRNAs showed a slight increase in expression levels upon NF-κB activation, eventually reaching the same steady state levels as the first scenario.These results suggest that the timing of cytokine stimuli is responsible for controlling the initial dynamics of microRNA expression, but does not significantly alter the overall expression profile or steady state levels.Experimental validation of initial microRNA dynamics could provide insight into the physiological relevance of differential cytokine signaling.

Re-Quiescence from the Activated HSC State
Finally, we explored the dynamics of microRNAs as the system returned from an activated to a quiescent phenotype.The simulation was allowed to reach an activated state phenotype followed by the termination of cytokine stimuli.During re-quiescence, both miR-146a and miR-21 levels decreased initially due to the rapid decline of NF-κB and STAT3 activity (Figure 7).Unlike miR-21, which steadily decreased to baseline levels as STAT3 and FG activity declines, miR-146a recovered from the drop in expression and eventually returned to baseline levels.Although the miR-146a expression level patterns were similar to those previously observed in vitro by Maubach et al., the return of miR-146a to baseline levels is insufficient to completely restore a quiescent HSC phenotype [22].This indicates the involvement of additional factors participating in the complex biological process of re-quiescence.Therefore, accounting for additional microRNAs and other pathways regulating the gene expression might be necessary to capture the multi-scale regulation of HSC phenotypes leading to re-quiescence.

Conclusions
Activation of HSCs, which relies on a variety of interconnected signaling pathways, plays a critical role in the development of liver fibrosis.There is increasing evidence for the role of microRNAs in regulating the transformation of quiescent into activated HSCs [23].Although many of the fibrogenic pathways have been previously studied, further investigation is required to better understand the role of crosstalk between pathways and their regulation mediated by the microRNAs.To address this gap in knowledge, we developed a Petri net system-based model to study regulatory interactions between the microRNAs miR-21 and miR-146a that may occur through as yet unknown direct or indirect mechanisms [22,36].
We performed model simulations to test whether known interactions between signaling and microRNAs can explain the differential regulation of miR-21 and miR-146a.Our results indicate that NF-κB and STAT3 signaling pathways are not sufficient to account for the anti-parallel regulation of these microRNAs.Activation of the TGF-β/SMAD pathway by STAT3 also appears to be insufficient to mimic the measured phenotype.Therefore, we introduced a phenomenological relationship between downstream fibrotic genes (FG) and microRNAs into our model as an additional component that upregulated miR-21 expression and downregulated miR-146a expression.The introduction of TGF-β/SMAD pathway-mediated differential microRNA regulation was necessary to induce a switch between microRNA levels.Upon termination of the cytokine stimuli, the system returned to its quiescent state.This ability of HSCs to switch between activated and quiescent states is critical for maintaining proper liver function, and dysregulation of this process can play a role in the development of liver fibrosis.
Timing of the activation of cytokines plays a significant role in the activation of inflammatory markers during liver regeneration.It is known that timing of cytokine signaling in HSCs is dependent on the type of damage.NF-κB and STAT3 pathways are likely to be activated simultaneously in the case of damage due to infection or alcohol intake.In contrast, in the case of drug-induced liver injury and liver regeneration, IL-1 stimulus is known to precede the IL-6-dependent activation of STAT3 pathway [51].Cytokine signaling in HSCs is the result of paracrine factors from other cells, and the timing of cell communication can determine HSC response to injury.In order to determine how differential activation of cytokine stimuli influences the dynamics of microRNA expression, we analyzed two possible scenarios.These are primarily applicable to the whole liver tissue, however testing these possibilities for HSCs using our model can lead to predictions reconciling the potential difference in temporal relationships between the activation of NF-κB and STAT3 pathways in in vivo vs. in vitro conditions.Our results show that the relative timing between cytokine stimuli in HSCs alters the initial dynamics of microRNA expression.Earlier activation of NF-κB resulted in an initial transient increase in expression of both microRNAs, whereas later activation of NF-κB versus STAT3 resulted in an initial monotonic decrease in miR-146a while miR-21 gradually increased and reach steady state levels.
There is a paucity of evidence in the literature regarding the kinetics of HSC activation.We have used empirical kinetic approximations in order to bridge this gap in knowledge.For parameters, we relied on a range of values obtained from previous modeling studies.However, most of our kinetic parameters were derived through an iterative process.This approach provided us with responses that qualitatively matched the expected kinetics.However, a formal global sensitivity analysis will enable us to fine-tune the model parameters to develop more accurate HSC specific predictions [52].We believe that such a model based on the HSC specific biology can lead to testable hypothesis on HSC activation that can be validated by experiments.For example, the use of exogenous cytokines for in vitro studies has been previously reported for HSCs, and could be easily adapted to investigate the possible role of cell-cell communication in microRNA expression dynamics using HSC cultures [53].Briefly, the cells would be pre-treated with IL-1 or IL-6, followed by addition of the alternate cytokine at a later time point, and microRNA levels would be quantified over this time series.Upregulation of microRNAs can be accomplished by addition of a precursor for the microRNA of interest, and could help determine, for example, the effect of miR-146a overexpression on re-quiescence [43].In addition to overexpression, knockout studies can be utilized both in vitro and in vivo to investigate the role of key components in HSC activation.For example, transfection of cultured HSCs with anti-miR sequences or siRNA against a gene can be used to investigate the relative importance of specific microRNAs and genes in HSC activation, as well as the consequences of their removal.This could be used to determine whether or not, for instance, TGF-β or miR-21 is necessary for HSC activation.These studies can be expanded to include transgenic knockout models, in the case of gene targets, or with the introduction of anti-miR pharmacological compounds, e.g., Vivo-Morpholino [43].In order to study the effects of these changes in activated HSCs, it is important to ensure that the cells achieve total differentiation from the quiescent state.To this end, it has been reported that HSCs require a stiff culture environment in order to undergo activation [54].However, it is not entirely clear how the biomechanical aspects specifically regulate signaling pathways and microRNAs considered in the present model.
Our model accounts for a representative marker of HSC activation representing the switch in miR-21 and miR-146a levels, and can be expanded to consider additional mechanisms of microRNA regulation of target genes, e.g., transcriptional degradation versus translational repression.Further avenues for model expansion include incorporation of co-regulated microRNAs as well as signaling and metabolic pathways involved in HSC activation and re-quiescence [2,22].Our model provides a new computational platform for studying HSC activation and the mechanisms where microRNAs act as fine-tuning markers modulating the system behavior, thereby contributing to potential therapeutic approaches.

Figure 2 .
Figure 2. Schematic view of the integrated network model of signaling pathways and microRNA regulatory interactions.Black lines represent activating interactions, while red lines represent inhibitory interactions.Hypothesized interactions are represented by dashed lines.

Figure 3 .
Figure 3. Graphical representation of the Petri net model.Light gray nodes represent "logical" elements that allow for crosstalk; dark gray nodes represent input stimuli.Edges with arrows correspond to mass flux.Edges with round ends ("read edges") correspond to influences without affecting the levels of the source nodes.Dashed edges ("modifier edges") correspond to inhibitory influences.

Figure 4 .
Figure 4. Dynamics of microRNA and transcription factor expression levels during HSC activation.(A) Expression levels considering the IL-1/NF-κB signaling pathway alone; (B) Expression levels considering IL-1/NF-κB and IL-6/STAT3 signaling pathways regulating the microRNAs.Solid arrowheads mark the start of input cytokine stimuli.

Figure 5 .
Figure 5. Dynamics of microRNA and transcription factor expression levels of the integrated model with all three signaling pathways.

Figure 6 .
Figure 6.Dynamics of microRNA expression levels in response to differential timing of cytokine stimuli.(A) Expression levels with synchronized cytokine production; (B) Expression levels with staggered cytokine production (IL-6 produced before IL-1).Solid arrowheads denote the introduction of input cytokine signals.

Figure 7 .
Figure 7. Dynamics of microRNA and transcription factor expression levels during re-quiescence.Solid arrowheads denote the introduction of input cytokine signals.Hollow arrowheads denote the cessation of cytokine signals.