RESEARCH ARTICLE

# A rule-based multiscale model of hepatic stellate cell plasticity: Critical role of the inactivation loop in fibrosis progression

Matthieu Bouguéon¹,², Vincent Legagneux², Octave Hazard³,⁴, Jérémy Bomo¹,², Anne Siegel¹, Jérôme Feret⁴,⁵*, Nathalie Théret¹,²*

1 Univ Rennes, Inria, CNRS, IRISA, UMR 6074, Rennes, France, 2 Univ Rennes, Inserm, EHESP, Irset, UMR S1085, Rennes, France, 3 École Polytechnique, Palaiseau, France, 4 DI-ENS (Inria, ÉNS, CNRS, PSL University), École normale supérieure, Paris, France, 5 Team Antique, Inria, Paris, France

* jerome.feret@ens.fr (JF); nathalie.theret@inserm.fr (NT)

<figure><img src="image_1.png" /></figure>

## Abstract

Hepatic stellate cells (HSC) are the source of extracellular matrix (ECM) whose overproduction leads to fibrosis, a condition that impairs liver functions in chronic liver diseases. Understanding the dynamics of HSCs will provide insights needed to develop new therapeutic approaches. Few models of hepatic fibrosis have been proposed, and none of them include the heterogeneity of HSC phenotypes recently highlighted by single-cell RNA sequencing analyses. Here, we developed rule-based models to study HSC dynamics during fibrosis progression and reversion. We used the Kappa graph rewriting language, for which we used tokens and counters to overcome temporal explosion. HSCs are modeled as agents that present seven physiological cellular states and that interact with (TGFβ1) molecules which regulate HSC activation and the secretion of type I collagen, the main component of the ECM. Simulation studies revealed the critical role of the HSC inactivation process during fibrosis progression and reversion. While inactivation allows elimination of activated HSCs during reversion steps, reactivation loops of inactivated HSCs (iHSCs) are required to sustain fibrosis. Furthermore, we demonstrated the model's sensitivity to (TGFβ1) parameters, suggesting its adaptability to a variety of pathophysiological conditions for which levels of (TGFβ1) production associated with the inflammatory response differ. Using new experimental data from a mouse model of CCl4-induced liver fibrosis, we validated the predicted ECM dynamics. Our model also predicts the accumulation of iHSCs during chronic liver disease. By analyzing RNA sequencing data from patients with non-alcoholic steatohepatitis (NASH) associated with liver fibrosis, we confirmed this accumulation, identifying iHSCs as novel markers of fibrosis progression. Overall, our study provides the first model of HSC dynamics in chronic liver disease that can be used to explore the regulatory role of iHSCs in liver homeostasis. Moreover, our model can also be generalized to fibroblasts during repair and fibrosis in other tissues.

<figure><img src="image_2.png" /></figure>

### OPEN ACCESS

**Citation:** Bouguéon M, Legagneux V, Hazard O, Bomo J, Siegel A, Feret J, et al. (2024) A rule-based multiscale model of hepatic stellate cell plasticity: Critical role of the inactivation loop in fibrosis progression. PLoS Comput Biol 20(7): e1011858. https://doi.org/10.1371/journal.pcbi.1011858

**Editor:** Philip K Maini, Oxford, UNITED KINGDOM OF GREAT BRITAIN AND NORTHERN IRELAND

**Received:** January 24, 2024

**Accepted:** July 5, 2024

**Published:** July 29, 2024

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1011858

**Copyright:** © 2024 Bouguéon 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.

**Data Availability Statement:** The 3 model family files are available in the GitHub repository: https://github.com/MBouguéon/HSC_model_Kappa.
Using the KAPPA user interface downloadable from
https://kappalanguage.org/download, the
simulations can be performed by uploading the file
of the model family reactMFB-with-inactivation.ka.
We used the version 4.2.12 of Kappa and 4.13.1 of
OCaml (required to run Kappa). The python scripts
use to launch parallelized simulation are available in
the GitHub repository: https://github.com/
MBougueon/simulation_pipeine.

**Funding:** This work was supported by the Institut National de la Sante et de la Recherche Medicale (Inserm)(NT, VL,MB), the Institut national de recherche en informatique et en automatique (INRIA)(JF), the Centre National de la Recherche Scientifique (AS) and the University of Rennes 1 (JB). 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.

Author summary

Chronic liver diseases (CLDs) are associated with the development of fibrosis which is characterized by an abnormal deposition of extracellular matrix (ECM) leading to severe liver dysfunction. Hepatic stellate cells (HSCs) are key players in liver fibrosis driving ECM remodeling. However numerous biological processes are involved including HSC activation, proliferation, differentiation and inactivation and novel computational modeling is necessary to integrate such complex dynamics. Here, we used the Kappa graph rewriting language to develop the first rule-based model describing the HSCs dynamics during liver fibrosis and its reversion. Simulation analyses enabled us to demonstrate the critical role of the HSC inactivation loop in the development of liver fibrosis, and to identify inactivated HSCs as potential new markers of fibrosis progression.

Introduction

Liver fibrosis is an excessive wound healing response induced by chronic injuries, mainly
caused by viral hepatitis (HCV, HBV), alcohol abuse and non-alcoholic steatohepatitis
(NASH). Fibrosis is characterized by an accumulation of extracellular matrix (ECM) which
increases the stiffness of tissues, leading to severe liver dysfunction. The final stage of fibrosis,
cirrhosis, leads to complications such as ascites, variceal hemorrhage, encephalopathy and
hepatocellular carcinoma, and is associated with a high mortality rate worldwide [1]. Activa-
tion of hepatic stellate cells (HSCs) is the main process underlying hepatic fibrosis [2]. In a
normal liver, HSCs are quiescent, store vitamin A and are located in the Disse space between
hepatocytes and endothelial cells that delineate sinusoids. Upon liver injury, HSCs are acti-
vated and transdifferentiate into ECM-secreting myofibroblasts (MFB) that contribute to tis-
sue repair. In chronic liver disease, the inflammatory signal persists, leading to sustained HSC
activation, ECM accumulation and fibrosis [2-4]. Importantly, reversibility or regression of
liver fibrosis after elimination of the inflammatory agent has been observed in experimental
models and clinical studies [5, 6]. This phenomenon is associated with the elimination of myo-
fibroblasts, involving processes of apoptosis, senescence and inactivation, but the contribution
of each of these mechanisms to the repair/fibrosis balance remains unclear.

Understanding the behavior of HSCs during liver injuries requires multiscale modeling approaches. Such approaches have been widely developed for modeling biological processes [7] and diseases [8, 9]. Because fibrosis is a pathological tissue repair activity that occurs in dif- ferent organs, unified approaches have been proposed by including common components such as the inflammatory response and extracellular matrix remodeling [10]. However, each tissue is characterized by specific cell microenvironments. Different multiscale models have therefore been developed for cardiac [11] and pulmonary [12] fibrosis, and various approaches have been used to model the progression of liver fibrosis. Among them, the agent-based (ABM) multiscale model developed by Dutta-Moscato et al. [13] successfully reproduced the experimental collagen deposition observed in carbon tetrachloride (CCl4)-induced fibrosis in rats. The model integrated parenchymal cells, inflammatory cells, collagen-producing cells and molecular regulatory agents, with rules defining the properties and interactions between all agents. This model was then modified by preventing the migration of collagen-producing cells, thus providing a more accurate dynamic of collagen deposition [14]. The Dutta-Moscato model has also been extended and modified by Wand and Jiang [15] to include information about lipid accumulation induced by CCl4 treatment, making it possible to study the progres- sion of liver fibrosis in presence or absence of steatosis. In addition to these agent-based
models, Friedman and Hao recently published a partial differential equation (PDE) model for
liver fibrosis that includes information on inflammation regulation and ECM remodeling,
enabling exploration of anti-fibrotic drugs [16]. Although instructive, all these models reduce
the dynamics of HSCs to their activation to become ECM-producing myofibroblasts, without
taking into account the plasticity of HSCs recently illustrated by large-scale single-cell analyses.
[17–19]. Among the complex behaviors of HSCs, the authors identified heterogenous cell phe-
notypes such as proliferating cells and ECM-producing cells. In addition the inactivation pro-
cesses remain poorly understood and information is lacking about the fate of inactivated HSCs
(iHSCs) during liver fibrosis and reversion. The present study aims to decipher HSC dynamics
by developing multiscale models which integrate for the first time all the HSC states during the
development of fibrosis and its reversion. The multi-state combinatorial nature of the cells
studied led us to choose a rule-based model (RBM) approach, much more appropriate for
dealing with this type of network. Among these formalisms, we selected the Kappa language
[20–22]. Using the Kappa language, our model describes the interactions between all cell pop-
ulations and molecules such as transforming growth factor β 1 (TGFβ1) which induces HSC
activation and promotes secretion of type I collagen, the major component of ECM [23]. We
validated the predicted ECM dynamics using new experimental data from a mouse model of
CCl4 induced liver fibrosis. Importantly, simulation studies demonstrated the regulatory role
of HSC inactivation loops during liver fibrosis and predicted the accumulation of inactivated
HSCs during chronic liver disease. By analyzing RNA sequencing data from 102 patients with
non-alcoholic steatohepatitis (NASH) associated with liver fibrosis, we confirmed this accu-
mulation, identifying iHSCs as novel markers of fibrosis progression.

**Results**

Three rule-based model families for HSCs dynamics

We integrated information based on literature, to describe HSC dynamics during liver fibrosis and its reversion however the lack of some information led us to propose hypotheses and finally to develop three different families of models. As shown in Fig 1A, the three families share common processes involving the different states of HSCs regulated by TGFβ1. In normal liver, HSCs are maintained in a non-proliferative quiescent state (qHSC) and store vitamin A. Upon liver injury (virus, toxic, etc.), the inflammatory response leads to the production of TGFβ1 which activates qHSCs into an activated state (aHSC). The fully activated HSCs are characterized by a myofibroblast (MFB) phenotype. Upon removal of liver injury, i.e in the absence of TGFβ1 in the models, MFBs are either eliminated through apoptosis and senescence pathways (apop_sene_MFB) or reverted to an inactivated HSC state (iHSC) close to but different from the quiescent state. Upon a new liver injury, i.e in the presence of TGFβ1 in the models, the iHSCs are reactivated (react_HSCs). In the same way that HSCs transform into MFBs, react_HSCs transform into reactivated MFBs (react_MFBs) however information is lacking about the fate of these cells. To overcome this issue, we developed two families of models in which the react_MFBs are either i) completely eliminated by the apoptosis/senescence pathways but not by inactivation, this family of models called *reactMFB-wo-inactivation* (wo for without) is described in Fig 1A₁, or ii) eliminated in the same way as MFBs, i.e. by the apoptosis/senescence pathways and/or the inactivation pathway, this family of models called *reactMFB-with-inactivation* is described in Fig 1A₂).

Similarly to react_MFB, information is lacking about the fate of iHSCs. While the activation
of iHSCs by TGFβ1 has been described using *in vivo* and *in vitro* experiments [24, 25], no data
have been reported about the fate of iHSCs in the absence of TGFβ1. However, Kisseleva *et al.*
[24] showed a decrease in the iHSC population during the reversion process. We therefore
<figure><img src="image_3.png" /><figcaption>**Fig 1.** **Three models to describe HSC dynamics.** A1 Schematic representation of models *reactMFB-wo-inactivation*; A2 Models *reactMFB-with-inactivation* and A3 Models *iHSC-reversion-to-qHSC*. (B) Schematic representation of Kappa-modeled processes in all models. HSCs are agents characterized by seven cell physiological states (qHSC, aHSC, MFB, iHSC, apop_sene_MFB, react-HSC, and react_MFB). Rules describe the processes: (1) qHSCs are activated into aHSCs, (2) aHSCs are transformed into MFBs, (3) MFBs are transformed into apop_sene_MFBs to be eliminated, (4) MFBs are inactivated into iHSCs, (5) iHSCs are reactivated into react_HSCs upon TGFβ1 stimulation, (6) react_HSCs are transformed into react_MFBs. Similarly to MFBs, react_MFBs are eliminated by apoptosis and senescence (3) or by inactivation (4). (7) aHSCs and react_HSCs proliferate, (8) MFBs and react_MFBs proliferate, (9) aHSCs and react_HSCs produce COL1, (10) MFBs and react_MFBs produce COL1, (11) apop_sene_MFBs are degraded, (12) iHSCs are degraded. (13) qHSCs self-renewal. The effect of TGFβ1 is represented by gray arrows, it induces the activation of qHSCs, the reactivation of iHSCs, and the production of COL1 in aHSCs, react_HSCs, MFBs and react_MFBs. The reversion of iHSCs into qHSCs are represented by a dotted arrow. (C) Kappa contact map. HSCs are agents with three sites: *cell_state* (qHSC, aHSC, MFB, iHSC, react_HSC, react_MFB, apop_sene_MFB), *TGFB1_bind*ing (free or bound) and *TGFBR* (membrane, internalized or degraded).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g001

developed a third family of models in which, in the absence of TGFβ1, iHSCs are either elimi-
nated or return to a quiescent state, thereby participating in the repopulation of the injured
liver with new qHSCs. This potential new source of qHSCs needs to be balanced with the self-
renewal of native qHSCs in order to restore HSC homeostasis in the repaired liver. This mod-
els in which iHSCs revert to a quiescent state constitutes the third family called *iHSC-rever-
sion-to-qHSC* which is illustrated in Fig 1A₃.

Finally, we developed three families of models: *reactMFB-wo-inactivation*, *reactMFB-with-inactivation*, *iHSC-reversion-to-qHSC*. The biological processes of all the models have been integrated in *Fig 1B*, where the interactions between TGFβ1 and aHSC, react_HSC, MFB and react_MFB cells leading to the induction of their proliferation and the production of COL1 are indicated.

The three families of models were built using the Kappa rule-based language (see Methods
section) and simulations were performed using the KaSim tool [26] based on the Gillespie sto-
chastic simulation algorithm (SSA) [27]. In these Kappa models, the hepatic stellate cells are
the agents which are characterized by three sites (Fig 1C) which are: i) the *cell_state* of HSCs
(qHSC, aHSC, MFB, iHSC, react_HSC, react_MFB, apop_sene_MFB), ii) the *TGFB1_bind*
that can be free or bound and iii) the state of the TGFβ1 receptor (*TGFBR*) that can be
localized at the membrane, internalized or degraded. In addition to these agents, we intro-
duced *counters* to scale the intermediate steps between cell states, and we used *tokens* to
describe TGFβ1 and COL1 quantities instead of detailing the behavior of each molecule,
thereby highly reducing the computational cost (see *Methods* section).

The *reactMFB-wo-inactivation* models contain 75 rules and 41 parameters. A rule for *react_MFB* inactivation and two parameters controlling the proportion between inactivated and eliminated *react_MFB* have been added to obtain the *reactMFB-with-inactivation* models, containing 76 rules and 43 parameters. An additional rule for iHSC reversion and two other parameters controlling the proportion between reverted and eliminated iHSC have been added to obtain the *iHSC-reversion-to-qHSC* models. These rules are grouped into 7 biological processes: 1) qHSC renewal, 2) binding of TGFβ1 to HSCs, 3) TGFBR dynamics, 4) HSC activation and differentiation, 5) HSC proliferation, 6) COL1 production and elimination and 7) MFB inactivation (detailed in *Methods* section).

Inactivation loops of react_MFBs are essential to maintain COL1 accumulation during chronic liver injury

To identify the families of models that predict HSC behaviors and collagen accumulation that
fits with the experimental data from Kiesseleva et al. [24], we conducted simulation studies
using the protocol described by the authors for a mouse model of liver fibrosis and reversion,
i.e. twice-weekly CCl4 (TGFβ1 in the model) injections for 8 weeks and an overall reversion
period of six months. Observations reported by Kiesseleva et al. [24] included a 1.43-fold
increase in cells expressing α-SMA and a ~ 14-fold increase in the percentage of COL1 depos-
its at 8 weeks and a return to initial values at 6 months after the last injection. In addition, the
authors assessed iHSCs/qHSCs ratios and identified a 50/50 distribution at one month and a
total number of HSCs similar to initial conditions at six months of reversion. As shown in Fig
2, most qHSCs (98%) were immediately activated upon the first stimulation of TGFβ1 across
all models, a small number (2±1%) of qHSCs persisting during the stimulation period. Activa-
tion of qHSCs led to the production of aHSCs that differentiated into MFBs. These cells
express α-SMA protein, and we tracked α-SMA positive cells by summing the occurrences of
aHSCs, MFBs, react_HSCs and react_MFBs. Note that even if the number of qHSCs and
iHSCs is low, they can be activated or reactivated by TGFβ1 during repeated stimulation for
two months. In line with this observation, Rosenthal et al. [28] recently suggested that iHSCs
remain at a low level even during the activation phase. Comparison of the models showed that
the react_MFB inactivation loop greatly affects model behavior. When inactivation of
react_MFBs was not allowed (models *reactMFB-wo-inactivation*) (Fig 2A), we observed a lower
amount of α-SMA expressing cells (0.8-fold±0.02) and of COL1 deposits (9.5-fold±0.16) than
expected at 8 weeks, with amounts even decreasing before stimulation was stopped. Con-
versely, when inactivation of react_MFBs was allowed (models *reactMFB-with-inactivation*)
(Fig 2B, 2C and 2D) we observed an accumulation of α-SMA expressing cells and COL1
deposits that varied according to the percentage of react_MFB inactivation. While 5% allows
the model to fit with experimental observations (red circle in Fig 2), increasing the percentage
to 10 and 50% led to higher amounts of α-SMA expressing cells (1.6±0.05 and 5-fold±0.03,
respectively) and COL1 deposits (16±0.36 and 40-fold±0.16, respectively) at 8 weeks. Only
the model with 5% of react_MFB inactivation led to a 50/50 distribution of the iHSCs/qHSCs
ratio at one month and a return to the initial total number of cells at six months of reversion,
as described in the experimental data [24].

Based on the models *reactMFB-with-inactivation* parameterized with 5% of *react_MFB* inactivation, we next evaluated the impact of a potential reversion of iHSCs to a quiescence
<figure><img src="image_4.png" /><figcaption>**Fig 2.** Time-course analyses of cells and COL1 in models without or with *react_MFB* inactivation. A series of simulations was performed using conditions of stimulation reported in Kisseleva *et al.* [24]. We used 10,000 TGFβ1 molecules per cell and 16 stimulations (twice a week) during 2 months. **A1** Models *reactMFB-wo-inactivation*. **B1** Models *reactMFB-with-inactivation* in which 5% of *react_MFBs* are inactivated. **C1** Models *reactMFB-with-inactivation* in which 10% of *react_MFBs* are inactivated. **D1** Models *reactMFB-with-inactivation* in which 50% of *react_MFBs* are inactivated. Stimulation of cells with TGFβ1 started on day 4 to allow the model to equilibrate. **A2,B2,C2** and **D2** display the variation in the number of cell occurrences for qHSCs, aHSCs, MFBs, iHSCs, *react_HSCs*, *react_MFBs* and α-SMA positive cells. The number of α-SMA cells is the sum of the number of aHSCs, MFBs, *react_HSCs* and *react_MFBs*. Simulations are expressed as the mean of 10 replicates. **A3, B3, C3** and **D3** display the variation of COL1% area. Simulations are expressed as arbitrary unit (a.u) and all 10 replicates are represented. The red circles indicate the biological observations reported in [24]. From left to right (**A2,B2,C2** and **D2**), the first circle is for an 1.43-fold increase in cells expressing α-SMA at 8 weeks, the</figcaption></figure>
second circle is for an equal proportion of iHSCs and qHSCs at one month of recovery and the third circle is for an equal number of cells at 6 months and of cells at T0. The circle for COL1 simulations (A3, B3, C3 and D3) is for a ~ 14-fold increase in the percentage of COL1 deposits at 8 weeks.

https://doi.org/10.1371/journal.pcbi.1011858.g002

state. As described in these new models *iHSC-reversion-to-qHSC*, we observed that regardless of the percentage of reversible cells, the overall behavior of the other model components was not affected (see S1 Fig). This may be due mainly to the slow dynamics of iHSCs elimination and quiescent cell renewal, which limit the impact of reversion. We have discarded this model, which was developed on the hypothesis of iHSC reversion but without any supporting experimental observations.

Taken together, our data demonstrate that inactivation of react_MFBs is required to sustain COL1 accumulation by maintaining the level of α-SMA expressing cells, thus validating the models *reactMFB-with-inactivation*. Simulation analyses allowed us to determine that a 5% ratio between reac_MFB inactivation and apoptosis_senescence pathways reflected the cells and COL1 behaviors seen experimentally. Regarding iHSCs, our observations suggest that complete reversion of iHSCs to a quiescent state is a potential pathway to eliminate iHSC in the absence of TGFβ1. Since the outcome generated by the models *iHSC-reversion-to-qHSC*, in which iHSCs can revert to qHSCs was similar to the one obtained in the models *reactMFB-with-inactivation* (S1 Fig) we chose to retain the latter for the rest of the analyses.

## Collagen I accumulation is sensitive to the dynamics of TGFβ1 stimulation

Our models are calibrated using experimental data from a mouse model of CCl4-induced fibrosis [24], a reproducible experimental model widely used to study fibrosis and repair [29]. Although the protocols vary somewhat in terms of the number, frequency and modalities of CCl4 injections, these experimental models induce similar liver fibrosis. Other experimental models of induction of liver fibrosis in rodents have been developed, including injection of other chemical agents (such as thioacetamide (TAA) and dimethylnitrosamine (DMN)), high-fat diets and surgical interventions (such as bile duct ligation) [30, 31]. Whatever the etiology, induced damages lead to TGFβ1-dependent activation of HSCs, but with variable dynamics. It can happen very quickly in a model induced by CCl4 or slower in a model induced by a fatty diet. To explore the sensitivity of our validated model *reactMFB-with-inactivation*, we varied the concentration of TGFβ1 and the number and periodicity of stimulations. As shown in Fig 3, the model was highly sensitive to the amount of TGFβ1, but this sensitivity varied with the number of stimuli and periodicity, particularly for TGFβ1 values lower than 10,000 molecules per cell. Beyond that, the number of molecules exceeded the number of receptors per cell, and the model became saturated, as shown by the superposition of kinetic curves at 10,000 (green) and 100,000 (yellow) molecules per cell. At lower TGFβ1 levels, collagen accumulation was much more subtle. When stimulations were spaced sufficiently far apart, collagen levels returned to baseline and, remarkably, we observed a much stronger induction upon second stimulation, particularly at doses of 10,000 and 100,000 molecules per cell (Fig 3A, periodicity 30 and 60 days). This induction of collagen 1 was less pronounced at low doses, but increasing the number of stimuli (Fig 3B, periodicity 180 hours and 15 days) showed that collagen continues to accumulate with each new addition of TGFβ1 at the dose of 1000 molecules per cell, whereas it decreased for doses 10,000 and 100,000 molecules per cell. This effect remained visible for a larger number of stimulations (Fig 3C, periodicity 180h and 15 days and Fig 3D, periodicity 90h and 180h).
<figure><img src="image_5.png" /><figcaption>**Fig 3.** Time course analysis of COL1 as a function of the variation in TGFβ1 concentration, number and periodicity of stimuli. Results are expressed as the number of occurrences of COL1 according to the number of TGFβ1 stimuli: (A) 2, (B) 4, (C) 8, (D) 16, (E) 32, the periodicity which varies from 22.5 hours to 60 days and the concentration of TGFβ1: 10 (Red), 100 (Blue), 1,000 (Purple), 10,000 (Green) and 100,000 (Yellow) molecules per cell. The first stimulation occurs at day 4. Simulations are expressed as the mean of 5 iterations. P, Periodicity.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g003

## Model predictions validated by new experimental data in mice

To evaluate the model predictions, we analysed COL1 accumulation in a mouse model of CCl4 induced liver fibrosis. For the quantification of collagen deposits, we used second harmonic generation (SHG) microscopy which we have previously demonstrated to quantify
<figure><img src="image_6.png" /><figcaption>**Fig 4.** Comparison of model predictions and experimental data obtained in a mouse model of CCl4 induced liver fibrosis. A1 Representative SHG microscopy images of collagen in mouse livers after 2, 4, 8 and 10 weeks of CCl4 treatment. A2 For the reversion experiments, treatment was stopped at week 4. A3 Collagen quantification was performed using ImageJ and results are expressed as percentage of SHG signal relative to total area (Mean ± SD). B1, Model simulations were performed using parameters corresponding to the CCl4 injection protocol (three TGFβ1 stimuli in the first week and one stimulus per week for ten weeks). Collagen accumulation was plotted and experimental data (from A1) are indicated with black bars. B2, Model simulations were performed using parameters corresponding to the reversion protocol (three TGFβ1 stimuli in the first week and one stimulus per week for 4 weeks). Collagen accumulation was plotted and experimental data (from A1) are indicated with blue bars.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g004

fibrillar collagen in human [32] and mouse [33] livers. Note that the global evaluation of collagen deposits using Sirius red [24] and SHG microscopy may overestimate the contribution of HSCs since other cells can contribute to liver fibrosis such as portal fibroblasts and bone marrow-derived cells [34]. However, HSCs have recently been shown to be the exclusive source of myofibroblasts in CCl4-treated liver [35] which is why we chose data from CCl4-induced liver fibrosis to calibrate our model. As shown in Fig 4A, we compared the signal quantification data obtained by SHG microscopy in our experimental model with the simulation data. Our experimental data showed maximum accumulation of collagen after 6 weeks and stabilization at a lower level at 8 and 10 weeks after treatment (Fig 4A₁). In tissue samples from mice that have undergone the reversion protocol, the SHG signal shows a decrease in COL1 at 8 weeks after CCl4 removal (Fig 4A₂). For the simulation studies, TGFβ1 parameters were adapted to the experimental protocol, i.e. three TGFβ1 stimuli in the first week and one stimulus per week for ten weeks. CCl4 induces rapid hepatocyte damages leading to inflammation and TGFβ1 production that we modeled using saturating TGFβ1 quantities (10,000 molecules/cell). As shown in Fig 4B₁, the predicted curves for COL1 accumulation reached a maximum value similar to the experimental data, then at 8 and 10 weeks, the stabilization plateau was slightly higher in our simulation. Regarding the reversion protocol, the predicted curves for COL1 are in full agreement with the experimental data (Fig 4B₂).

We also validated the prediction of COL1 dynamics using experimental data from other laboratories, including liver fibrosis models induced by dimethyl_nitrosamine [36] and thioacetamide [37]. Results are presented in (S2A and S2B Fig).
iHSC accumulation is associated with human liver fibrosis

A major observation by Kisseleva et al. [24] is that the fibrotic response is exacerbated in mice that have undergone an initial chronic aggression followed by a reversion period allowing a return to a healthy tissue phenotype. This was associated with the presence of inactivated stellate cells (iHSCs) which display an enhanced response to TGFβ1 when stimulated in cell culture, leading to faster and greater production of COL1 [24, 25, 38]. These iHSCs would therefore still be present in the “repaired” liver after six months of reversion. To take this into account, we calibrated our models *reactMFB-with-inactivation* on Kisseleva’s observations *et al.* [24] after the second aggression by introducing the hypothesis of a slow elimination of iHSCs. However, there is no experimental data describing the fate of these cells during aggression/reversion cycles. Using the relapse protocol described by Kisseleva *et al.* [24], we studied the behavior of iHSCs throughout injury-repair cycles. This protocol includes chronic injuries induced by eight CCl4 injections over one month, followed by a six-month recovery period, and then a new injury similar to the first. We repeated this cycle three times (Fig 5).

Our simulations show that the dynamics of α-SMA expressing cells and amounts of COL1 were consistent with Kisseleva's *et al.* [24] experimental observations, i.e. higher levels of these two markers after the second aggression compared to the first one (Fig 5A and 5B, cycles 1, 2). Importantly, the model allowed us to detail the dynamics of iHSCs that was not experimentally evaluated, and we observed an increased number of iHSCs (Fig 5C, cycles 1, 2). When the aggression/reversion cycles were repeated, we showed an accumulation of iHSCs that was associated with that of α-SMA expressing cells and COL1 (Fig 5A, 5B and 5C, cycles 3, 4). The increased number of iHSCs may result from the balance between the higher number of α-SMA expressing cells leading to the inactivation of more iHSCs and the half-life of iHSCs. On the basis of these observations, the model predicts a progressive accumulation of iHSCs during cycles of aggression/reversion, suggesting a critical role for iHSCs in the dynamics of fibrosis, with accumulated iHSCs constituting a source of COL1 overproducing cells during reactivation.

To validate this prediction, we searched for iHSC accumulation in samples from patients with liver fibrosis. To this end, we analyzed RNAseq data from a large multicenter study comprising 206 histologically characterized liver samples from patients with non-alcoholic fatty liver disease (NAFLD) including 102 patients with non-alcoholic steatohepatitis (NASH) associated with liver fibrosis [39]. Patients were divided in two groups: with low (F0-F1, n = 34) and high (F3-F4, n = 68) fibrosis grades. To identify iHSCs in human liver samples, we used an iHSC gene expression signature reported by Rosenthal *et al.* [28] in a mouse model of non-alcoholic steatohepatitis (NASH). This signature contained 39 genes that were significantly differentially expressed in iHSCs when compared to both qHSCs and aHSCs. In addition, we used 10 iHSC gene markers identified in a mouse model of CCl4-induced liver fibrosis [24] and validated in an *in vitro* reversion model of human activated cells [38]. The list of these genes is detailed in S1 Table.

As shown in Fig 6A, gene set enrichment analyses (GSEA) using the signature identified by Rosenthal et al. [28] showed that iHSC genes were collectively enriched in the high fibrosis group compared to the low fibrosis group (normalized enriched score = 1.38, p-value = 0.097). We performed a new analysis, adding the 10 marker genes identified by Kisseleva *et al.* and El Taghdouini *et al.* [24, 38] to the list of genes identified by Rosenthal *et al.* [28]. As shown in Fig 6B, this additional set of genes increased the statistical power of the GSEA analysis (normalized enriched score = 1.57, p-value = 0.048), suggesting that the expression of this iHSC-specific gene set is increased in patients with high fibrosis. These observations are in line with the increase in the number of iHSCs predicted by our model, during the progression of hepatic fibrosis.
A number of αSMA-cell
occurences (10^3)

<figure><img src="image_7.png" /></figure>

B number of COL1
occurences (a.u)

<figure><img src="image_8.png" /></figure>

C number of iHSC
occurences (10^3)

<figure><img src="image_9.png" /></figure>

Time (days)

<caption>**Fig 5.** iHSCs accumulate during injury-repair cycles. Simulations were performed using repeated cycles comprising a first injury cycle (8 injections twice a week for 1 month), followed by a 6-month recovery period and a second injury cycle as previously described by Kisseleva *et al.* [24]. Simulations are expressed as the mean of 50 replicates for cells (A), and all 50 replicates have been represented for COL1 (B) and iHSCs (C). Overlaps are shown as grey area.</caption>

https://doi.org/10.1371/journal.pcbi.1011858.g005
A

Enrichment Score (ES): 0.4256573

Normalized Enrichment Score (NES) 1.3889463

Nominal p-value: 0.097247705

FWER p-Value: 0.053

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

B

Enrichment Score (ES): 0.46911344

Normalized Enrichment Score (NES) 1.5760338

Nominal p-value: 0.04789272

FWER p-Value: 0.025

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig 6.** iHSC accumulation is associated with liver fibrosis. Gene set enrichment analyses of iHSC gene expression was performed in low fibrosis (F0-F1, n = 34) and high fibrosis (F3-F4, n = 68) samples from patients with non-alcoholic steatohepatitis (NASH) ([39]). A) Analysis performed with the iHSC gene signature identified by Rosenthal *et al.* [28]. B) Analysis performed with the iHSC gene signature identified by Rosenthal [28], supplemented by ten iHSC markers identified in [24, 38]. Results are presented in the form of correlation profiles between gene lists, enrichment scores (ES) which indicate to which phenotype (low or high fibrosis) the gene set is most positively correlated and heat maps visualizing clusters and gene distribution.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g006
Discussion

Regardless of the underlying etiology, all chronic liver diseases are associated with the develop-
ment of hepatic fibrosis, which is a major public health issue. Understanding the dynamics of
hepatic fibrosis and its reversion is essential for the clinical management of patients and the
adaptation of therapies. A pivotal process in the development of hepatic fibrosis is the
TGFβ1-dependent activation of hepatic stellate cells whose plasticity is now widely recognized
[40, 41]. Here, we investigated HSC plasticity by developing a multiscale model using a rule-
based language. We provided the first model to describe the dynamics of HSCs during liver
fibrosis and its reversion, demonstrating the importance of inactivated cells in these processes.

Kappa, a language adapted to multiscale modelling

While ordinary differential equations (ODEs)-based models have long been used to model bio-
logical dynamics, rule-based models (RBM) are much more appropriate for dealing with net-
works with combinatorial multi-state interactions as discussed by Chylek et al. [42]. Indeed,
ODE approaches require the implementation of an equation for each species. Agent-based
models (ABM), while more appropriate for bottom-up approaches than ODEs, require a
description of all possible behaviors of agents. Common ABM approaches are based on the
object-oriented paradigm [43] or process calculi [44].

In the case of our model, the single agent has 4 sites, some of which have 8 or 14 different states. Using an ABM approach would have exploded the number of rules (~ 300) in the model to describe the 39 possible agent states. RBM approaches were therefore the most appropriate to represent the known interactions operating on stellate cells during fibrosis development and reversion. Based on site graph rewriting, BioNetGen [45] and Kappa [21, 22] are the two major languages used for RBM implementation and simulation.

Kappa formalism has been widely used for modeling biological reaction networks such as cell signaling [46], gene regulation [47], epigenetic regulation [48], repressilator system [49] and DNA repair [50]. Using Kappa language, our team previously developed a model to describe the regulatory role of extracellular matrix networks in TGFβ1 activation [51]. The transition from molecular to cell-molecule interaction models has been a new challenge, and we have recently demonstrated the applicability of the Kappa language in multi-scale modeling [52]. However, to overcome the explosion in memory demand due to the declaration of molecules and cells as agents, we have used *tokens* and *counters*. The way *tokens* are stored in memory and the absence of explicit bond with agents considerably reduces the computational cost associated with the existence of different scales in the model. This involves the use of hybrid rules between agent and *token* whose applicability has been demonstrated in our model. *Counters* allowed a more detailed description of the activation and differentiation states of HSCs by creating intermediate steps and cell behaviors depending on their activation/differentiation stage, that better reflects experimental observations. Moreover, *counters* allow to model phenomena that the duration of which is more predictable than chemical interactions. With time-exponentially distributed event, there is only one parameter: standard deviation which is fully characterised by the average time. Erlang distributions provides a mean to narrow the distribution hence leading to events which are more predictable (and which corresponds from a mathematical point of views to the composition of several atomic events following time-exponential distributions).

In short, *counters* and *tokens* are key to multiscale modelling in Kappa. *Tokens* enable a more
compact representation of agent populations, *counters* enable a more compact representation of
site populations, with both facilitating the search and abstraction of potential rule applications.
To our knowledge, beside a model of pulmonary viral infection by William Waites [53], our
study is the first one where *counters* and *tokens* have been used in a Kappa model.
## A unique fibrosis model for different etiologies

While the role of TGFβ1 in HSC activation is well established, our model showed how its concentration, the number of stimuli and their periodicity influence the dynamics of liver fibrosis. These observations underline the adaptability of the model to the different pathophysiological conditions associated with the development of fibrosis. In murine models of liver fibrosis induced by chemical compounds such as CCl4, the toxic effect is massive, with strong hepatocyte necrosis, exacerbated inflammatory response and a saturating release of TGFβ1. In other models, such as diet-induced fibrosis, the immune response is different and fibrosis progresses more slowly. Sensitive to small variations in TGFβ1, our model reproduces the collagen accumulation observed in these different fibrosis mouse models. Similarly, the response to liver injury in human is etiology-dependent and may affect the dynamics of HSC activation. In contrast to chemical models, diet-induced fibrosis models developed to mimic the progression of NAFLD to non-alcoholic steatohepatitis (NASH) are characterized by a lower level of fibrosis. In addition, the dynamics of fibrosis are highly variable depending on the rodent strain and diet composition. Moreover, the optimal model of steatohepatitis leading to fibrosis (choline_deficient, amino_acid defined diet) is poorly reversible [54]. To overcome this heterogeneity, we chose a reversible diet model [55] which we compared with model predictions. We set the parameters of TGFβ1 to low levels (4500 molecules/cell, 6 stimuli and a periodicity of 14 days). As shown in supplementary Fig S2-C, the model captured the increase in COL1 but reversion dynamics were too fast, suggesting that for this condition the model requires additional regulators.

## A model to capture stellate cell heterogeneity and microenvironment remodeling

Recent single-cell studies have demonstrated the great heterogeneity of HSCs and MFBs during the development of fibrosis [17, 28, 56]. These works also show that, in addition to a diversity of phenotypes, the dynamics of HSCs during their activation and differentiation process are far from linear, as shown by pseudotime analyses [57, 58]. Our approach manages to model this heterogeneity by introducing intermediate steps in the activation and differentiation process, and allowing cells in each of these states to proliferate, progress through the activation process and produce COL1. The implementation of these rules combined with the SSA algorithm (Stochastic Simulation Algorithm) means that, at each point in the simulation, there are different phenotypes and a 'unique' history is kept for each agent.

However, modelling HSC dynamics does not capture the full complexity of the microenvironment that governs these dynamics such as the multitude of events regulated by immune cells. Other multiscale models include the role of Kupffer cells in HSC activation via TNF-α [13, 15] or antagonistic regulation by M1 and M2 macrophages [16]. Similarly, it is difficult to include all the components involved in extracellular matrix remodeling. This process is finely regulated by the balance between the production of matrix metalloproteinases (MMPs) and their inhibitors, Tissue Inhibitor of metalloproteinases (TIMPs), both secreted by HSCs, MFBs, but also by immune cells involved in the resolution of fibrosis [59-61]. In their model, Friedman & Hao [16] reduced ECM remodeling to the sole contribution of macrophages and did not take into account the known role of HSCs at all. To overcome this molecular and cellular complexity while keeping remodeling process into our own model, we introduced *tokens* to represent low and high remodeling COL1. In this way we modeled both collagen accumulation during TGFβ1 stimulation and its regression upon TGFβ1 withdrawal.
## Inactivation of reactivated MFBs, a critical factor in the progression of fibrosis

The process of inactivating HSC-derived MFBs has been presented as a key element in the reversion of fibrosis enabling the elimination of ~ 50% of MFBs [24, 25]. The authors have shown that the inactivated cells can be reactivated both *in vivo* and *in vitro* in response to a new stimulus, and that they are much more reactive and fibrogenic than initially activated cells [24, 38]. However, the fate of these reactivated cells remained unknown, leading us to hypothesize a potential re-inactivation of these reactivated cells. Our results showed that the inactivation of these cells was essential to maintain COL1 accumulation in chronic lesions, but that this inactivation had to be limited. Indeed, our models have shown that 5% of reactivated MFBs should be inactivated, experimental observations are no longer respected when using higher percentages of reactivated MFBs. This loop of inactivation and reactivation induces a change in the MFB population. As lesions accumulate, cells inactivated between aggressions are reactivated, resulting in a population of reactivated MFBs that becomes the majority and more fibrogenic cells. In support of these predictions, single-cell sequencing data identified these ECM-producing MFB phenotypes with increase fibrogenic activity [28].

## Reversion of inactivated stellate cells to a quiescent state as an alternative to cellular elimination

Although the existence of iHSCs has been supported by several studies, their behavior is not yet well understood. These cells have a phenotype close to that of quiescent cells, but retain the memory of their previous activation [24, 38]. However, the fate of iHSCs in the absence of reactivation stimuli remains unknown. The development of our models enabled us to observe that the reversion of these iHSCs to a quiescent state, regardless of the percentage of cells affected, did not impact the overall behavior of the other components. This may be due mainly to the time required for the elimination of iHSCs and the renewal of quiescent cells, which is particularly slow in our models. This dynamic limits the impact of reversion. It is also possible that these cells have a faster elimination dynamic, but that this is compensated by a proliferation capacity, thus slowing down their elimination and increasing the impact of potential reversion. As it stands, our models suggest that the complete reversion of iHSCs to a quiescent state is a potential means of eliminating them in the absence of TGFβ1. Experimental follow-up of these cells over reversion times of more than one month would provide more precise information on their fate.

## Inactivated HSCs as new markers of fibrosis progression?

The accumulation of iHSCs during fibrosis is a very important concept from a clinical perspective, as the progression of fibrosis is not at all linear. iHSCs within tissues could act as a kind of memory of an initial injury, making the tissue more sensitive to the next injury. Of course, this memory facilitates rapid repair, but the recurrence of these repair episodes, if too close together, accelerates the progression of fibrosis. Although our model is not fully adapted to follow the dynamics of the mouse NASH model due to the slow release of diet-induced TGFβ1, we showed that the observation of iHSC accumulation is supported by the enrichment of iHSC gene signature in NASH patients with increased fibrosis.

Identifying inactivated cells *in vivo* would provide a genuine marker of patient history. However, this identification remains complex due to the intermediate phenotype of these cells, between quiescent and activated state. There is as yet no established phenotypic signature for inactivated cells, and gene expression signatures vary according to the mouse models used
(NASH, CCl4). These phenotypic signatures depend both on the time of reversion and on the methods used to identify them, depending on whether their expression is compared with that of quiescent cells and/or activated cells. Recent advances in single-cell sequencing will enable us to refine these signatures by characterizing more specific groups of genes.

## Conclusion

The present study provides two major findings. First we have demonstrated the applicability of the Kappa language to multi-scale approaches by using: tokens, to represent the quantities of molecules in a compact way, and counters to abstract the dynamics of cell-specific phenotypes. Second, we developed the first model of HSC dynamics during liver fibrosis and reversion. Thus, our model provides biological predictions for the reactivation loop, the dynamics of inactivated cells and their accumulation during fibrosis progression. These predictions have been validated either experimentally in a mouse model of CCl4-induced fibrosis, or using RNA sequencing data from patients with fibrosis.

## Materials and methods

### Kappa syntax

Kappa is a site graph rewriting language [62–65] that uses a chemistry-inspired syntax to transparently describe the interactions between component occurrences. The syntax used by Kappa is that of site graphs, using rules to describe the behavior of agents over time. An agent describes species (e.g. cells, proteins) and defines the characteristics of these species. Rules define both the interactions between agents and their behavior. Agents are described by one or more sites which can have different states. These states allow agents to establish links between themselves. A rule describes a process that sets a condition (left-hand member) and an event (right-hand member). All the elements at the left of the → are the elements needed for the rule to be applied, the elements at the right are what will be produced. In the following example, two agents A and B interact to form a complex AB.

$$A(x\{p\}, y[.]), B(y[.]) \rightarrow A(x\{p\}, y[1]), B(y[1]) @ 'k' \quad (1)$$

The agents A and B and their respective sites are written as: Agent_name/site_name{state_value}[binding_state]). The agent A has its site x in a state p and its site y free (.), and the second agent B has it site y free. A and B will bind together via their site y. The symbol @ is followed by the expression that defines the rate of each potential application of the rule in the state of the system. That is to say the probability that a potential event effectively occurs within an infinitesimal interval of time.

Simulating models in which agents have to be explicitly bound to other agents can be computationally expensive. Instead, such bounds can be described implicitly. To this aim, the binding state of the binding site in an agent A can be encoded as an internal state which specifies whether the site if free, or bound. In the latter case, the site to which this site is bound to is not specified. Then the amount/quantity of agent B remaining is encoded by the means of a token. Consequently, Tokens constitute a pool of variables enabling the description of abundant molecular species in a continuous manner. Their use allows us to resolve the memory limitation caused by declaring molecular species as an agent. If a token such as ATP is added to the rule 1, the following rule is obtained (rule 2), where the binding of A and B requires ATP. What stands between | and @ defines the quantity of token to be deleted or added (2) when the rule is applied. Because the binding of A and B depends on the amount of ATP, the rate of the
rule has to be modified to include the quantity of ATP, written as $|ATP|$.

$$A(x\{p\}, y[.]), B(y[.]) \rightarrow A(x\{p\}, y[1]), B(y[1]) | - 10'ATP'@ |ATP|*k' \quad (2)$$

Models in which agents have numerous states require many rules to describe the changes between these states. Instead of describing all these states using numerous rules, it may be sufficient to know the number of steps to follow a process, requiring only a few rules, by increasing the step index. With *counters*, we can use discrete intervals to describe a special type of *site* which can be increased or decreased in the rules. In the following example (rule 3), A has a *site* called *n*, which is a *counter* that counts the number of bonds that can be formed by A during its lifetime in the system. This new rule allows A to bind to B as long as the value of *n* for A is greater or equal to one. When a bond is formed, the value of *n* decreases by one.

$$A(x\{p\}, y[.], n\{\ge 1\}), B(y[.]) \rightarrow A(x\{p\}, y[1], n\{-= 1\}), B(y[1]) @ 'k' \quad (3)$$

However, the current syntax of Kappa cannot allow to test whether the value of a *counter* is inferior to a given value (i.e., the "<" symbol is not supported). To overcome this limitation, it is necessary to use two *counters* whose sum remains constant. This syntactic limitation is overcome in the new version of Kappa (which internalises the use of pairs of *counters* to deal with less than inequality tests).

## Modeling HSC dynamics

Using Kappa rules, we described HSC dynamics by including various biological processes such as the activation of quiescent HSCs (qHSC) by TGFβ1 leading to activated HSCs (aHSC), the transdifferentiation of aHSCs into Myofribroblasts (MFB), the elimination of MFBs after removal of TGFβ1 either by apoptosis and senescence pathways (apop_sene_MFBs) or by reverting to an inactivated state (iHSC), the reactivation of iHSCs upon new TGFβ1 stimulation leading to reactivated HSCs (react_HSC) and the transdifferentiation of react_HSCs into reactivated MFBs (react_MFB). We also used rules to describe cell proliferation and production of collagen (except for qHSCs and iHSCs). In these models, HSCs are agents characterized by 3 sites: i) the cell state of HSCs (qHSC, aHSC, MFB, iHSC, react_HSC, react_MFB, apop_sene_MFBs), ii) the TGFB1 binding (free or bound) and iii) the state of the TGFβ1 receptor (TGFBR) that can be localized to the membrane, internalized or degraded. We introduced a pair of *counters* (called *intermediary-step* and *control_counter*) to scale the intermediate steps between cell physiological states, and we used 4 *tokens* to describe TGFβ1 and collagen molecules. The resulting families of models contained 75 to 77 rules and 37 to 41 parameters, depending on families. All rules are available in the GitHub repository: https://github.com/MBougueon/HSC_model_Kappa. We detailed here the representative rules of the biological processes. These rules are given in graphic form for easier reading.

**qHSC renewal.** qHSC renewal is modeled using two rules (Fig 7). These rules take into account the fact that the site *TGFB1 binding* on qHSCs must be *free*, since binding of TGFβ1 induces activation into proliferating HSC.

**TGFβ1 binding to cells.** To describe the interaction of TGFβ1 with HSCs, 7 different rules are used according to the cell state. Fig 8 shows the binding rule of TGFβ1 to qHSCs. Only the site *cell_state* of the agent differs between these rules, except for MFBs and react_MFBs, where the rule rate is decreased. The fact that each agent can bind TGFβ1 introduces competition between agents for the binding of TGFβ1.

**Turnover of TGFBR.** The endocytic pathways regulating *TGFBR* receptor signaling and turnover has been previously described by Di Guglielmo *et al.* [66] and modeled by Zi *et al.* [67]. During this process, TGFβ1 binding to TGFBR2 induces the recruitment of TGFBR1,
<figure><img src="image_14.png" /><figcaption>**Fig 7. Renewal of quiescent HSCs.** An agent *HSC* with its site *cell_state* in a state *quiescent* and its site *TGFB1_bind* in a state *free* disappears (degradation). Conversely, the second rule describes the creation of this same agent. The rate of application of each rule is shown on the right.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g007

<figure><img src="image_15.png" /><figcaption>**Fig 8. Binding of TGFβ1.** Binding of TGFβ1 (as a token *TGFB1_free*) on the site *cell_state* of an *HSC* agent in a *quiescent* state, its *TGFB1_bind* site is in a *free* state and its *TGFB1_bind* site is in a *Membrane* state. When the rule is fired, the value of the token *TGFB1_free* is decreased and the state of the site *TGFB1_bind* changes to *bound* and the state of the site *TGFB1_bind* changes to *Internalized*.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g008

leading to the activation of signaling pathways that regulate TGFβ1-dependent processes such as COL1 expression and cell proliferation. The trafficking of the TGFBR1/TGFBR2 complex between membrane and cytosol takes 30 min and is not affected by TGFβ1 binding. In the present study, we used *TGFBR* to represent the TGFBR1-TGFBR2 complex and we did not describe the intracellular signaling mechanisms. Upon TGFβ1 binding, the receptor *TGFBR* is internalized to transduce signals and can be either recycled (Fig 9A) or degraded (Fig 9B). To

<figure><img src="image_16.png" /><figcaption>**Fig 9. Degradation and recycling of the TGFβ1 receptor.** (A) The *TGFBR* site in an *Internalized* state changes into a *Membrane* state for the recycling rule, (B) and into a *Degraded* state for the degradation rule.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g009
<figure><img src="image_17.png" /><figcaption>**Fig 10.** HSC activation process. The rule specifies that as long as an agent *HSC*, whose site *cell_state* is in the state *activated* has its *counter_control_counter* with a value lower than or equal to eight then, the value of its *counter_control_counter* will be decremented by one and that of its *counter_intermediate_step* will be incremented by one. This simple rule describes the seven stages of activation.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g010

reduce model and time of calculation, the basal turnover of the receptor is not included in the absence of TGFβ1. To model TGFBR synthesis and localization to the membrane, a rule similar to the rule for degradation is used, changing the site *TGFBR* from the state *degraded* to *membrane*.

**HSC activation and differentiation.** During activation and differentiation, cells have the potential to proliferate and produce COL1 (in variable quantities depending on their activation and differentiation state). However, as each agent can only undertake one action at a time, different rules to describe cell activation, differentiation, proliferation and COL1 production are used.

To model HSC activation and differentiation into MFBs, *counters* are used ranging from 0 to 14, corresponding to the 14 days reported in *in vitro* experiments (0 to 7 corresponding to the change from qHSC to aHSC; 8 to 14 corresponding to the changes from aHSC to MFB). On one hand, two rules describe the change of cell state (qHSC to aHSC and aHSC to MFB) and on other hand, two rules describe the step-by-step dynamics by increasing *counter intermediate step* and decreasing *counter control counter*.

Increasing the *counter intermediate step* allows us to follow the activation and differentiation processes. The decrease of the *counter control counter* ensures that the value of the *counter intermediate step* remains in the range from 0 to 14 (needed because the “<” symbol is not supported in Kappa). Fig 10 describes the dynamics of the activation, defined in 7 steps, corresponding to the seven days of the transition from qHSC to aHSC, with the implementation of a single rule. For the reactivation process, rules similar to those of the activation process were used.

**Cell proliferation.** A proliferation rule is written for each *counter* value. If an *HSC* agent proliferates, all the sites of the newly created agent will have the same value as those of the parent cell, except for TGFBR. Fig 11 shows an example of proliferation rules for the agent *HSC* whose site *cell_state* is in the *activated* state and its *counter intermediate step* is equal to 7. The activation or differentiation steps do not affect the proliferation rate of agents and aHSCs and

<figure><img src="image_18.png" /><figcaption>**Fig 11.** Example of proliferation rules. For an agent *HSC* whose site *cell_state* is in the state *activated* and its *counter control_counter* equal to 7.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g011
<figure><img src="image_19.png" /><figcaption>**Fig 13.** Rules for eliminating MFBs by apoptosis/senescence and inactivation. A), The apoptosis/senescence rule consists of changing the site *cellular_state* of an agent *HSC* from a state *MFB* to *apoptosis_senescence*. B), For inactivation, two rules are necessary, the first changes the site *cellular_state* of an agent *HSC* from a state *MFB* to *inactivated* while reducing the value of the *counter intermediate_step* by 14. The second rule increases the values of *control_counter* by 13 and *intermediate_step* by 1.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011858.g013

(Fig 13A). An additional rule describes the process of elimination of apop_sene_MFB. The dynamics of the inactivation process is described in two steps using two rules (Fig 13B). The first rule changes the MFBs to an inactivated cell state (iHSC) and decrements the counter value named *intermediate_step* to 0. The second rule increments the value of this same *counters* by 1 and increments the value of the *control_counter* by 13.

## Calibration and parameters estimation

For all the biological processes included in the families of models, 13 parameters were obtained directly from the literature, 9 were calculated and 22 parameters were estimated from biological observations by comparing model simulations with experimental data, running hundreds of simulations (S2 Table). Special attention was paid to calibrate these 22 parameters because of the undocumented information about them, but also because of the interdependence between the parameters and the stochastic approach used.

To calibrate the models, we used a block-based parameter estimation method, i.e. a set of parameters for each subsystem (block) was estimated sequentially. Each set was defined in terms of a specific biological process such as cell proliferation, cell inactivation, collagen production and degradation. Block estimation reduced the impact of interdependence and enabled us to find a set of parameters corresponding to biological observations. As the evaluation of collagen was the experimental model best described, we first calibrated the set of parameters related to on collagen-producing cells. Next, we identified the set of parameters related to *react-MFB inactivation* which strongly affects model dynamics. The third set of parameters focused on iHSC and qHSC dynamics. The final set of parameters focused on collagen itself, encompassing its production and degradation. Following this approach, the estimation of the first set of parameters excluded the family of models called *reactMFB-wo-***
inactivation, as shown in Fig 2A. Finally, we performed a global verification to determine whether all the parameters matched the experimental observations (visual inspection).

In addition, the behavior of iHSCs remains unknown, and we developed 3 different families of models taking into account different hypotheses: 1) model family *reactMFB-wo-inactivation* where reactivated MFBs cannot be inactivated, 2) model family *reactMFB-with-inactivation* where reactivated MFBs can be inactivated, 3) model family *iHSC-reversion-to-qHSC* where iHSCs can revert into qHSCs, in the latter model family the *react_MFBs* can be inactivated. Here, we described the parameters based on biological observations. Note that the rule application rate is drawn according to an exponential distribution whose parameter is equal to the sum of the propensities of all the potential events in the system state. The rates of the rules must therefore be defined as the value $T_{1/2}$ of the biological process they describe.

**Parameters for activation and differentiation rules.** The activation process has been extensively documented using primary cultures of HSCs. Isolated qHSCs are spontaneously activated into aHSCs when cultured on plastic for 7 days, and complete transformation of qHSCs into MFBs takes place within 14 days, with a faster proliferation rate for aHSCs than for MFBs. [70, 71]. Based on these dynamics, the models was calibrated with 15 steps, and the proliferation rate of aHSCs was twice that of MFBs. *In vivo*, HSC activation is mainly driven by TGFβ1 which promotes proliferation and COL1 production. While qHSCs are highly sensitive to TGFβ1, MFBs are much less so [71, 72]. Consequently, we did not introduce rules for TGFβ1-dependent proliferation of MFBs in the models, and TGFβ1-dependent collagen production was five times lower in MFBs than in aHSCs.

**Parameters for TGFβ-dependent rules.** The kinetics of TGFβ1 binding to HSCs have already been reported using *in vitro* cell culture models [73] and a slower binding rate to MFBs has been described due to a decrease in the number of receptors [72]. The parameters of the receptor dynamics are those described by Vilar *et al.* [74] and the half-time required for TGFβ1 to induce a signal initiating the activation process was calculated using the model of Zi *et al.* [67]. This initialization step is necessary to activate the cellular machinery before the activation process begins. Note that this initialization time has been divided by 3.5 for iHSCs, as these cells are much more reactive than qHSCs in response to TGFβ1 [24, 25, 38]. Consistent with this, *react_HSCs* were allowed to proliferate immediately after being reactivated, whereas proliferation of aHSCs started at the third step of the activation process. This modeling choice was made to simulate the time required for qHSCs to initiate the cellular machinery of activation [75] whereas *react_HSCs*, due to their activation history, proliferated from the first step of their reactivation.

**Parameters for qHSC self-renewal rate and iHSC elimination rules.** In the absence of information on qHSC self-renewal rate and iHSC half-life, these parameters were estimated on the basis of data reported by Kisseleva *et al.* [24].

**Parameters for MFB and react_MFB elimination rules.** The half-life of MFBs and *react_MFBs* has not been documented, and was estimated on the basis of observations of the accumulation of these cells during TGFβ1 stimulation and their disappearance 1 to 2 months after the last TGFβ1 stimulation (Kisseleva *et al.* [24]). This accumulation represents an increase in the total number of HSCs after two months of CCl4 treatment from 10.6 to 14.3% of the total liver cell population, where qHSCs and *alpha*-SMA cells account for 10.6% (±0.8) and 14.3% (±1.5), respectively.

**Parameters for collagen 1 dynamic rules.** The rules governing the production, degradation, stabilization, and destabilization of *COL1_high_remodeling* and *COL1_low_remodeling* were calibrated using data from Kisseleva *et al.* [24]. Since COL1 is experimentally quantified as deposits resulting from global remodeling, total COL1 (COL1_tot) is defined as the sum of *COL1_remodeling_high* and *COL1_remodeling_low*. COL1_tot increased 12- and 14-fold after
1 and 2 months of TGFβ1 stimulation (two stimulation per week), respectively. aHSCs and react_HSCs produced 10 time less COL1_remodeling_high than MFBs and react_MFBs. However, TGFβ1-dependent COL1 production was increased 10-fold in aHSCs and react_HSCs and only 2-fold in MFBs and react_MFBs, since aHSCs and react_HSCs are more sensitive to TGFβ1 [24, 25, 38]. In the absence of TGFβ1 stimulation, the models were calibrated to observe a rapid decrease in COL1_tot over the first month, followed by a gradual return to its initial value around 6 months later.

**Parameters for the model families reactMFB-wo-inactivation and reactMFB-with-inactivation.** Upon removal of liver injury, i.e. in the absence of TGFβ1 stimulation in the models, 45% of MFBs are inactivated and the remaining 55% are eliminated through apoptosis and senescent pathways [25, 76]. iHSCs can be reactivated upon TGFβ1 stimulation to produce react_MFBs, but no information is available on their fate. Troeger *et al.* hypothesized that react_MFBs can undergo new inactivation/reactivation cycles, but should be limited for each cell [25]. Based on this hypothesis, two model families were created, including model family reactMFB-wo-inactivation in which react_MFBs cannot be inactivated and model family reactMFB-with-inactivation in which inactivation is possible. For the latter, we developed models by varying the percentage of reactivated MFBs that could be inactivated from 1 to 50%.

**Parameters for the model iHSC-reversion-to-qHSC.** Using mouse models of liver fibrosis induced by CCl4, Kisseleva *et al.* [24] assessed the amount of iHSCs compared with qHSCs after removal of CCl4. They reported that qHSCs and iHSCs were in equal proportion after one month of recovery, and that the sum of qHSCs and iHSCs was equal to the initial number of qHSCs after six months of recovery. However, there is no information about the behavior of iHSCs in absence of CCl4 (i.e. of TGFβ1 in the models). In the absence of information on how iHSCs behave in the absence of CCl4 (i.e. of TGFβ1 in the models), we developed a family of models with a possible return to the quiescent state by varying the ratio of iHSC that transform into qHSC from 0 to 100%.

## Experimental validation using a mouse model of CCl4-induced liver fibrosis

To validate the predictions of the model, collagen I was quantified during the course of fibrosis and reversion using a mouse model of CCl4 induced liver fibrosis. As previously described [77], seven-week-old female C57Bl/6 mice were treated with oral administration of CCl4 (Sigma-Aldrich, St. Louis, MO, USA) diluted in olive oil. A first dose of 2.4 g/kg of mouse weight was administered to mice three days before starting weekly treatment with a 1.6 g/kg dose for ten weeks. Control mice were treated with the vehicle only (olive oil). Mice were sacrificed at 24 hrs after the last CCl4 dose. For fibrosis reversal experiments, mice were injected for 4 weeks with CCl4 and sacrificed at 8 days (1 week of recovery), or 15 days (2 weeks of recovery) after the last injection. Each group (control, fibrosis and reversion) at each time of sacrifice contained 5 mice. Collagen quantification was performed by SHG microscopy on 20 µm frozen tissue sections as previously described [78].

## Computational validation using Gene set enrichment analysis of RNA-Seq data from patients with liver fibrosis

The enrichment for iHSC gene expression signature in human fibrotic samples was performed using Gene Set Enrichment Analysis (GSEA) [79, 80]. We used the iHSC gene expression signature previously identified by Rosenthal *et al.* [28] in a mouse model of non-alcoholic steatohepatitis (NASH) and we added iHSC markers previously identified in a mouse model of CCl4 induced liver fibrosis [24] and an *in vitro* reversion model of human activated HSCs [38]. The
list of genes is detailed in supplementary S1 Table. For clinical samples, RNAseq data from patients with non-alcoholic fatty liver disease (NAFLD) [39] were used. Note that three genes (KRT20, GABRA3 and GSTT1) were not identified in the RNAseq data due to a low expression level. Liver samples with fibrosis were selected and separated in two groupss: with low (F0-F1, n = 34) and high (F3-F4, n = 68) fibrosis grade. RNAseq data were filtered using the BIOMART database [81] to select genes encoding proteins and data were normalized using DESeq2 [82] with a threshold of 100 reads. The filtration steps reduced the number of genes from 21,595 to 16,021.

**Supporting information**

**S1 Fig. Temporal analysis of cells and COL1 in models with or without reversion of iHSC to a quiescent state.** Simulation series were performed using conditions of stimulation from Kisseleva et al. [24]. TGFβ1 parameters were as follows 10,000 molecules per cell, 16 stimulation (twice a week) during 2 months. (A1), models *reactMFB-with-inactivation* with 0% of the iHSCs reverting into qHSC, (A2) and (A3), models *iHSC-reversion-to-qHSC* with 50% and 100% of the iHSCs reverting into qHSC, the remaining iHSCs being eliminated. (B1, B2 and B3), variation in the number of cell occurrences for qHSCs, aHSCs, MFBs, iHSCs, react_HSCs, react_MFBs and α-SMA positive cells in the three families of models. The number of α-SMA cells is the sum of the number of aHSCs, MFBs, react_HSCs and react_MFBs. Simulations are expressed as the mean of 10 replicates. (C1, C2 and C3), variation in the number of COL1 occurrence in the three models. Data are expressed as arbitrary unit (a.u) and all 5 replicates are represented.
(TIF)

**S2 Fig. Comparison of model predictions and experimental data acquired in models of fibrosis induced by dimethyl_nitrosamine, thioacetamide and a high-fat diet.** Simulations were carried out using the stimulation protocols described in the different models. (A) dimethyl_nitrosamine model: 3 TGFβ1 stimuli per week for 6 weeks [36], (B) thioacetamide model: 2 simulations with TGFβ1 per week for 10 weeks [37], (C) high-fat model: 6 simulations with TGFβ1 every 14 days for 84 days, adapted to fit Farooq et al. results [55]. Collagen accumulation was plotted and experimental data from [36, 37, 55] are indicated with black cross.
(TIF)

**S1 Table. List of the genes used for the Gene Set Enrichment Analysis.** The first column contains the genes from Rosenthal et al. It includes 39 genes that are differentially expressed in inactivated stellate cells compared with activated stellate cells and quiescent stellate cells [28]. The second list contains the eight marker genes of inactivated stellate cells identified both in a mouse model of CCl4-induced liver fibrosis [24] and in a *in vitro* reversion model of human activated cells [38].
(XLSX)

**S2 Table. Table of parameters.** A.U: Arbitrary Unit corresponding to % areas of COL1 deposits, *Calculated*: parameter calculated directly from values found in the bibliography. *Estimated from biological observations*: parameters for which the values are estimated empirically in order to adapt the families of models to observations at a specific time. All the parameters describing a duration are calculated as the half-time of their reactions. Using the exponential law, the average dynamics of reactions corresponding to the experimental observations was ensured [83] using Block based parameter estimation.
(XLSX)
## Acknowledgments

The authors thank the excellent support of the GenOuest bioinformatics core facility, Dr G Baffet and F Ezan for biological experiments and Dr C Lavau (Inserm U1085, University of Rennes) for proofreading the manuscript.

## Author Contributions

**Conceptualization:** Matthieu Bouguéon, Vincent Legagneux, Octave Hazard, Jérôme Feret, Nathalie Théret.

**Data curation:** Matthieu Bouguéon, Vincent Legagneux, Jérôme Feret, Nathalie Théret.

**Formal analysis:** Matthieu Bouguéon, Jérôme Feret, Nathalie Théret.

**Funding acquisition:** Anne Siegel, Jérôme Feret, Nathalie Théret.

**Investigation:** Matthieu Bouguéon, Jérémy Bomo, Nathalie Théret.

**Methodology:** Matthieu Bouguéon, Vincent Legagneux, Jérôme Feret, Nathalie Théret.

**Project administration:** Anne Siegel, Jérôme Feret, Nathalie Théret.

**Resources:** Nathalie Théret.

**Software:** Matthieu Bouguéon.

**Supervision:** Anne Siegel, Jérôme Feret, Nathalie Théret.

**Validation:** Matthieu Bouguéon, Vincent Legagneux, Nathalie Théret.

**Visualization:** Matthieu Bouguéon, Vincent Legagneux, Nathalie Théret.

**Writing – original draft:** Matthieu Bouguéon, Nathalie Théret.

**Writing – review & editing:** Matthieu Bouguéon, Vincent Legagneux, Anne Siegel, Jérôme Feret, Nathalie Théret.

## References

1. Sepanlou SG, Safiri S, Bisignano C, Ikuta KS, Merat S, Saberifiroozi M, et al. The global, regional, and national burden of cirrhosis by cause in 195 countries and territories, 1990-2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet Gastroenterol Hepatol. 2020; 5(3):245-266. https://doi.org/10.1016/S2468-1253(19)30349-8

2. Tsuchida T, Friedman SL. Mechanisms of hepatic stellate cell activation. Nat Rev Gastroenterol Hepatol. 2017; 14(7):397-411. https://doi.org/10.1038/nrgastro.2017.38 PMID: 28487545

3. Cordero-Espinoza L, Huch M. The balancing act of the liver: tissue regeneration versus fibrosis. J Clin Invest. 2018; 128(1):85-96. https://doi.org/10.1172/JCI93562 PMID: 29293095

4. Kitto LJ, Henderson NC. Hepatic Stellate Cell Regulation of Liver Regeneration and Repair. Hepatol Commun. 2021; 5(3):358-370. https://doi.org/10.1002/hep4.1628 PMID: 33681672

5. Kisseleva T, Brenner D. Molecular and cellular mechanisms of liver fibrosis and its regression. Nature Reviews Gastroenterology & Hepatology. 2021; 18(3):151-166. https://doi.org/10.1038/s41575-020-00372-7 PMID: 33128017

6. Campana L, Iredale JP. Regression of liver fibrosis. In: Seminars in liver disease. vol. 58. Thieme Medical Publishers; 2017. p. 001-010.

7. Southern J, Pitt-Francis J, Whiteley J, Stokeley D, Kobashi H, Nobes R, et al. Multi-scale computational modelling in biology and physiology. Progress in biophysics and molecular biology. 2008; 96(1-3):60-89. https://doi.org/10.1016/j.pbiomolbio.2007.07.019 PMID: 17888502

8. Versypt ANF. Multiscale modeling in disease. Current opinion in systems biology. 2021; 27:100340. https://doi.org/10.1016/j.coisb.2021.05.001
9. Gondal MN, Chaudhary SU. Navigating Multi-Scale Cancer Systems Biology Towards Model-Driven Clinical Oncology and Its Applications in Personalized Therapeutics. Frontiers in Oncology. 2021; p. 4767. https://doi.org/10.3389/fonc.2021.712505 PMID: 34900668

10. Amar MB, Bianca C. Towards a unified approach in the modeling of fibrosis: A review with research perspectives. Physics of life reviews. 2016; 17:61–85. https://doi.org/10.1016/j.plrev.2016.03.005 PMID: 27079617

11. Zeigler AC, Richardson WJ, Holmes JW, Saucerman JJ. Computational modeling of cardiac fibroblasts and fibrosis. Journal of molecular and cellular cardiology. 2016; 93:73–83. https://doi.org/10.1016/j.yjmcc.2015.11.020 PMID: 26608708

12. Leonard-Duke J, Evans S, Hannan RT, Barker TH, Bates JH, Bonham CA, et al. Multi-scale models of lung fibrosis. Matrix biology. 2020; 91:35–50. https://doi.org/10.1016/j.matbio.2020.04.003 PMID: 32438056

13. Dutta-Moscato J, Solovyev A, Mi Q, Nishikawa T, Soto-Gutierrez A, Fox IJ, et al. A multiscale agent-based in silico model of liver fibrosis progression. Frontiers in bioengineering and biotechnology. 2014; 2:18. https://doi.org/10.3389/fbioe.2014.00018 PMID: 25152891

14. Yoshizawa M, Sugimoto M, Tanaka M, Sakai Y, Nishikawa M. Computational simulation of liver fibrosis dynamics. Scientific reports. 2022; 12(1):14112. https://doi.org/10.1038/s41598-022-18123-w PMID: 35982187

15. Wang Y, Jiang J. A two-dimensional (2d) systems biology-based discrete liver tissue model: A simulation study with implications for ultrasound elastography of liver fibrosis. Computers in biology and medicine. 2019; 104:227–234. https://doi.org/10.1016/j.compbiomed.2018.11.027 PMID: 30529712

16. Friedman A, Hao W. Mathematical modeling of liver fibrosis. Mathematical Biosciences & Engineering. 2016; 14(1):143–164.

17. Krenkel O, Hundertmark J, Ritz TP, Weiskirchen R, Tacke F. Single cell RNA sequencing identifies subsets of hepatic stellate cells and myofibroblasts in liver fibrosis. Cells. 2019; 8(5):503. https://doi.org/10.3390/cells8050503 PMID: 31137713

18. Terkelsen MK, Bendixen SM, Hansen D, Scott EAH, Moeller AF, Nielsen R, et al. Transcriptional Dynamics of Hepatic Sinusoid-Associated Cells After Liver Injury. Hepatology. 2020; 72(6):2119–2133. https://doi.org/10.1002/hep.31215 PMID: 32145072

19. Payen VL, Lavergne A, Alevra Sarika N, Colonval M, Karim L, Deckers M, et al. Single-cell RNA sequencing of human liver reveals hepatic stellate cell heterogeneity. JHEP Rep. 2021; 3(3):100278. https://doi.org/10.1016/j.jhepr.2021.100278 PMID: 34027339

20. Boutillier P, Maasha M, Li X, Medina-Abarca HF, Krivine J, Feret J, et al. The Kappa platform for rule-based modeling. Bioinformatics. 2018; 34(13):i583–i592. https://doi.org/10.1093/bioinformatics/bty272 PMID: 29950016

21. Danos V, Laneve C. Graphs for core molecular biology. In: International Conference on Computational Methods in Systems Biology. Springer; 2003. p. 34–46.

22. Danos V, Laneve C. Formal molecular biology. Theoretical Computer Science. 2004; 325(1):69–110. https://doi.org/10.1016/j.tcs.2004.03.065

23. Dewidar B, Meyer C, Dooley S, Meindl-Beinker AN. in Hepatic Stellate Cell Activation and Liver Fibrogenesis-Updated 2019. Cells. 2019; 8(11). https://doi.org/10.3390/cells8111419 PMID: 31718044

24. Kisseleva T, Cong M, Paik Y, Scholten D, Jiang C, Benner C, et al. Myofibroblasts revert to an inactive phenotype during regression of liver fibrosis. Proc Natl Acad Sci U S A. 2012; 109(24):9448–9453. https://doi.org/10.1073/pnas.1201840109 PMID: 22566629

25. Troeger JS, Mederacke I, Gwak GY, Dapito DH, Mu X, Hsu CC, et al. Deactivation of hepatic stellate cells during liver fibrosis resolution in mice. Gastroenterology. 2012; 143(4):1073–1083. https://doi.org/10.1053/j.gastro.2012.06.036 PMID: 22750464

26. Danos V, Feret J, Fontana W, Krivine J. Scalable simulation of cellular signaling networks. In: Asian Symposium on Programming Languages and Systems. Springer; 2007. p. 139–157.

27. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. The journal of physical chemistry. 1977; 81(25):2340–2361. https://doi.org/10.1021/j100540a008

28. Rosenthal SB, Liu X, Ganguly S, Dhar D, Pasillas MP, Ricciardelli E, et al. Heterogeneity of HSCs in a Mouse Model of NASH. Hepatology. 2021; 74(2):667–685. https://doi.org/10.1002/hep.31743 PMID: 33550587

29. Scholten D, Trebicka J, Liedtke C, Weiskirchen R. The carbon tetrachloride model in mice. Lab Anim. 2015; 49(1 Suppl):4–11. https://doi.org/10.1177/0023677215571192 PMID: 25835733

30. Wu S, Wang X, Xing W, Li F, Liang M, Li K, et al. An update on animal models of liver fibrosis. Front Med (Lausanne). 2023; 10:1160053. https://doi.org/10.3389/fmed.2023.1160053 PMID: 37035335
31. Yanguas SC, Cogliati B, Willebrords J, Maes M, Colle I, van den Bossche B, et al. Experimental models of liver fibrosis. Arch Toxicol. 2016; 90(5):1025–1048. https://doi.org/10.1007/s00204-015-1543-4 PMID: 26047667

32. Gailhouste L, Le Grand Y, Odin C, Guyader D, Turlin B, Ezan F, et al. Fibrillar collagen scoring by second harmonic microscopy: a new tool in the assessment of liver fibrosis. J Hepatol. 2010; 52(3):398–406. https://doi.org/10.1016/j.jhep.2009.12.009 PMID: 20149472

33. Bourd-Boittin K, Bonnier D, Leyme A, Mari B, Tuffery P, Samson M, et al. Protease profiling of liver fibrosis reveals the ADAM metallopeptidase with thrombospondin type 1 motif, 1 as a central activator of transforming growth factor beta. Hepatology. 2011; 54(6):2173–2184. https://doi.org/10.1002/hep.24598 PMID: 21826695

34. Iwaisako K, Brenner DA, Kisseleva T. What's new in liver fibrosis? The origin of myofibroblasts in liver fibrosis. J Gastroenterol Hepatol. 2012; 27 Suppl 2(Suppl 2):65–68. https://doi.org/10.1111/j.1440-1746.2011.07002.x PMID: 22320919

35. Yang W, He H, Wang T, Su N, Zhang F, Jiang K, et al. Single-Cell Transcriptomic Analysis Reveals a Hepatic Stellate Cell-Activation Roadmap and Myofibroblast Origin During Liver Fibrosis in Mice. Hepatology. 2021; 74(5):2774–2790. https://doi.org/10.1002/hep.31987 PMID: 34089528

36. Sato Y, Yoneda A, Shimizu F, Nishimura M, Shimoyama R, Tashiro Y, et al. Resolution of fibrosis by siRNA HSP47 in vitamin A-coupled liposomes induces regeneration of chronically injured livers. J Gastroenterol Hepatol. 2021; 36(12):3418–3428. https://doi.org/10.1111/jgh.15587 PMID: 34151462

37. Abramovitch S, Dahan-Bachar L, Sharvit E, Weisman Y, Ben Tov A, Brazowski E, et al. Vitamin D inhibits proliferation and profibrotic marker expression in hepatic stellate cells and decreases thioacetamide-induced liver fibrosis in rats. Gut. 2011; 60(12):1728–1737. https://doi.org/10.1136/gut.2010.234666 PMID: 21816960

38. El Taghdouini A, Najimi M, Sancho-Bru P, Sokal E, van Grunsven LA. In vitro reversion of activated primary human hepatic stellate cells. Fibrogenesis & tissue repair. 2015; 8(1):1–15. https://doi.org/10.1186/s13069-015-0031-z PMID: 26251672

39. Govaere O, Cockell S, Tiniakos D, Queen R, Younes R, Vacca M, et al. Transcriptomic profiling across the nonalcoholic fatty liver disease spectrum reveals gene signatures for steatohepatitis and fibrosis. Sci Transl Med. 2020; 12(572). https://doi.org/10.1126/scitranslmed.aba4448 PMID: 33268509

40. Trivedi P, Wang S, Friedman SL. The Power of Plasticity-Metabolic Regulation of Hepatic Stellate Cells. Cell Metab. 2021; 33(2):242–257. https://doi.org/10.1016/j.cmet.2020.10.026 PMID: 33232666

41. Kordes C, Bock HH, Reichert D, May P, ussinger D. Hepatic stellate cells: current state and open questions. Biol Chem. 2021; 402(9):1021–1032. https://doi.org/10.1515/hsz-2021-0180 PMID: 34008380

42. Chylek LA, Harris LA, Tung CS, Faeder JR, Lopez CF, Hlavacek WS. Rule-based modeling: a computational approach for studying biomolecular site dynamics in cell signaling systems. Wiley Interdiscip Rev Syst Biol Med. 2014; 6(1):13–36. https://doi.org/10.1002/wsbm.1245 PMID: 24123887

43. Dematté L, Priami C, Romanel A. The BlenX Language: A Tutorial. In: Bernardo M, Degano P, Zavattaro G, editors. Formal Methods for Computational Systems Biology: 8th International School on Formal Methods for the Design of Computer, Communication, and Software Systems, SFM 2008 Bertinoro, Italy, June 2-7, 2008 Advanced Lectures. Berlin, Heidelberg: Springer Berlin Heidelberg; 2008. p. 313–365. Available from: http://dx.doi.org/10.1007/978-3-540-68894-5_9

44. Ciocchetta F, Hillston J. Bio-PEPA: A framework for the modelling and analysis of biological systems. Theoretical Computer Science. 2009; 410(33–34):3065–3084. https://doi.org/10.1016/j.tcs.2009.02.037

45. Harris LA, Hogg JS, Tapia JJ, Sekar JA, Gupta S, Korsunsky I, et al. BioNetGen 2.2: advances in rule-based modeling. Bioinformatics. 2016; 32(21):3366–3368. https://doi.org/10.1093/bioinformatics/btw469 PMID: 27402907

46. Danos V, Feret J, Fontana W, Harmer R, Krivine J. Rule-based modelling of cellular signalling. In: International conference on concurrency theory. Springer; 2007. p. 17–41.

47. Petrov T, Igler C, Sezgin A, Henzinger TA, Guet CC. Long lived transients in gene regulation. Theoretical Computer Science. 2021; 893:1–16. https://doi.org/10.1016/j.tcs.2021.05.023

48. Krivine J, Danos V, Benecke A. Modelling Epigenetic Information Maintenance: A Kappa Tutorial. In: Bouajjani A, Maler O, editors. Computer Aided Verification. Berlin, Heidelberg: Springer Berlin Heidelberg; 2009. p. 17–32.

49. Wilson-Kanamori J, Danos V, Thomson T, Honorato-Zimmer R. Kappa rule-based modeling in synthetic biology. Methods Mol Biol. 2015; 1244:105–135. https://doi.org/10.1007/978-1-4939-1878-2_6 PMID: 25487095

50. Köhler A, Krivine J, Vidmar J. A Rule-Based Model of Base Excision Repair. In: Mendes P, Dada JO, Smallbone K, editors. Computational Methods in Systems Biology. Cham: Springer International Publishing; 2014. p. 173–195.
51. Theret N, Feret J, Hodgkinson A, Boutillier P, Vignet P, Radulescu O. Integrative models for TGF-b signaling and extracellular matrix. Extracellular Matrix Omics. 2020; p. 209–225. https://doi.org/10.1007/978-3-030-58330-9_10

52. Bougéon M, Boutillier P, Feret J, Hazard O, Théret N. The Rule-Based Model Approach: A Kappa Model for Hepatic Stellate Cells Activation by TGFB1. Systems Biology Modelling and Analysis: Formal Bioinformatics Methods and Tools. 2022; p. 69–126.

53. Waites W., Cavaliere M., Danos V., Datta R., Eggo R. M., Hallett T. B., Manheim D., Panovska-Griffiths J., Russell T. W. and Zarnitsyna V. I. Compositional modelling of immune response and virus transmission dynamics Philosophical Transactions of the Royal Society A, 2022, vol. 380, no 2233, p. 20210307. https://doi.org/10.1098/rsta.2021.0307 PMID: 35965463

54. Wei G, An P, Vaid KA, Nasser I, Huang P, Tan L, et al. Comparison of murine steatohepatitis models identifies a dietary intervention with robust fibrosis, ductular reaction, and rapid progression to cirrhosis and cancer. Am J Physiol Gastrointest Liver Physiol. 2020; 318(1):G174–G188. https://doi.org/10.1152/ajpgi.00041.2019

55. Farooq M, Hameed H, Dimanche-Boitrel MT, Piquet-Pellorce C, Samson M, Le Seyec J. Switching to Regular Diet Partially Resolves Liver Fibrosis Induced by High-Fat, High-Cholesterol Diet in Mice. Nutrients. 2022; 14(2). https://doi.org/10.3390/nu14020386 PMID: 35057565

56. Zhang W, Conway SJ, Liu Y, Snider P, Chen H, Gao H, et al. Heterogeneity of hepatic stellate cells in fibrogenesis of the liver: insights from single-cell transcriptomic analysis in liver injury. Cells. 2021; 10(8):2129. https://doi.org/10.3390/cells10082129 PMID: 34440898

57. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature biotechnology. 2014; 32(4):381–386. https://doi.org/10.1038/nbt.2859 PMID: 24658644

58. Street K, Risso D, Fletcher RB, Das D, Ngai J, Yosef N, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics. 2018; 19(1):477. https://doi.org/10.1186/s12864-018-4772-0 PMID: 29914354

59. Koda Y, Teratani T, Chu PS, Hagihara Y, Mikami Y, Harada Y, et al. CD8+ tissue-resident memory T cells promote liver fibrosis resolution by inducing apoptosis of hepatic stellate cells. Nature communications. 2021; 12(1):4474. https://doi.org/10.1038/s41467-021-24734-0 PMID: 34294714

60. Ramachandran P, Pellicoro A, Vernon MA, Boulter L, Aucott RL, Ali A, et al. Differential Ly-6C expression identifies the recruited macrophage phenotype, which orchestrates the regression of murine liver fibrosis. Proc Natl Acad Sci U S A. 2012; 109(46):E3186–E3195. https://doi.org/10.1073/pnas.1119964109 PMID: 23100531

61. Duffield JS, Forbes SJ, Constandinou CM, Clay S, Partolina M, Vuthoori S, et al. Selective depletion of macrophages reveals distinct, opposing roles during liver injury and repair. J Clin Invest. 2005; 115(1):65–65. https://doi.org/10.1172/JCI22675 PMID: 15630444

62. Danos V, Laneve C. Formal molecular biology. Theoretical Computer Science. 2004; 325(1):69–110. https://doi.org/10.1016/j.tcs.2004.03.065

63. Faeder JR, Blinov ML, Goldstein B, Hlavacek WS. Rule-based modeling of biochemical networks. Complexity. 2005; 10(4):22–41. https://doi.org/10.1002/cplx.20074

64. Andrei O, Kirchner H. A Rewriting Calculus for Multigraphs with Ports. Electr Notes Theor Comput Sci. 2008; 219:67–82. https://doi.org/10.1016/j.entcs.2008.10.035

65. John M, Lhoussaine C, Niehren J, Versari C. Biochemical Reaction Rules with Constraints. In: Programming Languages and Systems—20th European Symposium on Programming, ESOP 2011, Held as Part of the Joint European Conferences on Theory and Practice of Software, ETAPS 2011, Saarbrücken, Germany, March 26-April 3, 2011. Proceedings. vol. 6602 of Lecture Notes in Computer Science. Springer; 2011. p. 338–357. Available from: http://dx.doi.org/10.1007/978-3-642-19718-5

66. Di Guglielmo GM, Le Roy C, Goodfellow AF, Wrana JL. Distinct endocytic pathways regulate TGF-b receptor signalling and turnover. Nature cell biology. 2003; 5(5):410–421. https://doi.org/10.1038/ncb975 PMID: 12717440

67. Zi Z, Feng Z, Chapnick DA, Dahl M, Deng D, Klipp E, et al. Quantitative analysis of transient and sustained transforming growth factor-β signaling dynamics. Molecular systems biology. 2011; 7(1):492. https://doi.org/10.1038/msb.2011.22 PMID: 21613981

68. Krizhanovsky V, Yon M, Dickins RA, Hearn S, Simon J, Miething C, et al. Senescence of activated stellate cells limits liver fibrosis. Cell. 2008; 134(4):657–667. https://doi.org/10.1016/j.cell.2008.06.049 PMID: 18724938

69. Schnabl B, Purbeck CA, Choi YH, Hagedorn CH, Brenner D. Replicative senescence of activated human hepatic stellate cells is accompanied by a pronounced inflammatory but less fibrogenic phenotype. Hepatology. 2003; 37(3):653–664. https://doi.org/10.1053/jhep.2003.50097 PMID: 12601363
70. Bachem MG, Meyer D, Melchior R, Sell KM, Gressner AM. Activation of rat liver perisinusoidal lipocytes by transforming growth factors derived from myofibroblastlike cells. A potential mechanism of self perpetuation in liver fibrogenesis. J Clin Invest. 1992; 89(1):19-27. https://doi.org/10.1172/JCI115561 PMID: 1729271

71. Bachem MG, Meyer D, Schäfer W, Riess U, Melchior R, Sell KM, et al. The response of rat liver perisinusoidal lipocytes to polypeptide growth regulator changes with their transdifferentiation into myofibroblast-like cells in culture. Journal of hepatology. 1993; 18(1):40-52. https://doi.org/10.1016/S0168-8278(05)80008-X PMID: 8340608

72. Dooley S, Delvoux B, Lahme B, Mangasser-Stephan K, Gressner AM. Modulation of transforming growth factor β response and signaling during transdifferentiation of rat hepatic stellate cells to myofibroblasts. Hepatology. 2000; 31(5):1094-1106. https://doi.org/10.1053/he.2000.6126 PMID: 10796885

73. Friedman SL, Yamasaki G, Wong L. Modulation of transforming growth factor beta receptors of rat lipo-cytes during the hepatic wound healing response. Enhanced binding and reduced gene expression accompany cellular activation in culture and in vivo. Journal of Biological Chemistry. 1994; 269 (14):10551-10558. https://doi.org/10.1016/S0021-9258(17)34095-4 PMID: 8144642

74. Vilar JMG, Jansen R, Sander C. Signal processing in the TGF-β superfamily ligand-receptor network. PLoS computational biology. 2006; 2(1):e3. https://doi.org/10.1371/journal.pcbi.0020003 PMID: 16446785

75. Montosi G, Garuti C, Iannone A, Pietrangelo A. Spatial and temporal dynamics of hepatic stellate cell activation during oxidant-stress-induced fibrogenesis. Am J Pathol. 1998; 152(5):1319-1326. PMID: 9588900

76. Liu X, Xu J, Brenner DA, Kisseleva T. Reversibility of liver fibrosis and inactivation of fibrogenic myofibroblasts. Current pathobiology reports. 2013; 1(3):209-214. https://doi.org/10.1007/s40139-013-0018-7 PMID: 24000319

77. Marvie P, Lisbonne M, L'helgoualc'h A, Rauch M, Turlin B, Preisser L, et al. Interleukin-33 overexpression is associated with liver fibrosis in mice and humans. J Cell Mol Med. 2010; 14(6B):1726-1739. https://doi.org/10.1111/j.1582-4934.2009.00801.x PMID: 19508382

78. Gailhouste L, Le Grand Y, Odin C, Guyader D, Turlin B, Ezan F, et al. Fibrillar collagen scoring by second harmonic microscopy: a new tool in the assessment of liver fibrosis. Journal of hepatology. 2010; 52(3):398-406. https://doi.org/10.1016/j.jhep.2009.12.009 PMID: 20149472

79. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005; 102(43):15545-15550. https://doi.org/10.1073/pnas.0506580102 PMID: 16199517

80. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003; 34(3):267-273. https://doi.org/10.1038/ng1180 PMID: 12808457

81. Martin FJ, Amode MR, Aneja A, Austine-Orimoloye O, Azov AG, Barnes I, et al. Ensembl 2023. Nucleic Acids Res. 2023; 51(D1):D933-D941. https://doi.org/10.1093/nar/gkac958 PMID: 36318249

82. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 2014; 15(12):1-21. https://doi.org/10.1186/s13059-014-0550-8 PMID: 25516281

83. Denning PJ, Lewis TG. Exponential laws of computing growth. Communications of the ACM. 2016; 60 (1):54-65. https://doi.org/10.1145/2976758