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

HAL
open science

# A liver digital twin for in silico testing of cellular and inter-cellular mechanisms in regeneration after drug-induced damage

Jieling Zhao, Ahmed Ghallab, Reham Hassan, Steven Dooley, Jan Georg

Hengstler, Dirk Drasdo

► To cite this version:

Jieling Zhao, Ahmed Ghallab, Reham Hassan, Steven Dooley, Jan Georg Hengstler, et al. A liver digital twin for in silico testing of cellular and inter-cellular mechanisms in regeneration after drug-induced damage. iScience, 2024, 27 (2), pp.108077. 10.1016/j.isci.2023.108077. hal-03738207v2

HAL Id: hal-03738207

https://inria.hal.science/hal-03738207v2

Submitted on 8 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A liver digital twin for in silico testing of cellular and inter-cellular mechanisms in regeneration after drug-induced damage

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

Jieling Zhao,
Ahmed Ghallab,
Reham Hassan,
Steven Dooley, Jan
Georg Hengstler,
Dirk Drasdo

dirk.drasdo@inria.fr

Highlights

3D spatial-temporal digital
twin of liver regeneration
after APAP overdose

Digital twin implementing
a wide range of
hypothesized cell-cell
communications

Digital twin testing
alternative sets of
mechanisms against
regeneration data

Digital twin permitting
simulated predictions
guiding experimental
design
## Article

# A liver digital twin for in silico testing of cellular and inter-cellular mechanisms in regeneration after drug-induced damage

Jieling Zhao,<fnref n="1" />,<fnref n="2" /> Ahmed Ghallab,<fnref n="1" />,<fnref n="3" /> Reham Hassan,<fnref n="1" />,<fnref n="3" /> Steven Dooley,<fnref n="4" /> Jan Georg Hengstler,<fnref n="1" /> and Dirk Drasdo<fnref n="1" />,<fnref n="2" />,<fnref n="5" />,*

## SUMMARY

This communication presents a mathematical mechanism-based model of the regenerating liver after drug-induced pericentral lobule damage resolving tissue microarchitecture. The consequence of alternative hypotheses about the interplay of different cell types on regeneration was simulated. Regeneration dynamics has been quantified by the size of the damage-induced dead cell area, the hepatocyte density and the spatial-temporal profile of the different cell types. We use deviations of observed trajectories from the simulated system to identify branching points, at which the systems behavior cannot be explained by the underlying set of hypotheses anymore. Our procedure reflects a successful strategy for generating a fully digital liver twin that, among others, permits to test perturbations from the molecular up to the tissue scale. The model simulations are complementing current knowledge on liver regeneration by identifying gaps in mechanistic relationships and guiding the system toward the most informative (lacking) parameters that can be experimentally addressed.

## INTRODUCTION

Liver regeneration is an overly complex process involving many different cell types and factors.<fnref n="1" />⁻<fnref n="3" /> Upon inflicted liver damage, the injured main hepatic parenchyma cells, namely, hepatocytes start to release factors to trigger an inflammatory response.<fnref n="4" />,<fnref n="5" /> Liver resident macrophages, namely, Kupffer cells are activated to recruit neutrophils to initiate the death of injured hepatocytes.⁶ After liver injury is initiated, hepatic stellate cells (HSCs) migrate into the lesion and become activated through stimulating factors released from damaged hepatocytes, activated macrophages, the extracellular matrix, neighboring sinusoidal endothelial cells and platelets.<fnref n="3" />,⁷ After liver injury, there is also massive infiltration of bone marrow-derived (non-resident) macrophages into the liver to elicit liver impairment and to restore liver integrity at different stages.⁸

Despite extensive studies over decades, the precise orchestration is only partially understood, and the understanding is mainly qualitative. There are alternative hypotheses on how certain factors interplay (which cannot simultaneously apply), and it is not well known how sensitive the regeneration outcome is, regarding a modification of time scales or rate constants. This is difficult to study experimentally for several reasons: (1) It is not possible to control all parameters that may influence a certain observation in an experimental system, in particular, in vivo, such that the interpretation of experimental observations in the context of complex interactions is very limited. (2) Rate constants and time scales are hard to control in animal experiments. (3) Experiments are resource-intense with regards to personnel, technology, animal numbers, and consumables, limiting the number of experiments that can be performed. (4) The primary motivation for experiments is to understand the regeneration process in human, but experiments in human are per se not permitted and experiments in animal models are subject to close ethical control and do only partly reflect the processes in human.

Hence, a promising strategy to overcome these limitations may be to formulate alternative sets of mechanisms and implement them on the computer within alternative virtual/digital liver twins, each representing a different set of hypothesized mechanisms, whereby hypothesized mechanisms may emerge also from *in vitro* experiments or from other source. Two sets of mechanisms are considered as different if at least one mechanism is different. There is for far no unique definition of the term “digital twin.” Throughout this article, we distinguish between a digital twin (DT) and a digital twin candidate (DTC). By DTC we denote a computational (digital) model (not an animation or a set of digital images) implementing a defined set of mechanisms on the computer that is compatible with literature knowledge (including parameter ranges, physical laws and so forth), optionally complemented with some own experiments to specify the set of hypotheses that serves as input

<footnote n="1">¹Leibniz Research Centre for Working Environment and Human Factors, Technical University of Dortmund (IfADo), 44139 Dortmund, Germany</footnote>

<footnote n="2">²Group SIMBIOTX, INRIA Saclay, 91120 Palaiseau, France</footnote>

<footnote n="3">³Department of Forensic Medicine and Toxicology, Faculty of Veterinary Medicine, South Valley University, Qena 83523, Egypt</footnote>

<footnote n="4">⁴Molecular Hepatology Section, Department of Medicine II, Medical Faculty Mannheim, Heidelberg University, 68167 Mannheim, Germany</footnote>

<footnote n="5">⁵Lead contact</footnote>

<footnote>*Correspondence: dirk.drasdo@inria.fr
https://doi.org/10.1016/j.isci.2023.108077</footnote>
to construct the DTC. In our nomenclature, a DTC becomes a DT i.e., it is no "candidate" anymore, if it can explain a certain chosen set of "target" observations *in vivo* that are believed the DT needs to be able to explain (a more detailed definition is given in "STAR Methods"). The results of simulations with the DT (or its precursor(s) DTC(s)) do not need to be surprising but define a firm basis, mapping a set of hypotheses on a defined set of readouts. A DT can be used to provide predictions on perturbations and thereby guide experimental strategies. Showing whether any of the hypothesized set of mechanisms is able to explain a certain experimental readout or not, or whether a perturbation experiment is expected to modify the readout, and if it does, in which way, provides important information. Such an information can in most cases due to the complexity of the entire system not with certainty be replaced by schemes, experiments and so forth.

A recent example for a successful DT-model-guided experimentation in liver is ammonia detoxification after drug-induced liver dam-
age.⁹,¹⁰ The set of consensus mechanisms, implemented in a model was unable to quantitatively match the observed blood ammonia con-
centrations and demonstrated that adding a hypothetical ammonia sink mechanism would remove the discrepancy between data and model.
In subsequent experiments, guided by computational models implementing different specific alternative mechanism, an ammonia sink
mechanism could be identified, and gave rise to a potential new therapy approach in hyperammonemia. The mathematical models in the
aforementioned case were based on compartment models not resolving liver micro-architecture but instead assuming one well-mixed ho-
mogeneous reaction volume for each of three (periportal, midzonal, pericentral) liver zones.¹⁰ However, liver micro-architecture may play
an important role in the detoxification of blood from ammonia or other substances such as APAP (acetaminophen, paracetamol; Dichamp
et al.¹¹). Hence a more refined and realistic liver DT should integrate the liver micro-architecture, including the orchestration of events
such as division, death and migration of individual cells for the relevant cell types as well as cell-cell cross-talk, which can well be addressed
by single-cell (agent-) based models.

So far most models addressing liver do not fulfill the requirements of a so-defined DT resolving liver micro-architecture. E.g. Verma et al.¹² mimic the dynamics of molecular factors to regulate the distinct functional states of hepatocytes during regeneration after partial hepatectomy by ODEs (ordinary differential equations), Kuepfer et al.¹³ reproduced the patterns of acute drug-induced toxicity using ODEs, Remien et al.¹⁴ studied liver dysfunction by ODEs, and Naik et al.¹⁵ constructed a multi-level ODE-model, integrating multiple body compartments to investigate hepatic metabolism and its associated de-regulations. Hetherington et al.¹⁶ mimicked the behavior of glucose regulation in liver with diabetes with a PDE (partial differential equation) model, Friedman and Hao¹⁷ explored the efficiency of potential drugs to block the liver fibrosis progress by PDEs, and Schwen et al.¹⁸,¹⁹ studied drug perfusion in the liver using PDEs, too. These differential equation-based models did not address damage or regeneration of the liver at micro-architectural levels, i.e., at cellular resolution, and hence did not use single-cell-based models. Lambers et al.²⁰ studied liver perfusion at the level of lobules representing their statistical shape but not resolving cells individually. Pellicer-Valero et al.²¹ addressed the biomechanical properties of the liver at the entire organ scale in a continuum approach. Some modeling works do consider cells individually, but either in 2D, addressing different questions, focusing on one cell type only, considering a static scaffold structure, or choosing rule-based models for which the representation of a proper biomechanics is questionable. For example, Dutta-Moscato et al.²² used an agent-based model to mimic liver fibrosis formation and integrated different cell types. However, their model is two dimensional, hence quasi-one-dimensional manifolds as capillaries cannot be crossed by a migrating cell without breaking the manifold for topological reasons, while in 3D a cell can in principle easily migrate around a 1D structure. The interactions of the cells in this model were rule based, which may well qualitatively capture certain physical observations, but have a significant risk to miss out the quantitative dynamics. A number of recent works recently studied the consequence of APAP (acetaminophen, paracetamol) overdose on liver, which is known to kill Cytochrom-P450 positive hepatocytes, which are located close to the central vein of each liver lobule. Means and Ho²³ addressed APAP detoxification in a microarchitecture model that was built on a fixed cell scaffold in 2D. Adhyapok et al.²⁴ studied liver damage by overdose of APAP and its subsequent repair using a 2D cellular automaton model on a hexagonal lattice where each cell is modeled as one hexagon. Heldring et al.²⁵ model again APAP toxicity using a 2D multiscale mechanistic model involving various cell types and integrating a drug detoxification pathway in a 2D Cellular Potts model of an individual liver lobule. Dichamp et al.¹¹ established a multicompartment model of APAP toxicity, which embeds a 3D quantitative liver lobule model of APAP toxicity at liver microarchitectural level in an extrahepatic compartment model, which mimics the transport of APAP into, within, and out of the liver lobule capillary network. The model has been used for in vitro-to-in vivo extrapolation and directly confronted with toxicity data on mouse. As in Heldring et al.²⁵ the APAP detoxification pathway is executed in each individual hepatocyte; Sluka et al.²⁶ addressed a liver-centric multiscale modeling framework for xenobiotics in 2D comprising an individual perfused sinusoid with its adjacent hepatocytes. Hepatocytes are modeled agent-based (Cellular Potts model), addressing hepatotoxicity, but not regeneration. Wambaugh and Shah²⁷ studied chemical metabolization in the liver, focusing on hepatocyte function and representing each individual hepatocyte within a rigid vascular scaffold. A former work of some of us addressed 3D liver microarchitecture predicting and subsequently validating an order principle, namely the alignment of hepatocyte division along its neighboring sinusoids (HSA), during regeneration after drug induced damage.²⁸ HSA could also explain tumor phenotypes in drug-induced hepatocellular carcinoma.²⁹ The same model, refined and extended, could quantitatively explain multimodal data in liver regeneration after partial organ removal (partial hepatectomy; PHx) in mouse, representing an entire multilayer slice of a liver lobe of a mouse composed of many lobules, and correctly predicted the proliferation pattern after PHx in pig.³⁰ However, the models in those works were limited to hepatocytes as only cell type, which is not sufficient to address the complex orchestration of different cell types at work during regeneration, degeneration, or disease progression. Many of the experiments performed to understand liver regeneration or disease progression today include the depletion of modification of certain liver cell types beyond hepatocytes, which requires models representing those cell types.

In this article, we create a set of DTs of liver regeneration after acute liver damage by APAP overdose in mouse. Different from the work of Hoehme et al.²⁸ that has focused on hepatocytes only, our model integrates many, but at that stage not all cell types involved, and their
possible orchestration to permit performing realistic *in silico* hypothesis testing to accompany wet-lab experiments. As a proof of concept, the consequence of many alternatively discussed cell-cell cross-talks in 3D microarchitecture is studied, representing each cell type as individual physical entity with a realistic representation of its shape and biomechanics within a single-cell-based model framework. In the first step, a set of DTCs are constructed from literature knowledge and own experiments. Each of the DTCs is constructed to both implement the interaction network between the cell types by molecular signals and physical forces, as well as to be physically realistic. This facilitates parametrization by measurable parameters so that their ranges can be identified. Taken together, the DTCs represented here differ from former work in its main objective as well as in its composition and richness. In the next step, we identify those DTs as that subset of DTCs for which the simulation readout data agrees to a certain set of experimental data and observations that characterizes a normal regeneration process such as the occurrence and size of the dead cell lesion, the spatial distributions of certain cell types in time, or as an additional observable, the concentration of signaling molecules. DTCs for which the simulation results disagree with the experimental target observations are considered to correspond to sets of hypotheses that are either incomplete or false and are therefore ruled out. One of the DTs is then considered as a reference. In the follow-up step, virtual perturbation experiments are run with this specific (reference) DT, such as abrogating a mechanism or predicting the influence of cell type depletion (e.g., neutrophils, HSCs, or macrophages) on liver regeneration. Those perturbations, which are most "informative" in the sense that these perturbations change the result of the regeneration process in the simulation, can be proposed for "wet-lab" approval. Their result can then be used to refine the respective DT and run the next virtual experiment. Such an iterative strategy stepwise improves the liver DT and reduces the number of necessary experiments, thereby speeding up sequences of experiments and reducing costs. Finally, it will result in a digital liver twin that captures many different processes.

Our liver DT is not thought to replace experiments per-se but to reduce their number. The objective is to establish liver DTCs/DTs to test hypotheses with regard to their expected consequence on specific observables in virtual experiments as a pillar complementary to wet-lab experiments.

## RESULTS

### Modeling approach

#### Construction of a lobule: Geometry of lobule, choice of cell types, number and distribution

The chosen tissue unit in this work is the individual hepatic lobule. Following previous work (Figure 1B, Hoehme et al. 28), modeling was performed in a statistically representative lobule of hexagonal shape, with the central vein in the center of the lobule and portal veins in three of its 6 corners. The precise micro-architecture of the lobule was constructed from parameter distributions obtained from 3D liver tissue reconstructions. Upon the injection of a hepatotoxin, such as $CCl_4$ (or APAP), the liver parenchymal cells, namely hepatocytes, in the pericentral region of the lobule die, resulting in a pericentral dead cell area. The starting state of the simulation was defined by labeling pericentral hepatocytes as injured and potential candidates to die. According to the current understanding of hepatotoxicity upon $CCl_4$ or APAP, the injured hepatocytes are those that are cytochrome P450-enzymes positive and received a critical dose of $CCl_4$ (or APAP).31,32

#### Identification of relevant cell types and cell-cell interactions

In a next step, relevant cell types were identified and incorporated in the lobule. These comprise hepatocytes, sinusoidal endothelial cells organizing the sinusoidal network, hepatic stellate cells (HSCs), macrophages (both Kupffer cells and infiltrating macrophages), platelets, and neutrophils (Figure 1C). Two types of interactions were considered: Mechanical interactions by adhesive and repulsive forces as well as mechanical friction forces, and intercellular communication via signaling molecules, including DAMPs (Damage-associated molecular pattern), PDGF (Platelet-derived growth factor), TGFβ (Transforming growth factor β), CXCL1 (C-X-C motif chemokine ligand 1), and CCL2 (C-C motif chemokine ligand 2).

The interaction network of the different cell types was classified into highly probable "cl-1" (with direct experimental support) or probable "cl-2" (with indirect experimental support) (Figure 1C).

There is a huge body of experimental data on specific aspects of the regeneration process that partially looks as alternative or backup mechanisms, which do not significantly modulate the liver regeneration model. We here focus on the plausible selection of those mechanisms emerging from the data that have been implemented in the model. We first detail these mechanism before discussing possible alternatives.

In the model, DAMPs are released by the injured hepatocytes,33 e.g., to activate Kupffer cells (cl-1). This is based on the report that injured hepatocytes can release HMGB1 and HSP-79¹,⁵,³⁴ (both are DAMP signals) to activate Kupffer cells³⁵ (cl-1). Platelets are recruited to the site of injury and bind with sinusoidal endothelial cells during the early stage of liver injury36,37 (cl-1) to promote the proliferation of hepatocytes³⁸ and synthesize PDGF³⁹ (cl-1). PDGF acts on activated HSCs, which display upregulated PDGFR expression to induce proliferation and migration⁴⁰,⁴¹ (cl-1). In the model, the source of PDGF was simplified as the positions of sinusoids inside the dead region, of activated Kupffer cells and of infiltrating macrophages. TGFβ is synthesized by activated Kupffer cells based on the report that TGFβ is predominantly expressed in Kupffer cells⁴² (cl-1). TGFβ acts on HSCs, which are highly responsive to that cytokine, and in most cases, TGFβ is described as prominent driver of activation⁴³,⁴⁴ (cl-1); there are also cases described, where TGFβ is more involved in HSC survival than activation, whereby another factor secreted from infiltrating macrophages is a more prominent driver of activation⁴⁵ (cl-2). In the model, two different factors to activate HSCs were assumed: one is TGFβ that is provided by the activation of ECM deposited Latent TGFβ,⁴⁴ subsequently synthesized by the liver non parenchymal cells, that is, activated Kupffer cells, HSCs, and LSECs. The other factor is an assumed alternate factor generated by activated infiltrating macrophages. Then, CXCL1 is synthesized by both activated Kupffer cells and HSCs to attract neutrophils⁶,⁴⁶ (cl-1). CCL2 is
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Digital Twin (Candidate) of multi-cellular lobule system based on experiments</figcaption></figure>

(A) APAP-induced liver regeneration on mouse and the spatiotemporal pattern of the different liver cell types. We have corresponding experimental stainings as reference.

(B) A lobule in hexagonal shape is constructed to study the liver regeneration process upon toxin-induced acute damage following our previous work.²⁸ In this study, $CCl_4$ is used to induce lobular damage, where the necrotic lesion size reaches the maximum area 2 days after the injection of $CCl_4$.

(C) The network includes relevant cell types and intercellular signals contributing to liver regeneration. “Positive” means to activate or to attract a certain type of cell; “Negative” means to kill, to eliminate, or to deactivate a certain type of cell. The number in the yellow circle marks the related reference (Table S2). The number in the cyan square indicates its classification, (cl-”)1” indicates the relationship has literature support; (cl-”)2” indicates the relationship is discussed as option by literature without data-based evidence.

(D) For the different liver cell types, specific geometric objects are used to represent their shape (gray: hepatocytes; red: sinusoids; cyan: HSCs; brown: macrophages; green: neutrophils). Cell movement is updated by solving a force-velocity equation according to the respective cell-mechanical properties. The signal gradient is updated by solving a related reaction-diffusion equation.

(E) Illustration of different cell types at different days after the injury.

synthesized by activated HSCs, Kupffer cells, infiltrating macrophages and sinusoids in the lesion region to attract infiltrating macrophages (cl-1). This is based on the report that activated stellate cells, Kupffer cells, macrophages, and endothelial cells secrete CCL2 to control the macrophage infiltration.⁴⁷,⁴⁸

In the model, a subpopulation of Kupffer cells are assumed to be able to eliminate dead hepatocytes based on the observation that dead hepatocytes are engulfed by Kupffer cells⁴⁹ (cl-1). Moreover, infiltrating macrophages are assumed to also eliminate dead hepatocytes as previously described by Boulter et al.⁵⁰ (cl-2).

We have three additional hypotheses on the function of infiltrating macrophages. The first is to contribute to the activation of HSCs e.g., based on reports that infiltrating macrophages secrete at least one HSC activating factor⁴⁵ (cl-2); the second is to induce death of the activated HSCs based on the report that infiltrating macrophages might trigger apoptosis of activated HSCs⁸ (cl-2); the third is to revert the activated HSCs to the quiescent phenotype⁵¹⁻⁵³ (cl-2). This is probably depending on different macrophage phenotypes/subpopulations, as becomes more and more evident from recent scRNASeq experiments.⁵⁴
<figure><img src="image_5.png" /></figure>

One advantage of in silico modeling is to permit straightforward testing of all different selected regeneration scenarios, which may then serve as a guide to exclude those ones for experimental validation that by the simulations have to be assumed *apriori* to fail. This is demonstrated later by starting from a reference model that has been demonstrated to reproduce the experimentally observed regeneration scenario (Figure 2) and then running simulations where interactions labeled as “cl-2” are switched on or off, therewith testing the influence of each of the following hypothesized interactions (and the combination of some of them) on the process of liver regeneration by implementing each time one of them (Figure 3).

## Definition of a reference model (digital twin candidate)

We define a reference DTC by the hypotheses carrying a "cl-1" in Figure 1C and the following hypotheses (Hypoth.) carrying a "cl-2" in Figure 1C (the 'X' in 'no. X' refers to the reaction identifier in Figure 1C).

(1) **Hypoth. 1.** Infiltrating macrophages eliminate dead hepatocytes (Figure 1C, no. 22, cl-2);

(2) **Hypoth. 2.** DAMPs activate infiltrating macrophages (Figure 1C, no. 3, cl-2);

(3) **Hypoth. 3.** HSCs are activated by TGFβ, e.g., from Kupffer cells (no. 21, cl-2);

(4) **Hypoth. 4.** HSCs attract Kupffer cells to migrate (no. 24, cl-2);

(5) **Hypoth. 5.** Infiltrating macrophages revert activated HSCs to a quiescent phenotype (no. 19, cl-2).

If these hypothesized interactions (Hypoth 1–5) were all implemented, regeneration succeeded in the model (Figure 2), hence the respective DTC qualified by its agreement with the experimental data as DT.

We compare the spatial-temporal dynamics of the reference DT with five different models by switching each time one of the above hypotheses to one of the following alternative hypotheses (AltHypoth.) i.e., Hypoth X - > AltHypoth X, whereby all other hypotheses ($Y \neq X$) do not change (this switch is straightforward to generate one DTC from another):

(1) **AltHypoth. 1.** Infiltrating macrophages do not eliminate dead hepatocytes (Figure 1C, no. 22, cl-2);

(2) **AltHypoth. 2.** DAMPs do not activate infiltrating macrophages (Figure 1C, no. 3, cl-2);

(3) **AltHypoth. 3.** HSCs are activated by an alternate factor produced by infiltrating macrophages (no. 18, cl-2);

(4) **AltHypoth. 4.** HSCs do not attract Kupffer cells to migrate (no. 24, cl-2);

(5) **AltHypoth. 5.** Infiltrating macrophages induce cell death of activated HSCs (included in no. 19, cl-2).

Hence the liver DT approach permits to study in how far alternative assumptions on cell-cell interactions impact on the regeneration of the peri-central drug-induced lesion.

## Timing of events during drug-induced liver injury and regeneration as implemented in the reference digital twin candidate

Right after the injection of the hepatotoxic substance (CCl₄/APAP), hepatocyte injury is expected to occur in less than an hour, supported by fact that in mice, for example, the highest APAP blood concentration ("Cmax") is reached ~30 min after injecting the drug into the peritoneum.<sup>31,55,56</sup> Our recent *in vitro* data has shown that there exists direct effect of APAP on Cyp450 positive hepatocytes via NAPQI generation.<sup>11</sup> The timing of the interaction processes not directly depending on NAPQI is explained in Figure S2 as following: (1: 1h) Up to about 1 h after drug administration, injured hepatocytes secrete DAMPs<sup>35</sup> (Figure S2A) that activate Kupffer cells (Figure S2B), making them secreting CXCL1, TGFβ and CCL2 (Figure S2C). At the same time, sinusoidal endothelial cells/platelets localized in the lesion are secreting PDGF and also CCL2 (Figure S2A). PDGF attracts HSCs (Figure S2B), which then causes the migration of Kupffer cells toward them (Figure S2C). (2: 2h) CXCL1 attracts neutrophils, which about 2 h after drug administration are observed to migrate toward the lesion<sup>57,58</sup> (Figure S2D). Neutrophils in contact with injured hepatocytes are assumed to initiate their death (Figure S2D). (3: 1day) The dead hepatocytes are eliminated by macrophages after about 1 day: Kupffer cells are reported to eliminate the bodies of dead hepatocytes<sup>49</sup> (Figure S2E). The population size of Kupffer cells dropped after activation until day 2 and recover thereafter as observed by Dragomir et al.,<sup>59</sup> Graubardt et al.<sup>60</sup> and in our experimental data (Figure 1A). About the same time, TGFβ (experimentally hypothesized to be in a first step derived from LTGFβ deposits in the ECM adjacent to the damaged hepatocyte area, and activated by the dying hepatocytes<sup>44</sup>) activates HSCs making them secret CCL2 and CXCL1 amplifying and/or backing up the effect of Kupffer cells and platelets (Figure S2E). Activated HSCs proliferate and secrete extracellular matrix (ECM).<sup>2</sup> At the same time, infiltrating macrophages are attracted by CCL2 (Figure S2F) approaching the dead hepatocytes (Figure S2G). (4: about 1.5days) Around day 1.5–2, hepatocytes outside of the dead lesion start to enter S-phase to eventually replace the dead hepatocytes.<sup>28</sup> Sinusoids located in the lesion are among the prime candidates for HGF and EGF secretion, which are both mitogens for hepatocytes<sup>61,62</sup> (Figure S2H). The infiltrating macrophages are initiated as monocytes with Ly6C-high phenotype. (5: 2-3days) After a period of time between 2 and 3 days, the infiltrated monocytes transform into macrophages with Ly6C-low phenotype.<sup>60,63</sup> They can phagocytose the dead hepatocytes when they are at Ly6C-high phenotype. Once they adopted the Ly6C-low phenotype, they can phagocytose or revert activated HSCs to quiescence.<sup>8,64,65</sup> NK (natural killer) and NKT (natural killer T) cells also contribute to induce death of activated HSCs<sup>66</sup> (Figure S2I).

Summarizing, the state changes of cell types concern attributes such as activation and deactivation, or initialization for differentiation, as we now know into multiple different functional phenotypes, for example resident and infiltrating macrophages, and among the latter, the Ly6C-high and Ly6C-low phenotype. The cellular consequences can be, among others, migration, reversion from activation, killing, proliferation and death, depending on the specific cell type (specified in Table S2), and other parameters.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Pattern of liver lobule regeneration (reference model)</figcaption></figure>

(A) Simulated lobule, consisting of different cell types and the sinusoids over time, taking into account different simulation scenarios. Gray spheres are healthy hepatocytes, black spheres are dead hepatocytes; sinusoids are visualized as red lines. The software ParaView is used for visualization.

(B) The distribution of HSCs (dark red: activated; gray: quiescent) and macrophages (brown) over time. The distribution of Ly6C-high monocytes on days 1 and 2 are shown in the blue boxes. Their precise spatial pattern depends on the timepoint of secretion and the range of CCL2.

(C and D) Lesion area and number of hepatocytes over time from both, simulations and experiments.

(E) Spatial distribution of activated HSCs, macrophages, and neutrophils in the lobule, measured by counting the corresponding cell number over the distance to the central vein (CV). As illustrated, there are six regions considered: Region I (<21.4 µm to the CV), II (between 21.4 and 42.8 µm to the CV), III (between 42.8 and 64.2 µm to the CV), IV (between 64.2 and 85.6 µm to the CV), V (between 85.6 and 107 µm to the CV), VI (between 107 and 128.4 µm to the CV). (F to H) Number of activated HSCs, macrophages, and neutrophils in all regions over time. (In C, the points specify the mean, the bar the standard deviation. In D, the data has been represented in a box plot displaying the median (black bar), the first (Q₁) and third (Q₃) quantile (delineating the grey shadded area), and wiskers representing Q₁−1.5 * IQR and Q₃+1.5 * IQR, respectively, whereby IQR denotes the interquatile range IQR = Q₃−Q₁. In F–H, the bold lines represent averages, the error bars the standard deviation over four simulation runs, using different random seeds.)
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Pattern of liver lobule regeneration resulting from alternative hypothesized interactions</figcaption></figure>

(A) The regenerating lobule over time under reference and five alternative hypothesized interactions: (1) Infiltrating macrophages do not eliminate dead hepatocytes; (2) DAMPs do not activate infiltrating macrophages; (3) HSCs are not activated by TGFβ but by one factor from infiltrating macrophages; (4) Kupffer cells do not migrate toward HSCs; (5) Kupffer cells and infiltrating macrophages eliminate activated HSCs instead of reverting them to quiescent mode.
(B) The distribution of activated HSCs (antibody: αSMA) and macrophages (antibody: F4/80) over time under reference and five perturbed interactions.
(C and D) Lesion area and hepatocyte density over time under reference and five alternative hypothesized interactions. (The data representation is the same as in Figures 2C, 2D, and 2F-2H.)
## Digital twin candidate at cell and tissue level: cell geometry

Hepatocytes, hepatic stellate cells, Kupffer cells and infiltrating macrophages, as well as neutrophils are modeled, each as individual entity. The endothelial cells are not modeled individually, but as part of a network of sinusoids. Each hepatocyte has been approximated by an isotropic, elastic, and adhesive sphere (Figure 1D, gray objects), named "center-based model" (CBM), capable of interacting with other cells or blood vessels by mechanical forces or chemical signals. The sphere can be thought as specifying the region in space where the hepatocyte is localized with overwhelming probability. The CBM has been parameterized by material and cell-kinetic parameters, which permits to readily identify the physiological parameter ranges. Macrophages and neutrophils were equivalently mimicked by a CBM, but with different cell parameters (Figure 1D, green and brown objects). Different from the former cell types, the HSCs were approximated by an isotropic, elastic, and adhesive sphere with chains of elastic springs emanating from their body to capture their long protrusions (Figure 1D, cyan objects). The sinusoidal network is modeled as semi-flexible chains of spheres that are connected by springs (Figure 1D, red objects). This design accounts for the fact that vessels resist bending and stretching. The parameters of the sinusoids comprise volume, density and branching orientation, which are sampled from the scanning data of the real liver sinusoidal system (see more details in Hoehme et al.²⁸).

## Digital twin candidate at cell and tissue level: Force balance and cell movement

Cell and sinusoid movement is computed based on force balance. The DTC takes into account passive and active forces. Passive forces are friction forces of each cell with its environment (other cells, sinusoids, intercellular medium), deformation and compression forces experienced by a cell, as well as cell-cell and cell-sinusoidal adhesion forces. Cell migration, as it occurs by the anchoring of cells in the extracellular matrix, e.g., in the space of Disse, is mimicked as an active force. The precise form of the forces (as detailed in the STAR Methods) has been chosen to directly represent cell material parameters. For example, the cell-cell interaction force has been approximated by the "Johnson-Kendall-Roberts" (JKR)- force model for homogeneous elastic sticky spheres, which has been shown by micro-pipette experiments⁶⁷ to quantitatively reproduce the force-distance relation of two cells brought in contact and pulled apart. In conjunction with the friction force, the emerging behavior at the tissue level is viscoelastic. During liver regeneration, cell proliferation causes cell compression, which cannot be properly addressed by standard pairwise forces (as JKR, Hertz and so forth).⁶⁸ To correct for this shortcoming, the JKR-force was modified for small cell-cell distances by a term that accounts for volume compression forces arising from large cell deformation, calibrating with the "Deformable Cell Model (DCM)."⁶⁹

The motion of cells and sinusoids is updated by solving an overdamped stochastic equation of motion, a Langevin equation, which summarizes all forces exerted on them: $\sum F_{fri} + F_{def} + F_{adh} = 0$, where $F_{fri}$, $F_{def}$, $F_{adh}$ are friction force with the environment, deformation force, and adhesive force with other cells or elements, respectively (Figure 1D, see more details of each force term in STAR Methods).

## Digital twin candidate at cell and tissue level: Molecular interactions

Cell types may secrete signal molecules that are sensed by other cells. The spread of the molecules is mimicked by a diffusion equation. The cells generating the molecules are represented as source terms in the diffusion equation, turning the diffusion equation into a reaction-diffusion equation. Moreover, the equation contains a general first-order kinetics decay term (molecule degradation term in Figure 1D, see more details of each equation term in STAR Methods).

## Digital twin candidate parameterization and initial conditions

The model parameters characterizing the biomechanical properties of hepatocytes and sinusoids, and the liver lobule micro-architecture have been chosen as previously developed²⁸ (Table S1). The densities of non-parenchymal cells not considered in that reference has been estimated from other published references, as for neutrophils for example, from McDonald et al.,⁵⁸ Marques et al.⁷⁰ (Table S1).

In addition, for hepatocytes, the direct effect of APAP on Cyp450 positive hepatocytes via NAPQI generation was included. This was realized by constructing a dose-dependent cell death induction in this subpopulation of hepatocytes, and the time (point) at which death occurs. Both parameters are based on experimental data in that the NAPQI-pathway displays a time and dose-dependent (in vitro: concentration dependent).¹¹

Both HSCs and Kupffer cells are initially (at the time point of drug administration) distributed homogeneously in the lobule (Figure 2B, quiescent/activated HSCs are colored in gray/red and Kupffer cells are colored in brown). As explained above, the population size of Kupffer cells has been modeled to drop after activation until day 2 and recover thereafter. Death of Kupffer cells has been mimicked by a death rate such that the decay in the F4/80 (macrophage marker) staining data in Figure 1A, could be qualitatively reproduced. (We also considered a constant population of Kupffer cells to test the sensitivity of the regeneration scenario with regard to the Kupffer cell population.) The recovery has been modeled by adding new Kupffer cells with a certain rate. The new Kupffer cells can, as the Ly6C-low phenotype, phagocytose the dead hepatocytes. The infiltrating macrophages are initiated as monocytes with Ly6C-high phenotype, which do not count as the macrophage population. Those monocytes that are in the range of CCL2 (i.e., where the concentration of CCL2 is higher than a threshold value specified in Table S1) move during the two days toward the hepatocyte damage region. After a period of time between 2 and 3 days, the infiltrated monocytes transform into macrophages with Ly6C-low phenotype. They can phagocytose the dead hepatocytes when they have adopted the Ly6C-high phenotype. Once they adopted the Ly6C-low phenotype, they can phagocytose or revert activated HSCs.
<figure><img src="image_9.png" /></figure>

Neutrophils in contact with injured hepatocytes are assumed to initiate their death reaching a state in which they become "flagged" for a potential elimination by macrophages after about 1 day.

For molecular signals, the diffusion rates of all molecular signals are taken from literature (Table S1).

## Simulated results of regeneration scenarios with the reference and alternative interaction models

The timing of the processes depicted in Figure 1C follows the assumed logical order of those subprocesses defining the reference DTC (later in discussion shown do qualify as a DT) and their possible alternative DTCs described above (emerging from Figure S2). So, our set of liver DTCs can reproduce each step of the processes at each time point (e.g., after about 1 h, Kupffer cells are activated, after 3 days, infiltrating macrophages switch from Ly6C-high to Ly6C-low phenotype), and permits to study the multicellular organization process emerging from the complex interactions in time and space.

### Simulation result of the reference model

After 1 day, the lesion is filled by dead hepatocytes (black in Figure 2A, day 1). HSCs migrate and accumulate in the lesion and the number of Kupffer cells is decreased (Figure 2B, day 1). After 2 days, most HSCs in the lesion are activated and the monocytes with Ly6C-high phenotype are accumulated in the lesion (Figure 2B, day 2). This leads to a clearance of the lesion from dead hepatocytes (Figure 2A, day 2). After about 3 days, macrophages begin to deactivate HSCs (Figure 2B, HSC, day 4 & 6). At day 6, the lesion is closed by healthy hepatocytes having re-placed the removed (dead) hepatocytes (Figure 2A, day 6). In different simulation runs, we find sporadic dead hepatocytes that are later phagocytosed, while in some simulations the dead hepatocytes have been entirely eliminated by day 6 and the lesion was closed (Figure 2A, day 6). Infiltrating macrophages have already disappeared. The HSCs become deactivated, and both HSCs and Kupffer cells redistribute in the lobule (Figure 2B, day 6), the latter accompanied and driven by a change of fate from Ly6C-high phenotype to Ly6C-low phenotype (see details in STAR Methods).

For the lesion area and hepatocyte density over time, for which we had quantitative experimental values,²⁸ the simulation results show a perfect agreement to the experimentally observed values (Figures 2C and 2D). In Hoehme et al.,²⁸ the lesion area was experimentally defined as the area not containing hepatocyte nuclei anymore, while signs of cell death occurred already earlier, but were not considered in the experimental curve (SI to that reference). Referring to the direct effect of NAPQI, the lesion may alternatively be assessed regarding the lobule space not occupied by "healthy" hepatocytes, in which case the dead cell area shows up earlier as in Figure 2C, different from the previous model that did not take into account the direct effect by NAPQI-detoxification.²⁸ Following the same line of argument, the hepatocyte density may be defined as the number of hepatocytes over the area of the lobule (as in Hoehme et al.²⁸) or the number of healthy hepatocytes over the area of the lobule. Both measures result in different spatial profiles before 2 days after the injury. In Figure 2D we depicted the number of healthy hepatocytes, which may be difficult to demarcate from those already damaged and about to die in the experiment.

We lacked quantitative values for the other cell types, so estimated their density from published references (Bouwens et al.,⁷¹ Wake,⁷² McDonald et al.,⁵⁸ Zigmond et al.⁶³) and measured in the simulation the change of their profile, which may be seen as a mainly qualitative model prediction. In order to simulate the dynamic change of the distributions of HSCs, macrophages (both Kupffer cells and infiltrating macrophages), and neutrophils, the lesion region was divided evenly into 6 sub-regions according to their distance to the CV (Figure 2E), and the number of each type of sinusoidal cell type in each sub-region was counted over time. As shown in Figure 2F (experimental: Figure S1D), the number of activated HSCs starts to rise after day 1. This parallels with the time point when Kupffer cells begin to eliminate the bodies of dead hepatocytes. The number of activated HSCs is highest in region I, which is located most close to the CV. The number of activated HSCs peaks around day 2, then dropping dramatically due to the interventional effect of Ly6C-low expressing infiltrating macrophages. At around day 6, there are no activated HSCs left in the lobule. As shown in Figure 2G, the number of macrophages in the pericentral area decreases before day 2, due to toxic damage upon phagocytosing the dead hepatocytes. After day 3, the Kupffer cell number quickly rises again, due to the repopulation of the pericentral area by the proliferation and migration of the resident Kupffer cells and/or differentiation from the massively infiltrating macrophages into the lesion. After day 4, the macrophage population in the lesion drops again, this time due to the depletion of infiltrating macrophages. After around day 6, the number almost resumes to initial values, e.g., resulting from the Kupffer cells relocating back to their initial distribution. As shown in Figure 2H, the number of neutrophils rises quickly after day 0, attracted by a CXCL1 gradient and in order to deplete the injured hepatocytes. Recruited neutrophils accumulate in region I and the number rise until day 1. Until day 3, the neutrophil population disappears.

Our simulation of the reference case can achieve the expected pattern of normal liver regeneration and agree with the experimentally observed dynamics for the necrotic lesion size and hepatocyte density. Next, we explored the communication between different cell types and factors to achieve liver regeneration after NAPQI induced damage by testing type "2" interactions (i.e. those "probably" happening in vivo) in the network.

### Constant Kupffer-cell population

As the degree of staining may not directly and linearly relate to the Kupffer cell population size, we studied the impact of a variation of the Kupffer cell population size by considering the simpler case of a constant population size and could not find any effect (see details in SI). Alternatively, one may model the Kupffer population size kinetics as an appearance (birth) and death process (death at days 1–2, appearance and spread from day 2 on), but this was not the purpose of this work and no changes in the regeneration scenarios were expected.
Alternative hypothesized cell-cell interaction network models (AltHypoth)

In a next step variations of the reference DT ("R-DT") were considered by individually modifying the following five interactions of the reference network: (1) Infiltrating macrophages cannot (R-DT: can; no. 22) eliminate the dead hepatocytes; (2) Infiltrating macrophages are not activated by DAMPs (R-DT: are activated; no. 3); (3) HSCs are activated by another factor, besides TGFβ, which is produced by activated macrophages derived from infiltrating monocytes (R-DT: TGFβ is activating HSC, no. 21; we assumed that this factor has about the same diffusion constant as TGFβ); (4) HSCs cannot attract Kupffer cells to migrate (R-DT: they do, no. 24); (5a) Infiltrating macrophages deplete activated HSCs (Figures S6A–S6E); (5b) 50% of activated HSCs are depleted/50% are reverted (Figure S6G) (R-DT: infiltrating macrophages revert activated HSCs to quiescent phenotype, no. 19).

Later in discussion, we refer to a DTC, in which the assumption (k; k = 1, 2, 3, 4, 5) has been modified partially as "model". As shown in Figures 3A, 3C, and 3D, the alternative models (2) and (5) do not significantly differ from the result of the R-DT regarding a complete liver regeneration, while the modifications (1), (3), and (4) influence the lesion recovery in that regeneration after 16 days is still incomplete. Figure 3A displays a single simulation scenario, where dead hepatocytes are colored in black. Figure 3C shows the average over 4 runs with different random seeds. Among the alternative DTCs (1), (3) and (4), the area of the unhealed lesion is largest for DTC (1), due to the inhibition of phagocytosis of infiltrating macrophages. The Kupffer cells alone are not able to clear dead hepatocytes in time. However, a parameter sensitivity analysis suggests that Kupffer cells equipped with a stronger phagocytosis capacity, by a shorter elimination duration, would alone be able to clear the necrotic area from dead hepatocytes in the experimentally observed time. i.e., by model-guided parameter sensitivity analysis the phagocytosis capacity of the resident macrophages is identified to be a critical sensitive parameter required for a quantitative understanding of the regeneration process (Figure S6H; matching with the R-DT). The importance of infiltrating macrophages in clearing up the lesion is also indicated by the modification implemented as DTC (3), where HSC activation is depleted (Figure 3B; Figure S4D). Due to the lack of activated HSCs in the lesion, the expression level of CCL2 is not sufficient to attract infiltrating macrophages into the lesion to clear up dead hepatocytes (Figure 3B; Figure S4C; Figure S5B). Similar to DTC (3), in DTC (4) the number of Kupffer cells in the lesion is much smaller than that in the reference case, hence less TGFβ is produced in the lesion resulting in decreased numbers of activated HSCs (Figure 3B; Figures S4C and S4D). Again, the expression level of CCL2 is reduced, and as a consequence, fewer infiltrating macrophages are attracted into the lesion, as compared to the reference case (Figure 3A; Figure S4C).

We further tested a DTC (6: 2&5) simultaneously implementing assumptions (2) and (5), which did not have significant impact on liver regeneration (Figures S6F–S6H).

Our computational study of these five individually modified interactions suggests two main conclusions: (1) There is a positive feedback loop between Kupffer cells and HSCs, where the presence of HSCs promotes the migration of Kupffer cells into the lesion, while increasing numbers of Kupffer cells in the lesion increase the number of activated HSCs, which attract more infiltrating macrophages to help clear up dead hepatocytes; (2) The way to diminish activated HSCs by macrophages, either by engulfment or by phenotype reversion has no significant impact on lesion recovery. Note however, that in case HSCs are depleted, the remaining HSCs would have to proliferate to repopulate the lobule to its original population size. It is now obvious from experimental data that half of the population of activated HSCs is depleted and half of it is reverted, which by itself would more or less lead to the repopulation of the liver to the original healthy state².

With regard to the digital twin classification, besides the R-DT, DTCs (2), (5), (6) qualify as DTs as for each of them, the simulation results agree with the experimental data.

In a next step, we test how far simulations of perturbation scenarios are suited to guide experiments in a way to permit pinpointing differences between interaction mechanisms at the tissue level, and hence serve to validate the model predictions.

Simulated perturbation experiments: Depletion of non-parenchyma cells

After the simulated testing of the reference digital twin (R-DT; all mechanism given by cl-1 in Figure 1C and Hypoth. 1–5) and of variations that differ from the R-DT by one interaction (all mechanisms given by cl-1 in Figure 1C and the replacement of one Hypoth. X by AltHypoth. X; X either 1, 2, 3, 4, 5), additional perturbation simulations were performed to predict the impact of depletion of cells typically found along the sinusoidal network spaces on liver regeneration. In each of four perturbation simulations (I-IV), one of the four cell types HSCs, Kupffer cells, infiltrating macrophages and neutrophils was depleted individually.

A further advantage of the DT is that it permits to test hypothetical cases and see how much contribution to an observed effect can be attributed to sub-processes. Both will be discussed later in discussion.

Perturbation I: The digital twin predicts that the depletion of hepatic stellate cells results in an unhealed lesion

The depletion of HSCs (Perturb. 1; Figure 4A, scenario 1) did not change the size of the lesion induced by the drug and characterized by dead or dying hepatocytes. This turns out to be the same for all cell type depletions, indicating that the lesion size is mainly controlled by the cell death caused downstream of the NAPQI pathway. In the case of HSC depletion, however, the dead hepatocytes are not cleared up. Due to the lack of activated HSCs, CCL2 levels are not sufficient to attract infiltrating macrophages to clean up dead hepatocytes. In addition, due to the absence of HSCs, Kupffer cells do not migrate. Consequently, the number of Kupffer cells in the lesion is smaller than in the reference case. This is similar to the perturbation that disables HSCs to attract Kupffer cells to migrate (scenario 4, Figure 3A). In general, the depletion of HSCs leads to reduced numbers of neutrophils, Kupffer cells, and infiltrating macrophages in the lesion. As a consequence, a large number of dead hepatocytes remains uncleared in the lobule.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Pattern of liver lobule regeneration upon depletion of different cell types in the reference digital twins (Figure 2)</figcaption></figure>

(A) The regenerating lobule over time as reference state and upon depletion of four sinusoidal cell types: (1) Depletion of HSCs; (2) Depletion of Kupffer cells; (3) Depletion of infiltrating macrophages; (4) Depletion of neutrophils.

(B and C) Lesion area and hepatocyte density over time in the reference state and upon depletion of different cell types as indicated. (The data representation is the same as in Figures 2C, 2D, and 2F-2H.)

Our prediction is consistent with a previous study, where the depletion of HSCs resulted in decreased expression levels of CXCL1 and reduced numbers of infiltrating neutrophils.⁷³

**Perturbation II/III: The digital twin predicts that the depletion of macrophages abrogates clearance of lesion from dead hepatocytes, independent on whether it concerns the Kupffer cells or the infiltrating macrophages**

As shown in Figure 4A (scenario 2), after the depletion of Kupffer cells (Perturb. 2), the necrotic lesion is only marginally cleared up from dead hepatocytes. Due to the lack of Kupffer cells, there is a reduced TGFβ signal to fully activate HSCs, leading to qualitatively the same scenario as above after the depletion of HSCs. As Kupffer cells are totally missing, the phagocytosis of dead hepatocytes by Kupffer cells does not take place hence the lesion is not cleared from dead hepatocytes. As now in addition to infiltrating macrophages (cf. scenario 1, these are attracted by activated HSCs) also Kupffer cells are lacking, the lesion is even slightly larger than in scenario 1. We further tested if HSCs can be activated
even after the depletion of Kupffer cells, we found that the lesion can be recovered in time (Fig. S6G). These results indicated the importance of activation of HSCs in regulating lesion recovery.

After the depletion of infiltrating macrophages (Perturb. 3, scenario 3 in Figure 4A), again a large part of the dead hepatocytes remains in the lobule not being cleared. This is similar to the perturbation of disabling infiltrating macrophages to eliminate dead hepatocytes (Figure 3A, scenario 1).

Our prediction is consistent with a previous study of depleting infiltrating macrophages in acute liver injury, which resulted in delayed tissue recovery.⁷⁴

Because of the (expected) dramatic effect of macrophage depletion on the regeneration result, we refined the assumed perturbation by replacing the dynamic kinetics of Kupffer cell numbers in the course of the regeneration process by assuming a constant population size after the 6-h time point.

For this case (Figure S8), we found no consequence on the readout parameters at day 6, interestingly indicating that the precise dynamics of the Kupffer cell population may not be a critical determinant of the regeneration process, as long as Kupffer cells are present for being activated by damaged hepatocytes through DAMPs, and to facilitate HSC activation early (about 1 h) after drug administration.

## Perturbation IV: The digital twin predicts that the depletion of neutrophils has no significant influence on liver regeneration

As shown in Figure 4A, after simulation the depletion of neutrophils (Perturb. 4), there is neither a significant effect on the generation of damage, nor on the regeneration of the necrotic lesion in comparison to the reference scenario. Therefore, neutrophils might be considered as providing backup mechanisms to induce damage toward epithelial cells, to phagocytose the damaged debris, and to induce cell death in HSC.² This is in-line with Krenkel et al.,⁴⁸ stating that cell death induced by neutrophils in APAP-induced liver injury is not very obvious in vivo, even though it can be observed in vitro.⁷⁵

This still needs further debate, as an older study of depleting neutrophils during acute liver injury reported protection of the lobule against hepatotoxicity with a significantly reduced centrilobular necrosis area.⁷⁶ As the latter was not the case in our simulations for the hypotheses studied in the main body of this work, it prompted us to investigate a further highly speculative case, that may not be expected in acute liver damage, but perhaps could be obtained by proper experimental manipulations (Figure S7). Here, we consider regeneration from a necrotic lesion that has not been generated by cell-death due to the NAPQI-pathway. The NAPQI pathway is for example not activated in hepatocytes that lack Cytochrom-P-450 enzymes, as occurring during periportalization after the repetitive administration of APAP (or CCl₄; Ghallab et al.⁷⁷), and another hepatotoxic event would be responsible for the multicellular necrotic lesion. In that simulation, for simplicity, we assumed also a circular necrotic lesion that initiates phagocytosis of dead hepatocyte debris by neutrophils and macrophages. Here, we could observe a protective effect against hepatotoxicity upon neutrophil depletion (Figure S7). The simulations demonstrate that our framework may permit in silico tests of hypotheses that may not be present as such in nature, but may be realized by proper engineered manipulations.

## DISCUSSION

In our study, we developed a set of multi-level computational liver lobule models each integrating a slightly different inter-cellular network of major hepatic cells and takes into account cell-cell communication based on a significant selection of published knowledge. By its construction, each model defines a digital twin candidate (DTC). We applied each DTC to simulate the liver regeneration process after APAP-induced acute hepatotoxicity, thereby testing documented and hypothesized interactions between cells and signals. Furthermore, we made testable predictions on perturbations that were already performed in previous experimental studies.

The starting model of our digital twin (DT) is as simple as possible and implements those basic mechanisms that can be verified (not validated) to yield a successful regeneration in the regeneration time as in the experiment. The second step is knocking down mechanisms in the DT and studying its consequences, giving a prediction, which may turn out to be not correct, followed ideally by validation experiments. If the validation experiments show a deviation from the regeneration process in absence of the knock-down, for example a delay, this would indicate the action of further mechanisms that then need to be integrated to arrive at a new, refined DT.

Our liver DT illustrates that testing or perturbing the inter-cellular communication between cells and signals allows exploring the mechanism behind liver regeneration and guiding the design of relevant experiments. The in silico liver considered here includes reported and hypothesized interactions between cell types, mediated both mechanically and by chemical signals. It permits to study alternative hypotheses and their consequences on the regeneration process. A plausible DTC implemented the following temporal events: (1) Injured hepatocytes produce DAMPs and are killed via cell stress from APAP detoxification by NAPQI (Figure S2). Platelets in the lesion produce PDGF; (2) Moreover, Kupffer cells and HSCs are also activated by DAMPs. The so activated HSCs are assumed in agreement with our data to be not detected through αSMA, and migrate up the gradient of PDGF; (3) activated non-parenchymal cells (and damaged hepatocytes) produce CXCL1, while neutrophils migrate up the gradient of CXCL1 and induce death in injured hepatocytes; (4) Kupffer cells migrate toward the cluster of HSCs; (5) Activated Kupffer cells produce TGFβ to activate another function of HSCs, namely the production of CCL2 to attract infiltrating macrophages into the lesion. We assume that those HSCs are labeled by αSMA; (6) infiltrating macrophages along with Kupffer cells engulf and eliminate the bodies of dead hepatocytes; (7) the healthy hepatocytes surrounding the lesion proliferate and collectively migrate toward the lesion to recover the lost hepatocytes, following a tug-of-war mechanism²⁸,³⁰,⁷⁸; (8) Infiltrating macrophages switch from Ly6C-high phenotype to Ly6C-low phenotype to revert the activated HSCs into quiescent HSCs, whereas activated Kupffer cells switch back to quiescent Kupffer cells; (9) To restore the original cell distribution, we simulated alternative scenarios. Either HSCs previously localized in the lesion may die and be phagocytosed by macrophages, and those outside the lesion could proliferate. Another, perhaps more hypothetical scenario was the
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5.** Summary of the construction and application of the liver digital twin to study liver regeneration after acute damage</figcaption></figure>

migration of quiescent Kupffer cells and HSCs. In a computer simulation comparing both scenarios (not shown) the regeneration of the necrotic hepatocyte lesion turned out to be insensitive to the precise mechanism of restorage of the original cell distribution. During the entire process, different types of hepatic cells interact and collaborate to achieve a perfect liver regeneration in time.

An important advantage of this liver digital twin is that it permits to investigate the role of individual cell types and different signals by simulations. For example, the presence of one or several factors, or the normal or abnormal function of these factors, (e.g., the normal function of TGFβ, or the absence of normal function as in Figure 3, scenario (3)) as well as its interaction with other factors can either be depleted completely, weakened or amplified, and the consequence of this modification on the regeneration process and its final outcome can be studied (in Figure 3, scenario (3) leading to a lack or large absence in lesion closure). This may be impossible in an in vivo experimental model due to several reasons. Firstly, some modifications are technically not feasible or accompanied with huge effort. Secondly, possible backup mechanisms may hide the effect of a given change such that the direct and indirect consequences are not detangled (e.g., a wet-lab validation experiment of Figure 4, scenario (1), suppressing HSCs may not lead to a failure in a closure of the lesion as a consequence of an unrecognized backup-mechanism that so far is not part of the DT). Moreover, DT simulations may show that certain manipulations are likely not resulting in a significant change of outcome, and therefore are not sufficiently informative to justify the effort of experimental realization (e.g., depletion of neutrophils, Figure 4, scenario 4), especially in regard to the 3R for animal experimentation. On the other hand, differences in the outcome of experiments and DT simulations may indicate so far unrecognized mechanisms that would remain hidden without such simulations of a reference situation (e.g., if suppression of HSCs as in Figure 4, scenario (1) does not result in a failure or regeneration, there may be backup mechanisms that may then be added to the existing DT as hypothetical mechanisms). Hence, such a liver DT may help to systematically classify the mechanisms, hypothetical or not, with regard to their importance in closing the dead cell lesion.

To demonstrate opportunity and power provided by the digital liver, we simulated a number of perturbation experiments in regeneration after drug-induced pericentral damage (two of them indicated above already) and classify their result as one of three possible basic scenarios, as follows (Figure 5).

## Incomplete regeneration (Figure 5, blue scenario)

An incomplete regeneration with dead hepatocytes remaining in the drug-induced lesion is predicted if either the infiltrating macrophages are lacking phagocytotic activity (scenario 1 in Figure 3A), or are completely depleted from the lobule (scenario 3 in Figure 4A). This suggests that the Kupffer cells, which are the resident macrophages, are insufficient to clear the lesion. This is still the case, even if their dead body elimination time is reduced from 3 h⁷⁹ to 6 min (Figure S6H). The Kupffer cells are observed to engulf the dead hepatocytes⁴⁹ while the
infiltrating macrophages are hypothesized to remove the dead hepatocytes in vivo.⁵⁰ Our simulation (scenario 1 in Figure 3) strongly support this hypothesis as depleting the phagocytosis ability of infiltrating macrophages leads to incomplete regeneration. However, one may argue that if the Kupffer cells have a sufficiently strong ability to phagocytize dead cells, they might in principle clear the lesion alone. Our DT permits to specify the parameters under which this would be the case.

Seki et al.⁸⁰ observed *in vitro* that stimulated HSCs attract Kupffer cells to migrate toward them. Within our DT we could show the prospective effect of such a mechanism *in vivo*: by depleting Kupffer cells or inhibiting Kupffer cells to migrate toward the concentrated HSCs in the lesion, there are also uncleared dead hepatocytes remaining in the lesion (scenario 4 in Figure 3A and scenario 2 in Figure 4A), suggesting that the relation observed by Seki et al.⁸⁰ *in vitro* should indeed be present *in vivo*, not requiring an *in vivo* study in the first place, where the conditions are more complex and difficult to control. I.e., the simulation indicates that in such a case as for the *in vitro*-finding by Seki et al.⁸⁰ *an in vivo* validation experiment (as it had been done) is promising to perform hence the resources (e.g., time, money, material, personnel) are likely well invested. In that sense, simulations with our DT can guide the experimental strategy. When we depleted the HSCs from the lobule or disabled the activation of HSCs, for example through TGFβ produced by Kupffer cells, a similar unhealed lobule resulted (scenario 3 in Figure 3A and scenario 1 in Figure 4A). Together, the simulations suggest that the key to clear up the lesion is to attract a large enough number of macrophages in time to eliminate the dead hepatocytes. This requires the help of HSCs, e.g., to guide the migration of both Kupffer cells and infiltrating macrophages into the lesion.

### Reduced or no lesion (Figure 5, yellow scenario)

We suggested that the majority of cells dying after APAP-induced liver injury is due to the effect of NAPQI-induced cell death.⁷⁷ The role of neutrophils is complex and controversial, some studies have shown that a lack of neutrophils does not affect the outcome or severity of APAP-induced liver injury,⁴⁸ while other studies indicated that the neutrophils can directly mediate hepatocyte death in APAP-induced liver injury.⁷⁵ Therefore, the cell death induced by neutrophils in APAP-induced liver injury is not very obvious *in vivo*,⁴⁸ but can be observed *in vitro*.⁷⁵ To test this assumption, we abrogated NAPQI-induced cell death, in which we found the majority of the hepatocytes in the lesion are killed by neutrophils. This is in agreement with the indication that neutrophils can mediate hepatocyte death. Furthermore, we abrogated NAPQI-induced cell death and depleted neutrophils. The lobule was surprisingly protected with no dead hepatocytes present (scenario 4 in Figure S7A). This is another example of how our DT can be used to test if or/and within which parameter ranges certain hypothesized mechanism have an observed effect. The simulation on the influence of neutrophils exactly demonstrates where the models can be very valuable: they permit to verify whether a certain set of mechanistic hypotheses is able to explain a certain set of experimental findings, because they constitute a sort of mapping of the set of hypotheses to the result space. This is a useful information for experimentalists: if one is sure about a certain set of mechanisms, the failure in obtaining the observed experimental results and so forth may have been erroneous. If one is sure about the experimental result but did not obtain it in the DT simulation, then the hypothesized mechanisms may be incomplete or false. Quantitative biology, and even more quantitative medicine will require removal of inaccuracies and uncertainties as much as possible, or, if these remain, a better understanding of their origin.

### No effect on lesion generation nor on regeneration (Figure 5, black scenario):

As mentioned in the last point, experimental studies on the role of neutrophils in regeneration have reported conflicting observations, constituting a typical case where a DT can help to identify the expected result under well-defined hypotheses and within defined parameter ranges. Some studies report that a lack of neutrophils does not affect the outcome or severity of APAP-induced liver injury,⁸¹,⁸² while others indicated that the neutrophils can directly mediate hepatocyte death in APAP-induced liver injury⁷⁵ and lack of neutrophils does affect the severity of APAP-induced liver injury.⁷⁶ Simulations performed with the DT predict that the depletion of neutrophils has no significant impact on the regeneration of the dead cell lesion (scenario 4 in Figure 4A), supporting the conclusion of Cover et al.⁸¹ and Williams et al.⁸² This result indicates that the neutrophil-induced cell death is supplementary to the NAPQI-induced cell death. Even if NAPQI-induced cell death is disabled in the simulation, the neutrophils can still kill the Cyp450 hepatocytes to generate room for the dependent new proliferation generated hepatocytes (Figure S7). In addition, the repopulation of the regenerated liver with quiescent HSCs, subsequent to (1) reversion of activated HSCs back to the quiescent phenotype by macrophages, and (2) inducing cell death of activated HSCs by macrophages or neutrophils, has no significant effect on the lesion regeneration (scenario 5 in Figure 3A). Furthermore, we also tested the case if HSCs are activated in the early stage (right after the injury). As shown in Figure S6G, when HSCs are activated first and migrate afterward or in the case when KCs are depleted, HSCs can still be activated, both have no significant effect on lesion regeneration.

In conclusion, we have demonstrated the possible benefit and opportunities of a liver DT to simulate, in this case, liver regeneration after acute damage in control and perturbation cases, at the level of a lobule and its constituting cell types, in time and space. This may also include situations that can only be attained upon complex engineered manipulations, e.g., to explore potential therapy or protection effects. The DT can be further developed and specified to simulate physiological and pathophysiological scenarios of liver and liver diseases. For example, in order to capture complex hepatocyte shapes as it occurs during fatty liver disease, fibrosis, cirrhosis or hepatocellular cancer, the spheroidal hepatocyte (cell) model may be replaced by a cell model at higher spatial resolution.⁶⁹,⁸³ In summary, we show the potential of our *in-silico* liver to successfully simulate liver physiology, and therewith present a promising strategy toward a full liver DT, that permits to test perturbations from the molecular up to cell, tissue and body scales. We are not expecting it to replace experiments, but guiding toward the most informative experiments by identifying gaps in mechanistic knowledge. Such a liver DT is a key milestone on the route to guide diagnosis and therapy if fed with patient data. In so far, our work responds on the question of systems complexity. Our DTs explain data on a set of
subjects (here: mice) but have not been calibrated to explain data of an specific, individual subject (a specific mouse). This could be achieved
by a final calibration step of model, which is out of the scope of this work.

Limitations of the study

1) Our liver DT is based only punctually on proteomics/transcriptomics data (e.g., to characterize the cell type) and mainly on spatial-temporal staining data, as its objective was to set up a model focusing on the spatially temporally resolved interplay of cell types during regeneration after drug-induced pericentral liver damage that so far has not been addressed. In a further step, the intracellular processes shall be inte- grated. For example, if the human situation is to be considered, more multidimensional human data should be integrated to adapt the animal liver DT to the human case. We refrained from integrating molecular pathways here as this would have significantly increased the number of parameters with difficulties in parametrisation and enormous prolongation of simulation times as we experienced in another recent work. ¹¹ To arrive at a quantitative liver DT, a stepwise build-up module by module is more promising than an all-in-once DT(C) as the former promotes understanding of the different components of such a complex system. 2) Many cell types and signals have not been integrated into the DT yet. Therefore, their influence on liver regeneration has not been studied and discussed in this work. However, these groups of cells and their behaviors such as portal fibroblast activation, endothelial-mesenchymal transition, or hepatocyte myofibroblast transdifferentiation can be integrated one by one in the further process. 3) A precise determination of the threshold concentrations in the simulation was not possible as the data were not present, and as the precise orchestration of molecular factors was and could not be mimicked in great detail. The un- derlying parametrisation concept was to choose plausible values out of accessible parameter ranges and infer parameters by comparison of their effect on the regeneration process. The robustness of the so determined parameters was studied in simulated sensitivity analyses. Nevertheless, it cannot be fully excluded that some parameters could in reality be outside of the ranges. However, this is not critical: the pa- rameters/mechanisms that have been found to critically modify the regeneration process compared to the experimental observation should be challenged by additional experiments, whose outcome would then serve to re-calibrate the model parameter if necessary. Such an iter- ative procedure – identification of critical parameters and/or mechanisms by the model simulation – then experimental testing – then re-cali- bration of the liver DT, will ultimately lead to a full quantitative DT of liver generation.

STAR★METHODS

Detailed methods are provided in the online version of this paper and include the following:

● KEY RESOURCES TABLE

● RESOURCE AVAILABILITY

○ Lead contact

○ Materials availability

○ Data and code availability

● EXPERIMENTAL MODEL AND SUBJECT DETAILS

○ Mice

● **METHOD DETAILS**

○ Mouse experiments

○ Digital twin notions

○ Agent-based modeling of cells and elements

○ Gradient of signals to regulate cell behaviors

○ Scenario of lobule regeneration

● QUANTIFICATION AND STATISTICAL ANALYSIS

SUPPLEMENTAL INFORMATION

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2023.108077.

ACKNOWLEDGMENTS

The authors acknowledge financial support from BMBF LiSyM FKZ: 031L0045, BMBF LiSyM-Krebs: 031L0257D, ANR iLite: ANR-16-RHUS-0005, ANR STEDI-NASH: ANR-20-CE19-0005. A.G. was funded by the German Research Foundation (DGF; Project Ids 517010379 & 457840828). S.D. was funded by LiSyM Grant PTJ-FKZ: 031L0043.

AUTHOR CONTRIBUTIONS

The modeling part was performed by D.D. and J.Z., the experimental part by A.G., R.H., J.G.H.

**Conceptualization:** D.D.; **Data curation:** J.Z., A.G., R.H.; **Formal analysis:** J.Z., A.G., D.D., **Funding acquisition:** A.G., S.D., J.G.H., D.D.; **Investigation:** D.D., J.Z.; **Methodology:** D.D., J.Z. (modeling); A.G., J.G.H. (experimental); **Project administration:** D.D.; **Resources:** D.D., J.G.H.; **Software implementation:** J.Z.; **Supervision:** D.D., A.G.; **Validation:** D.D., J.Z., A.G., S.D., J.G.H.; **Visualization:** J.Z.; **Writing, Original Draft:** J.Z. and D.D.; **Writing, Review and Editing:** J.Z., A.G., S.D., J.G.H, D.D.
## DECLARATION OF INTERESTS

The authors declare no competing interests.

## INCLUSION AND DIVERSITY

We support inclusive, diverse, and equitable conduct of research.

Received: May 9, 2022

Revised: February 22, 2023

Accepted: September 25, 2023

Published: September 28, 2023

## REFERENCES

1. Kang, L.I., Mars, W.M., and Michalopoulos, G.K. (2012). Signals and cells involved in regulating liver regeneration. Cells 1, 1261-1292.

2. Kisseleva, T., and Brenner, D. (2021). Molecular and cellular mechanisms of liver fibrosis and its regression. Nat. Rev. Gastroenterol. Hepatol. 18, 151-166.

3. Michalopoulos, G.K., and Bhushan, B. (2021). Liver regeneration: biological and pathological mechanisms and implications. Nat. Rev. Gastroenterol. Hepatol. 18, 40-55.

4. Calderwood, S.K., Gong, J., and Murshid, A. (2016). Extracellular HSPs: the complicated roles of extracellular HSPs in immunity. Front. Immunol. 7, 159.

5. Li, H., Zhou, Y., Wang, H., Zhang, M., Qiu, P., Zhang, M., Zhang, R., Zhao, Q., and Liu, J. (2020). Crosstalk between liver macrophages and surrounding cells in nonalcoholic steatohepatitis. Front. Immunol. 11, 1169.

6. Marra, F., and Tacke, F. (2014). Roles for chemokines in liver disease. Gastroenterology 147, 577-594.e1.

7. Puche, J.E., Saiman, Y., and Friedman, S.L. (2013). Hepatic stem cells and liver fibrosis. Compr. Physiol. 3, 1473-1492.

8. Tacke, F., and Zimmermann, H.W. (2014). Macrophage heterogeneity in liver injury fibrosis. J. Hepatol. 60, 1090-1096.

9. Schliess, F., Hoehme, S., Henkel, S.G., Ghallab, A., Driesch, D., Böttger, J., Guthke, R., Pfaff, M., Hengstler, J.G., Gebhardt, R., et al. (2014). Integrated metabolic spatial-temporal model for the prediction of ammonia detoxification during liver damage and regeneration. Hepatology 60, 2040-2051.

10. Ghallab, A., Celliere, G., Henkel, S.G., Driesch, D., Hoehme, S., Hofmann, U., Zellmer, S., Godoy, P., Sachinidis, A., Blaszkewicz, M., et al. (2016). Model-guided identification of a therapeutic strategy to reduce hyperammonemia in liver diseases. J. Hepatol. 64, 860-871.

11. Dichamp, J., Cellière, G., Ghallab, A., Hassan, R., Boissier, N., Hofmann, U., Reinders, J., Sezgin, S., Zühlke, S., Hengstler, J.G., and Drasdo, D. (2023). In vitro to in vivo acetaminophen hepatotoxicity extrapolation using classical schemes, pharmaco-temporal models and a multiscale spatial-temporal liver twin. Front. Bioeng. Biotechnol. 11, 1049564.

12. Verma, B.K., Subramaniam, P., and Vadigepalli, R. (2019). Model-based virtual patient analysis of human liver regeneration predicts critical perioperative factors

controlling the dynamic mode of response to resection. BMC Syst. Biol. 13, 9.

13. Kuepfer, L., Clayton, O., Thiel, C., Cordes, H., Nudischer, R., Blank, L.M., Baier, V., Heymans, S., Caiment, F., Roth, A., et al. (2018). A model-based assay design to reproduce in vivo patterns of acute drug-induced toxicity. Arch. Toxicol. 92, 553-555.

14. Remien, C.H., Adler, F.R., Waddoups, L., Box, T.D., and Sussman, N.L. (2012). Mathematical modeling of liber injury and dysfunction after acetaminophen overdose: Early discrimination between survival and death. Hepatology 56, 727-734.

15. Naik, A., Rozman, D., and Belić, A. (2014). SteatoNet: The first integrated human metaolic model with multi-layered regulation to investigate liver-associated pathologies. PLoS Comput. Biol. 10, e1003993.

16. Hetherington, J., Sumner, T., Seymour, R.M., Li, L., Rey, M.V., Yamaji, S., Saffrey, P., Margoninski, O., Bogle, I.D.L., Finkelstein, A., and Warner, A. (2012). A composite computational model of liver glucose homeostasis. I. building the composite model. J. R. Soc. Interface 9, 689-700.

17. Friedman, A., and Hao, W. (2017). Mathematical modeling of liver fibrosis. Math. Biosci. Eng. 14, 143-164.

18. Schwen, L.O., Krauss, M., Niederalt, C., Gremse, F., Kiessling, F., Schenk, A., Preusser, T., and Kuepfer, L. (2014). Spatio-temporal simulation of first pass drug perfusion in the liver. PLoS Comput. Biol. 10, e1003499.

19. Schwen, L.O., Kuefer, L., and Preusser, T. (2016). Modeling approaches for hepatic spatial heterogeneity in pharmacokinetic simulations. Drug Discov. Today Dis. Model. 22, 35-43.

20. Lambers, L., Suditsch, M., Wagner, A., and Ricken, T. (2021). A multiscale and multiphase model of function-perfusion growth processes in the human liver. Proc. Appl. Math. and Mech. 20, e20200290.

21. Pellicer-Valero, O.J., Rupérez, M.J., Martínez-Sanchis, S., and Martín-Guerrero, J.D. (2020). Real-time biomechanical modeling of the liver using machine learning models trained on finite element method simulations. Expert Syst. Appl. 143, 113083.

22. Dutta-Moscato, J., Solovyev, A., Mi, Q., Nishikawa, T., Soto-Gutierrez, A., Fox, I.J., and Vodovotz, Y. (2014). A multiscale agent-based in silico model of liver fibrosis progression. Front. Bioeng. Biotechnol. 2, 18.

23. Means, S.A., and Ho, H. (2019). A spatial-temporal model for zonal hepatotoxicity of

acetaminophen. Drug Metabol. Pharmacokinet. 34, 71-77.

24. Adhyapok, P., Fu, X., Sluka, J.P., Clendenon, S.G., Sluka, V.D., Wang, Z., Dunn, K., Klaunig, J.E., and Glazier, J.A. (2020). A computational model of liver tissue damage and repair. PLoS One 15, e0243451.

25. Heldring, M.M., Shaw, A.H., and Beltman, J.B. (2022). Unraveling the effect of intra- and intercellular processes on acetaminophen-induced liver injury. NPJ Syst. Biol. Appl. 8, 27.

26. Sluka, J.P., Fu, X., Swat, M., Belmonte, J.M., Cosmanescu, A., Clendenon, S.G., Wambaugh, J.F., and Glazier, J.A. (2016). A liver-centric multiscale modeling framework for xenobiotics. PLoS One 11, e0162428.

27. Wambaugh, J., and Shah, I. (2010). Simulating microdosimetry in a virtual hepatic lobule. PLoS Comput. Biol. 6, e1000756.

28. Hoehme, S., and Drasdo, D. (2010). A cell-based simulation software for multi-cellular systems. Bioinformatics 26, 2641-2642.

29. Hoehme, S., Bertaux, F., Weens, W., Grasl-Kraupp, B., Hengstler, J.G., and Drasdo, D. (2018). Model prediction and validation of an order mechanism controlling the spatiotemporal phenotype of early hepatocellular carcinoma. Bull. Math. Biol. 80, 1134-1171.

30. Hoehme, S., Hammad, S., Boettger, J., Begher-Tibbe, B., Bucur, P., Vibert, E., Gebhardt, R., Hengstler, J.G., and Drasdo, D. (2023). Digital twin demonstrates significance of biomechanical growth control in liver regeneration after partial hepatectomy. iScience 26, 105714.

31. Sezgin, S., Hassan, R., Zuehlke, S., Kuepfer, L., Hengstler, J.G., Spiteller, M., and Ghallab, A. (2018). Spatio-temporal visualization of the distribution of acetaminophen as well as its metabolites and adducts in mouse livers by MALDI MSI. Arch. Toxicol. 92, 2963-2977.

32. Ghallab, A., Myllys, M., Friebel, A., Duda, J., Edlund, K., Halilbasic, E., Vucur, M., Hobloss, Z., Brackhagen, L., Begher-Tibbe, B., et al. (2021). Spatio-Temporal Multiscale Analysis of Western Diet-Fed Mice Reveals a Translationally Relevant Sequence of Events during NAFLD Progression. Cells 10, 2516.

33. Brenner, C., Galluzzi, L., Kepp, O., and Kroemer, G. (2013). Decoding cell death signals in liver inflammation. J. Hepatol. 59, 583-594.

34. Huebener, P., Hernandez, C., and Schwabe, R.F. (2015). HMGB1 and injury amplification. Oncotarget 6, 23048-23049.
<figure><img src="image_14.png" /></figure>

35. Martin-Murphy, B.V., Holt, M.P., and Ju, C. (2010). The role of damage associated molecular pattern molecules in acetaminophen-induced liver injury in mice. *Toxicol. Lett.* 192, 387–394.

36. Nowatari, T., Murata, S., Fukunaga, K., and Ohkohchi, N. (2014). Role of platelets in chronic liver disease and acute liver injury. *Hepatol. Res.* 44, 165–172.

37. Ramadori, P., Klag, T., Malek, N.P., and Heikenwalder, M. (2019). Platelets in chronic liver disease, from bench to bedside. *JHEP Rep.* 1, 448–459.

38. Meyer, J., Lejmi, E., Fontana, P., Morel, P., Gonelle-Gispert, C., and Bühler, L. (2015). A focus on the role of platelets in liver regeneration: Do platelet-endothelial cell interactions initiate the regenerative process? *J. Hepatol.* 63, 1263–1271.

39. Pinzani, M., Milani, S., Grappone, C., Weber, F.L., Jr., Gentilini, P., and Abboud, H.E. (1994). Expression of platelet-derived growth factor in a model of acute liver injury. *Hepatology* 19, 701–707.

40. Yang, C., Zeisberg, M., Mosterman, B., Sudhakar, A., Yerramalla, U., Holthaus, K., Xu, L., Eng, F., Afdhal, N., and Kalluri, R. (2003). Liver fibrosis: insights into migration of hepatic stellate cells in response to extracellular matrix and growth factors. *Gastroenterology* 124, 147–159.

41. Melton, A.C., and Yee, H.F. (2007). Hepatic stellate cell protrusions couple platelet-derived growth factor-BB to chemotaxis. *Hepatology* 45, 1446–1453.

42. De Bleser, P.J., Niki, T., Rogiers, V., and Geerts, A. (1997). Transforming growth factor-β gene expression in normal and fibrotic rat liver. *J. Hepatol.* 26, 886–893.

43. Cai, X., Li, Z., Zhang, Q., Qu, Y., Xu, M., Wan, X., and Lu, L. (2018). CXCL-EGFR-induced kupffer cells secrete TGF-β1 promoting hepatic stellate cell activation via the SMAD2/BRD4/C-MYC/EZH2 pathway in liver fibrosis. *J. Cell Mol. Med.* 22, 5050–5061.

44. Fan, W., Liu, T., Chen, W., Hammad, S., Longerich, T., Hauser, I., Fu, Y., Li, N., He, Y., Liu, C., et al. (2019). ECM1 prevents activation of transforming growth factor β, hepatic stellate cells, and fibrogenesis in mice. *Gastroenterology* 157, 1352–1367.e13.

45. Imamura, M., Ogawa, T., Sasaguri, Y., Chayama, K., and Ueno, H. (2005). Suppression of macrophage infiltration inhibits activation of hepatic stellate cells and liver fibrogenesis in rats. *Gastroenterology* 128, 138–146.

46. Kisseleva, T., and Brenner, D.A. (2007). Role of hepatic stellate cells in fibrogenesis and the reversal of fibrosis. *J. Gastroenterol. Hepatol.* 22, S73–S78.

47. Baeck, C., Wehr, A., Karlmark, K.R., Heymann, F., Vucur, M., Gassler, N., Huss, S., Klussmann, S., Eulberg, D., Luedde, T., et al. (2012). Pharmacological inhibition of the chemokine CCL2 (MCP-1) diminishes liver macrophage infiltration and steatohepatitis in chronic hepatic injury. *Gut* 61, 416–426.

48. Krenkel, O., Mossanen, J.C., and Tacke, F. (2014). Immune mechanisms in acetaminophen-induced acute liver injury. *Hepatobiliary Surg. Nutr.* 3, 331–343.

49. Canbay, A., Feldstein, A.E., Higuchi, H., Werneburg, N., Grambihler, A., Bronk, S.F., and Gores, G.J. (2003). Kupffer cell engulfment of apoptotic bodies stimulates

death ligand and cytokine expression. *Hepatology* 38, 1188–1198.

50. Boulter, L., Govaere, O., Bird, T.G., Radulescu, S., Ramachandran, P., Peliccoro, A., Ridgway, R.A., Seo, S.S., Spee, B., Van Rooijen, N., et al. (2012). Macrophage derived Wnt signaling opposes Notch signalling in a Numb mediated manner to specify HPC fate in chronic liver disease in human and mouse. *Nat. Med.* 18, 572–579.

51. Kisseleva, T., Cong, M., Paik, Y., Scholten, D., Jiang, C., Benner, C., Iwaisako, K., Moore-Morris, T., Scott, B., Tsukamoto, H., et al. (2012). Myofibroblasts revert to an inactive phenotype during regression of liver fibrosis. *Proc. Natl. Acad. Sci. USA* 109, 9448–9453.

52. Troeger, J.S., Mederacke, I., Gwak, G.Y., Dapito, D.H., Mu, X., Hsu, C.C., Pradere, J.P., Friedman, R.A., and Schwabe, R.F. (2012). Deactivation of hepatic stellate cells during liver fibrosis resolution in mice. *Gastroenterology* 143, 1073–1083.e22.

53. Hassan, R. (2017). Mechanisms of activated hepatic stellate cell removal in acute and chronic liver injury. *Thesis, Justus-Liebig-Universität Giessen*.

54. Willemsen, L., and de Winther, M.P. (2020). Macrophage subsets in atherosclerosis as defined by single-cell technologies. *J. Pathol.* 250, 705–714.

55. Schneider, K.M., Candels, L.S., Hov, J.R., Myllys, M., Hassan, R., Schneider, C.V., Wahlström, A., Mohs, A., Zühlke, S., Liao, L., et al. (2021). Gut microbiota depletion exacerbates cholestasis liver injury via loss of FXR signaling. *Nat. Metab.* 3, 1228–1241.

56. Schuran, F.A., Lommetz, C., Steudter, A., Ghallab, A., Wieschendorf, B., Schwinge, D., Zuehlke, S., Reinders, J., Heeren, J., Lohse, A.W., et al. (2021). Aryl hydrocarbon receptor activity in hepatocytes sensitizes to hyperacute acetaminophen-induced hepatotoxicity in mice. *Cell. Mol. Gastroenterol. Hepatol.* 11, 371–388.

57. Ramaiah, S.K., and Jaeschke, H. (2007). Role of neutrophils in the pathogenesis of acute inflammatory liver injury. *Toxicol. Pathol.* 35, 757–766.

58. McDonald, B., Pittman, K., Menezes, G.B., Hirota, S.A., Slaba, I., Waterhouse, C.C.M., Beck, P.L., Muruve, D.A., and Kubes, P. (2010). Intravascular danger signals guide neutrophils to sites of sterile inflammation. *Science* 330, 362–366.

59. Dragomir, A.C.D., Sun, R., Choi, H., Laskin, J.D., and Laskin, D.L. (2012). Role of Galectin-3 in classical and alternative macrophage activation in the liver following acetaminophen intoxication. *J. Immunol.* 189, 5934–5941.

60. Graubardt, N., Vugman, M., Mouhadeb, O., Caliari, G., Pasmanik-Chor, M., Reuveni, D., Zigmond, E., Brazowski, E., David, E., Chappell-Maor, L., et al. (2017). Ly6Chi monocytes and their macrophage descendants regulate neutrophil function and clearance in acetaminophen-induced liver injury. *Front. Immunol.* 8, 626.

61. Michalopoulos, G.K. (2010). Liver regeneration after partial hepatectomy: critical analysis of mechanistic dilemmas. *Am. J. Pathol.* 176, 2–13.

62. Michalopoulos, G.K. (2017). Hepatostat: Liver regeneration and normal liver tissue maintenance. *Hepatology* 65, 1384–1392.

63. Zigmond, E., Samia-Grinberg, S., Pasmanik-Chor, M., Brazowski, E., Shiboleit, O., Halpern, Z., and Varol, C. (2014). Infiltrating

monocyte-derived macrophages and resident Kupffer cells display different ontogeny and functions in acute liver injury. *J. Immunol.* 193, 344–353.

64. Ritz, T., Krenkel, O., and Tacke, F. (2018). Dynamic plasticity of macrophage functions in diseased liver. *Cell. Immunol.* 330, 175–182.

65. Fischer, R., Cariers, A., Reinehr, R., and Häussinger, D. (2002). Caspase 9-dependent killing of Hepatic Stellate Cells by activated Kupffer cells. *Gastroenterology* 123, 845–861.

66. Gao, B., and Radaeva, S. (2013). Natural killer and natural killer T cells in liver fibrosis. *Biochim. Biophys. Acta* 1832, 1061–1069.

67. Chu, Y.S., Dufour, S., Thiery, J.P., Perez, E., and Pincet, F. (2005). Johnson-Kendall-Roberts theory applied to living cells. *Phys. Rev. Lett.* 94, 028102.

68. Van Liedekerke, P., Palm, M.M., Jagiella, N., and Drasdo, D. (2015). Simulating tissue mechanics with agent-based models: concepts, perspectives, and some novel results. *Comput. Part. Mech.* 2 (4), 401–444.

69. Van Liedekerke, P., Neitsch, P., Johann, T., Alessandri, K., Nassoy, P., and Drasdo, D. (2019). Quantitative agent-based modeling reveals mechanical stress response of growing tumor spheroids is predictable over various growth conditions and cell lines. *PLoS Comput. Biol.* 15, e1006273.

70. Marques, P.E., Oliveira, A.G., Pereira, R.V., David, B.A., Gomides, L.F., Saraiva, A.M., Pires, D.A., Novaes, J.T., Patricio, D.O., Cisalpino, D., et al. (2015). Hepatic DNA deposition drives drug-induced liver injury and inflammation in mice. *Hepatology* 61, 348–360.

71. Bouwens, L., Knook, D.L., and Wisse, E. (1986). Local proliferation and extrahepatic recruitment of liver macrophages (Kupffer cells) in partial-body irradiated rats. *J. Leukoc. Biol.* 39, 687–697.

72. Wake, K. (2006). Hepatic stellate cells: Three-dimensional structure, localization, heterogeneity and development. *Proc. Jpn. Acad. Ser. B Phys. Biol. Sci.* 82, 155–164.

73. Stewart, R.K., Dangi, A., Huang, C., Murase, N., Kimura, S., Stolz, D.B., Wilson, G.C., Lentsch, A.B., and Gandhi, C.R. (2014). A novel mouse model of depletion of stellate cells clarifies their role in ischemia/reperfusion- and endotoxin-induced acute liver injury. *J. Hepatol.* 60, 298–305.

74. You, Q., Holt, M., Yin, H., Li, G., Hu, C.J., and Ju, C. (2013). Role of hepatic resident and infiltrating macrophages in liver repair after acute injury. *Biochem. Pharmacol.* 86, 836–843.

75. Marques, P.E., Amaral, S.S., Pires, D.A., Nogueira, L.L., Soriani, F.M., Lima, B.H.F., Lopes, G.A.O., Russo, R.C., Avila, T.V., Melgaço, J.G., et al. (2012). Chemokines and mitochondrial products activate neutrophils to amplify organ injury during mouse acute liver failure. *Hepatology* 56, 1971–1982.

76. Liu, Z.X., Han, D., Gunawan, B., and Kaplowitz, N. (2006). Neutrophil depletion protects against murine acetaminophen hepatotoxicity. *Hepatology* 43, 1220–1230.

77. Ghallab, A., Myllys, M., Holland, C.H., Zaza, A., Murad, W., Hassan, R., Ahmed, Y.A., Abbas, T., Abdelrahim, E.A., Schneider, K.M., et al. (2019). Influence of liver fibrosis on lobular zonation. *Cells* 8, 1556.

78. Trepat, X., Wasserman, M.R., Angelini, T.E., Millet, E., Weitz, D.A., Butler, J.P., and Fredberg, J.J. (2009). Physical forces during
collective cell migration. Nat. Phys. 5, 426-430.

79. Haecker, H., Fuermann, C., Wagner, H., and Haeker, G. (2002). Caspase-9/-3 activation and apoptosis are induced in mouse macrophages upon ingestion and digestion of Escherichia coli bacteria. J. Immunol. 169, 3172-3179.

80. Seki, E., De Minicis, S., Osterreicher, C.H., Kluwe, J., Osawa, Y., Brenner, D.A., and Schwabe, R.F. (2007). TLR4 enhances TGF-β signaling and hepatic fibrosis. Nat. Med. 13, 1324-1332.

81. Cover, C., Liu, J., Farhood, A., Malle, E., Waalkes, M.P., Bajt, M.L., and Jaeschke, H. (2006). Pathophysiological role of the acute inflammatory response during acetaminophen hepatotoxicity. Toxicol. Appl. Pharmacol. 216, 98-107.

82. Williams, C.D., Bajt, M.L., Farhood, A., and Jaeschke, H. (2010). Acetaminophen-induced hepatic neutrophil accumulation and inflammatory liver injury in CD18-deficient mice. Liver Int. 30, 1280-1292.

83. Van Liedekerke, P., Neitsch, J., Johann, T., Warmt, E., González-Valverde, I., Hoehme, S., Grosser, S., Kaes, J., and Drasdo, D. (2020). A quantitative high-resolution computational mechanics cell model for growing and regenerating tissues. Biomech. Model. Mechanobiol. 19, 189-220.

84. Holland, C.H., Ramirez, F., Myllys, M., Hassan, R., Edlund, K., Hofman, U., Marchan, R., Cadenas, C., Reinders, J., Hoehme, S., et al. (2022). Transcriptomic cross-species analysis of chronic liver disease reveals consistent regulation between humans and mice. Hepatol. Commun. 6, 161-177.

85. Campos, G., Schmidt-Heck, W., De Smedt, J., Widera, A., Ghallab, A., Pütter, L., González, D., Edlund, K., Cadenas, C., Marchan, R., et al. (2020). Inflammation-associated suppression of metabolic gene networks in acute and chronic liver disease. Arch. Toxicol. 94, 205-217.

86. Gianmoena, K., Gasparoni, N., Jashari, A., Gabrys, P., Grgas, K., Ghallab, A., Nordström, K., Gasparoni, G., Reinders, J., Edlund, K., et al. (2021). Epigenomic and transcriptional profiling identifies impaired glyoxylate detoxification in NAFLD as a risk factor for hyperoxaluria. Cell Rep. 36, 109526.

87. Schneider, K.M., Elfers, C., Ghallab, A., Schneider, C.V., Galvez, E.J.C., Mohs, A., Gui, W., Candels, L.S., Wirtz, T.H., Zuehlke, S., et al. (2021). Intestinal dysbiosis amplifies acetaminophen-induced acute liver injury.

Cell. Mol. Gastroenterol. Hepatol. 11, 909-933.

88. Wright, L., and Davidson, S. (2020). How to tell the difference between a model and a digital twin. Adv. Model. Simul. Eng. Sci. 7, 1-13.

89. Ju, C., and Tacke, F. (2016). Hepatic macrophages in homeostasis and liver diseases: from pathogenesis to novel therapeutic strategies. Cell. Mol. Immunol. 13, 316-327.

90. Rohrschneider, M., Scheuermann, G., Hoehme, S., and Drasdo, D. (2007). Shape characterization of extracted and simulated tumor samples using topological and geometric measures. 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 6271-6277.

91. Popov, V.L. (2010). Contact Mechanics and Friction (Springer).

92. Arndt, D., Bangerth, W., Blais, B., Fehling, M., Gassmöller, R., heister, T., Heltai, L., Köcher, U., Kronbichler, M., Maier, M., et al. (2021). The deal. II library, version 9.3. J. Numer. Math. 29, 171-186.

93. Goodhill, G.J. (1997). Diffusion in axon guidance. Eur. J. Neurosci. 9, 1414-1421.

94. Al-Alwan, L.A., Chang, Y., Rousseau, S., Martin, J.G., Eidelman, D.H., and Hamid, Q. (2014). CXCL1 inhibits airway smooth muscle cell migration through the decoy receptor duffy antigen receptor for chemokines. J. Immunol. 193, 1416-1426.

95. De Donatis, A., Comito, G., Buricchi, F., Vinci, M.C., Parenti, A., Caselli, A., Camici, G., Manao, G., Ramponi, G., and Cirri, P. (2008). Proliferation versus migration in platelet-derived growth factor signaling. J. Biol. Chem. 283, 19948-19956.

96. Gouwy, M., Struyf, S., Verbeke, H., Put, W., Proost, P., Opdenakker, G., and Van Damme, J. (2009). CC chemokine ligand-2 synergizes with the nonchemokine G protein-coupled receptor ligand fMLP in monocyte chemotaxis, and it cooperates with the TLR ligand LPS via induction of CXCL8. J. Leukoc. Biol. 86, 671-680.

97. Jube, S., Rivera, Z.S., Bianchi, M.E., Powers, A., Wang, E., Pagano, I., Pass, H.I., Gaudino, G., Carbone, M., and Yang, H. (2012). Cancer cell secretion of the DAMP protein HMGB1 supports progression in malignant mesothelioma. Cancer Res. 72, 3290-3301.

98. Sieber, P., Schäfer, A., Lieberherr, R., Le Goff, F., Stritt, M., Welford, R.W.D., Gatfield, J., Peter, O., Nayler, O., and Lüthi, U. (2018). Novel high-throughput myofibroblast assays identify agonists with therapeutic

potential in pulmonary fibrosis that act via EP₂ and EP₄ receptors. PLoS One 13, e0207872.

99. Mihm, S. (2018). Danger-associated molecular patterns (DAMPs): Molecular triggers for sterile inflammation in the liver. Int. J. Mol. Sci. 19, 3104.

100. Maruyama, R., Takemura, G., Aoyama, T., Hayakawa, K., Koda, M., Kawase, Y., Qiu, X., Ohno, Y., Minatoguchi, S., Miyata, K., et al. (2001). Dynamic process of apoptosis in adult rat cardiomyocytes analyzed using 48-hour videomicroscopy and electron microscopy. Am. J. Pathol. 159, 683-691.

101. Erwig, L.P., Gordon, S., Walsh, G.M., and Rees, A.J. (1999). Previous uptake of apoptotic neutrophils or ligation of integrin receptors downmodulates the ability of macrophages to ingest apoptotic neutrophils. Blood 93, 1406-1412.

102. Naito, M., Hasegawa, G., Ebe, Y., and Yamamoto, T. (2004). Differentiation and function of Kupffer cells. Med. Electron. Microsc. 37, 16-28.

103. Tangkijvanich, P., Tam, S.P., and Yee, H.F., Jr. (2001). Wound-induced migration of rat hepatic stellate cells is modulated by endothelin-1 through Rho-kinase-mediated alterations in the acto-myosin cytoskeleton. Hepatology 33, 74-80.

104. Grabher, C., Cliffe, A., Miura, K., Hayflick, J., Pepperkok, R., Rørth, P., and Wittbrodt, J. (2007). Birth and life of tissue macrophages and their migration in embryogenesis and inflammation in medaka. J. Leukoc. Biol. 81, 263-271.

105. Ramachandran, P., Pellicoro, A., Vernon, M.A., Boulter, L., Aucott, R.L., Ali, A., Hartland, S.N., Snowdon, V.K., Cappon, A., Gordon-Walker, T.T., et al. (2012). Differential Ly-6C expression identifies the recruited macrophage phenotype, which orchestrates the regression of murine liver fibrosis. Proc. Natl. Acad. Sci. USA 109, E3186-E3195.

106. De Filippo, K., Dudeck, A., Hasenberg, M., Nye, E., Van Rooijen, N., Hartmann, K., Gunzer, M., Roers, A., and Hogg, N. (2013). Mast cell and macrophage chemokines CXCL1/CXCL2 control the early stage of neutrophil recruitment during tissue inflammation. Blood 121, 4930-4937.

107. Stravitz, R.T., Ellerbe, C., Durkalski, V., Reuben, A., Lisman, T., and Lee, M.W.; Acute Liver Failure Study Group (2017). Thromobyclopenia is associated with multi-organ system failure in patients with acute liver failure. Clin. Gastroenterol. Hepatol. 14, 613-620.
## STAR★METHODS

### KEY RESOURCES TABLE

<table><thead><tr><th>REAGENT or RESOURCE</th><th>SOURCE</th><th>IDENTIFIER</th></tr></thead><tbody><tr><td><strong>Antibodies</strong></td><td></td><td></td></tr><tr><td>Anti-human smooth muscle actin monoclonal antibody, mouse</td><td>Dako</td><td>Cat#M0851; RRID: AB_10920231</td></tr><tr><td>Anti-mouse F4/80 monoclonal antibody, rat</td><td>BioRad</td><td>Cat#MCA497; RRID: AB_2335599</td></tr><tr><td>Ultra-Map anti-mouse HRP</td><td>Roche</td><td>NA</td></tr><tr><td>Ultra-Map anti-rat</td><td>Roche</td><td>NA</td></tr><tr><td><strong>Chemicals, peptides, and recombinant proteins</strong></td><td></td><td></td></tr><tr><td>Acetaminophen</td><td>Sigma-Aldrich</td><td>Cat#A7085-500G; RRID: SCR_003986</td></tr><tr><td>ROTI Histofix 4%</td><td>ROTH</td><td>Cat#P087.3</td></tr><tr><td><strong>Experimental models: Organisms/strains</strong></td><td></td><td></td></tr><tr><td>C57BL/6NRj Mouse</td><td>JANVIER LABS</td><td>NA</td></tr><tr><td><strong>Software and algorithms</strong></td><td></td><td></td></tr><tr><td>R</td><td>R core</td><td>https://www.r-project.org/</td></tr><tr><td>ParaView</td><td>Kitware</td><td>https://www.paraview.org/</td></tr><tr><td>TiSim</td><td>This paper</td><td>https://zenodo.org/record/8252951</td></tr></tbody></table>

### RESOURCE AVAILABILITY

#### Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Dirk Drasdo (dirk.drasdo@inria.fr).

#### Materials availability

This study did not generate new unique reagents.

#### Data and code availability

* All data reported in this paper will be shared by the lead contact upon request.

* The simulation source code used in this study and the R code to analyze the data are available online via zenodo: https://doi.org/10.5281/zenodo.8252951.

* Any additional information required to reanalyse the data reported in this paper is available from the lead contact upon request.

### EXPERIMENTAL MODEL AND SUBJECT DETAILS

#### Mice

Male, 10-week-old, C57BL6/n mice (Janvier Labs, France). The animals were housed in standard environmental conditions with 12-h dark/light cycles and were fed ad libitum on a normal rodent diet (Ssniff, Soest, Germany) with free access to water.⁸⁴ All experiments were approved by the local animal welfare committee (LANUV, North Rhine-Westphalia, Germany, application number: 84–02.04.2016.A279).

### METHOD DETAILS

#### Mouse experiments

##### Induction of acute liver injury by acetaminophen

Acetaminophen (APAP)-induced acute liver injury was done in mice. Briefly, a dose of 300 mg/kg APAP was administered intraperitoneally into overnight fasted mice. APAP was dissolved in phosphate-buffered saline (PBS) with an application volume of 30 mL/kg. The mice were fed ad libitum after APAP administration. Liver tissue samples were collected time-dependently, as indicated in the results section, and were processed for immunohistochemistry.⁸⁵,⁸⁶ Briefly, liver tissue samples of ~7 × 10 mm area were collected from the left liver lobe and immediately fixed in paraformaldehyde (PFA; ROTI Histofix 4%) for two days at 4°C. Subsequently, the tissue samples were washed in PBS and were prepared for paraffin embedding using a spin tissue processor STP 120 embedding automate (Microm, walldorf, Germany). Finally, solid paraffin blocks were prepared using EC 350-modular tissue embedding center (Microm, walldorf, Germany).
## Immunohistochemistry

Immunostaining of macrophages and activated hepatic stellate cells was performed in 4 µm-thick PFA-fixed paraffine-embedded liver tissue sections using an autostainer (Discovery Ultra Automated Slide Preparation System, Roche, Germany), using the technique described in Schneider et al.⁸⁷ Briefly, rat anti-mouse F4/80 monoclonal antibody (dilution 1: 50) and mouse anti-human smooth muscle actin monoclonal antibody (dilution 1: 100) were used for binding to macrophages and activated stellate cells, respectively (key resources table). Appropriate Ultra-Map anti-rat or Ultra-Map anti-mouse HRP were used as secondary antibodies (key resources table). Representative images were acquired from the stained tissues using a brightfield microscope (Olympus, Hamburg, Germany).

## Digital twin notions

### Virtual twin, digital twin, digital twin candidate

We here consider a "virtual twin" of a liver a(n) (abstracted) as a "copy" of the real liver on the computer, which prospectively may converge to an object on the computer representing all liver functions. As the computer internally represents the virtual twin in terms of digits (i.e., the digital twin is an executable virtual twin of the physical "liver" or part of it⁸⁸), we do here not distinguish between a virtual and a digital twin, and use within this work the term "digital twin". A (computational or, simulation) "model" may just be a mapping from a certain input to an output without necessarily reflecting the alphabet of components of the liver (so may not be a "twin" in our nomenclature), while in our case we attempt to build upon the same components and mechanistic relationships between them. Candidates for a digital twin (DT) (here called "digital twin candidates" (DTCs) to stress the objective of using the same alphabet as the reality (while it is still a model, not an animation) and it is (finally) processed on the computer), that are disproved by data, are still models but disqualify for a "twin". The selection criterion is hence the comparison to data: those digital twin candidates (DTCs) that agree with data may all be considered as digital twins (DTs). Finally, the DT emerges from comparison with data from many mice in a statistical sense i.e., average values associated with an error bar, hence represents a population of mice (as opposed to, for example, a specific mouse).

### Digital twin candidates and digital twins: Modeling strategy

Our liver DTCs implement alternative hypothesized mechanisms of cellular interplay during liver regeneration and studies the consequences on readout parameters that partially are already experimentally tracked, either quantitatively or qualitatively (as e.g., the occurrence and size of the dead cell lesion, or spatial distributions of certain cell types in time), or may be informative to track, e.g., the concentration of signaling molecules. Moreover, they permit studies of the sensitivity of readout parameters of interest on experimental manipulations, for example the suppression of neutrophils. Theoretical variation of parameters and variables in the DTCs can help to identify those model parameters (or mechanisms), for which a modification would have the largest influence on the observables (e.g., hepatic stellate cells). We first identify the DTs as the subset of DTCs that agree to a certain set of data describing a normal regeneration process. One of them is considered as reference. We then apply different types of perturbations to this specific (reference) DT. Furthermore, we run perturbation simulations to predict the influence of cell type depletions on liver regeneration. Our DTCs resolve liver microarchitecture, representing each cell individually as a basic modeling element with realistic cell-biomechanics. By comparison with experimental readout parameters, we were able to identify hypotheses that were incompatible with data, which eliminates the respective DTC as candidate for a DT, and to identify perturbation experiments that are expected to critically modify the regeneration process (schematized in Figure 5). For example, in case of neutrophil suppression, Cover et al.,⁸¹ Williams et al.,⁸² reported no significant effect on severity of APAP-induced liver injury, while Liu et al.⁷⁶ claimed a significant influence. Marques et al.⁷⁵ observed hepatocyte death triggered by neutrophils in vitro. Simulations with our liver DT suggest that the impact of neutrophils in a normal liver should be negligible, unless in the purely hypothetical cases where APAP would not be metabolized by NAPQI, or where hepatocytes would not express Cytochrome P450 enzymes. Another example is the migration of Kupffer cells. A previous study has reported that Kupffer cells "are not suited" to migrate to the injury site.⁸⁹ However, Seki et al.⁸⁰ reported from in vitro experiments that the presence of HSCs does promote the migration of Kupffer cells toward HSCs. We tested both alternative mechanisms, whether the presence of HSCs promotes or does not promote the migration of Kupffer cells in vivo, in DTC simulations. Our results show for the implemented cell-cell-interaction network that if the presence of HSCs does not promote the migration of Kupffer cells, then regeneration is not complete. I.e. only the mechanism in which HSCs promote migration of Kupffer cells qualifies for a realistic DT. Hence a follow-up series of validation experiments could focus on exactly those two items.

## Agent-based modeling of cells and elements

To capture their approximate shape, hepatocytes, macrophages, and neutrophils are approximated as spheres within a "center-based model" (CBM)-approach, which mimics the forces between cells as forces between their centers.⁶⁸ Within the CBM-concept, the spherical shape does not represent the precise cell shape but region in space where the cell is located at with overwhelming probability. An HSC is modeled as a sphere forming the cell's core body, connected to four semi-flexible chains of spheres; sinusoids are modeled as semi-flexible chains of spheres. Such spheres can for example readily be constructed from inscribing spheres of maximal radius into image volume datasets stained for endothelial cells (e.g., by CD31).⁹⁰ During each time step $\Delta t$, the velocity $\vec{v}_i$ of each sphere $i$ is calculated according to a Langevin equation (an overdamped stochastic equation) of motion. The new position of $i$ is updated to $\vec{p}_i(t) + \vec{v}_i(t)\Delta t$, where $\vec{p}_i(t)$ is the position of $i$ at time $t$. Since HSCs, macrophages and neutrophils are either all located inside sinusoids, wrap around, or migrate along sinusoids, we assume in the simulations that they always remain in contact with a sinusoid edge (for HSC, the core body is sticking on the sinusoid edge). Therefore, the
vector $\vec{v}_i \Delta t$ for an HSC, macrophage, or neutrophil $i$ is projected onto the line of sinusoid edge to which $i$ is in contact with. The position of $i$ is then updated as $\vec{p}_i + \vec{e}_i(\vec{v}_i \Delta t \cdot \vec{e}_i)$, where $\vec{e}_i$ is the unit orientation vector of the local tangent to the sinusoid to which $i$ is connected with. The equations of motion for each type of cell and element are defined below.

### Equation of motion for hepatocytes

Each hepatocyte is represented as a homogeneous isotropic, elastic, adhesive sphere. It can migrate, grow, divide, and interact with other cells or sinusoids. The position of hepatocyte $i$ is updated from:

$$ \Gamma_{ECM,i} \vec{v}_i + \sum_j \Gamma_{i,j} (\vec{v}_i - \vec{v}_j) = \sum_j \vec{F}_{ij} + \vec{F}_{mig,i} \quad (\text{Equation 1}) $$

where $\Gamma_{ECM,i}$ is the friction coefficient with the extracellular matrix (ECM, which is not explicitly modeled in this study), $\Gamma_{i,j} = \gamma_{\perp} (\vec{e}_{ij} \otimes \vec{e}_{ij}) + \gamma_{\parallel} (I - \vec{e}_{ij} \otimes \vec{e}_{ij})$ is the friction tensor between cell $i$ and the sphere $j$ of other cell type (e.g., hepatocyte $j$) or sinusoid sphere (as sinusoids are made of spherical elements), and $\vec{e}_{ij}$ is the unit vector from $i$ toward $j$, $\vec{F}_{ij}$ is the corresponding central repulsion/adhesion interaction force, $\vec{F}_{mig,i}$ is an (active) migration force of cell $i$. The central force is computed by$^{91}$:

$$ \vec{F}_{ij} = \left( \frac{4\hat{E}}{3\hat{R}} [a(d_{ij})]^3 - \sqrt{8\pi\sigma\hat{E}} [a(\delta_{ij})]^3 \right) \vec{e}_{ij} \quad (\text{Equation 2}) $$

where the contact radius $a(\delta_{ij})$ allows to compute hepatocyte-hepatocyte contact area, and can be obtained by $\delta_{ij} = \frac{\hat{e}_i^2}{\hat{R}^2} - \sqrt{\frac{2\pi\sigma}{E}} \hat{E}$ and $\hat{R}$ are defined as $\hat{E}_0 = (\frac{1-r_i^2}{E_i} + \frac{1-r_j^2}{E_j})^{-1}$, $\hat{d}_{ij} = 1 - d_{ij}/(r_i + r_j)$, $\hat{E} = \begin{cases} \hat{E}_0, 0 \le \hat{d}_{ij} \le 0.08 \\ a_0 + a_1 \hat{d}_{ij} + \dots + a_6 \hat{d}_{ij}^6, 0.08 \le \hat{d}_{ij} \end{cases}$ (This is to consider the limited cell volume compressibility in a pairwise cell-cell interaction force. See more details of choosing the values of $\alpha_0, \dots, \alpha_6$ in Van Liedekerke et al.$^{69}$) and $\hat{R} = (\frac{1}{R_i} + \frac{1}{R_j})^{-1}$, with $E_i$ and $E_j$ being the Young's moduli, $v_i$ and $v_j$ the Poisson ratios, and $R_i$ and $R_j$ the radii of $i$ and $j$, $d_{ij}$ is the distance between $i$ and $j$. The migration force is computed by $\vec{F}_{mig,i} = f_{dir} \vec{e}_i + \sqrt{2\Gamma_{ECM,i}^2 D_i} \cdot \vec{\eta}_i$, where $f_{dir}$ is a constant force magnitude, $\vec{e}_i$ is the unit vector from $i$ toward the central vein, $D_i$ is the diffusion constant of $i$, $\vec{\eta}_i$ is an uncorrelated noise term with amplitude $(\eta_{in}(t)\eta_{jm}(t')) = \delta_{ij}\delta_{mn}\delta(t - t')$, $t$ and $t'$ denote times, and $m, n \in (x, y, z)$ denote the coordinates (see more details in Hoehme et al.$^{28}$).

### Equation of motion for macrophages (Kupffer cells and infiltrating macrophages) and neutrophils

Macrophages and neutrophils $i$ are represented as point objects not interacting with any other structure but capable of migrating along the sinusoids (if not otherwise stated). The position of $i$ is updated by solving the following equation:

$$ \Gamma_{ECM,i} \vec{v}_i + \sum_j \Gamma_{i,j} (\vec{v}_i - \vec{v}_j) = \sum_j \vec{F}_{ij} + \vec{F}_{mig,i} \quad (\text{Equation 3}) $$

where $\Gamma_{ECM,i}$ is the friction coefficient with the extracellular matrix, $\Gamma_{i,j}$ is the friction tensor between $i$ and $j$ (same as Equation 1), $\vec{F}_{ij}$ is the corresponding central repulsion/adhesion interaction force, $\vec{F}_{mig,i}$ is the migration force to drive $i$ to migrate. The interaction force in Equation 3 does not play an important role as compared to the densely distributed hepatocytes, macrophages and neutrophils are distributed much more sparsely in the lobule, and they are much smaller than hepatocytes.

### Equation of motion for HSCs

The core body of HSCs is modeled as a homogeneous isotropic, elastic, adhesive sphere (mainly representing the HSC's nucleus) with several semi-flexible chains of spheres as "arms" (to mimic the long HSC's protruding branches). The position of an HSC $i$ (Figure S1C) is updated by solving the following equation of motion:

$$ \Gamma_{ECM,i} \vec{v}_i + \sum_j \Gamma_{i,j} (\vec{v}_i - \vec{v}_j) = \sum_j \vec{F}_{ij} + \vec{F}_{mig,i} + \sum_k \vec{F}_{ela,ik} \quad (\text{Equation 4}) $$

where $\Gamma_{ECM,i}$ is the friction coefficient with ECM, $\Gamma_{i,j}$ is the friction tensor between $i$ and another cell or element $j$ (same as Equation 1), $\vec{F}_{ij}$ is the interaction force between $i$ and another cell or element $j$ (same as Equation 2). $\vec{F}_{ela,ik}$ represents elastic force between the head sphere of $i$ and its connected arm spheres $k$ (see definition in Figure S1C). $\vec{F}_{mig,i}$ is the migration force to drive $i$ to migrate. The equation of motion for an arm sphere $i'$ of HSC $i$ (blue sphere in Figure S1C) is approximated by:

$$ \Gamma_{ECM,i'} \vec{v}_{i'} = \sum_{k'} \vec{F}_{ela,i'k'} \quad (\text{Equation 4b}) $$

where $k'$ denote the connected spheres of arm sphere $i'$. The arms are represented to permit representation of the HSC shape as well as of potential direct cell-cell communication through the arms. However, in reality, these arms are so small that within this work their interaction with other cells is neglected.
### Equation of motion for sinusoids

Sinusoids are modeled as semi-flexible chain of spheres.²⁸ For each sinusoid sphere i, the position of i is updated by solving the following equation of motion:

$$ \Gamma_{ECM,i} \vec{v}_i = \sum_j (\Gamma_{\parallel,SE}(\vec{w}_{ij} - \vec{e}_{ij}(\vec{w}_{ij} \cdot \vec{e}_{ij})) + \vec{F}_{ij}) + \sum_k (\Gamma_{\parallel,SS}(\vec{w}_{ik} - \vec{e}_{ik}(\vec{w}_{ik} \cdot \vec{e}_{ik})) + \vec{F}_{ik}) + \vec{F}_{i,ela} \quad (\text{Equation 5}) $$

where $\Gamma_{ECM,i}$ is the friction coefficient with environment, $\Gamma_{\parallel,SE}$ denotes the longitudinal friction between sinusoid sphere i and its interacting sphere, for example hepatocyte j, $\vec{w}_{ij} = \vec{v}_j - \vec{v}_i$ is the difference of velocity between i and j, $\vec{e}_{ij}$ is the unit direction vector from i toward j, $\vec{F}_{ij}$ is the interaction force between i and j (same as Equation 2), $\Gamma_{\parallel,SS}$ denotes the longitudinal friction between two sinusoid spheres i and k, $\vec{w}_{ik} = \vec{v}_k - \vec{v}_i$ is the difference of velocity between i and k, $\vec{e}_{ik}$ is the unit direction vector from i toward k, $\vec{F}_{ik}$ is the interaction force between i and k, $\vec{F}_{i,ela}$ is the spring force that arises from the chain connections between spheres belonging to the same sinusoid.

### Gradient of signals to regulate cell behaviors

In our model, each type of molecular signal is produced by one or more certain types of cells. The dynamic of signal i released by cell j is governed by a partial differential equation (PDE):

$$ \frac{\partial \varphi_i}{\partial t} = \nabla(D_i \nabla \varphi_i) + s_i \left( \sum_j \delta(x - x_j) \right) - \gamma_i \varphi_i \quad (\text{Equation 6}) $$

where $\varphi_i$ is the density of i, $D_i$ is the diffusion coefficient of i, $s_i$ is the production rate of i, $x_j$ denotes the position of cell j, which produces i, $\gamma_i$ is the decay rate of i. The center of mass of cell j is set as the origin of the source. The simulation domain $\Omega$ is set as a box large enough to contain the entire lobule. The diffusion process of i inside $\Omega$ is assumed to be isotropic and homogeneous. We assume Dirichlet boundary condition $\varphi_i(x) = 0$ for $x \in \partial\Omega$.

### The cubic system to approximate the signal gradient

We used a simple cubic system to approximate the concentration of all signals. The entire lobule is located in a big cube which is divided into $N^3$ small element cubes (Figure S1A, N is number of element cubes on each axis). The concentration of signal i at cube i, j, k at time t is denoted as $c_{i,j,k}^t$. After the time lapse of $\Delta t$, the concentration is denoted as $c_{i,j,k}^{t+\Delta t}$. Then the solution of Equation 6 can be explicitly approximated as:

$$ \frac{c_{i,j,k}^{t+\Delta t} - c_{i,j,k}^{t}}{\Delta t} = $$

$$ D_i \left( \frac{c_{i:1-j,k}^t - 2c_{i:j,k}^t + c_{i:j+1,k}^t}{\Delta x^2} + \frac{c_{i:j-1,k}^t - 2c_{i:j,k}^t + c_{i:j+1,k}^t}{\Delta y^2} + \frac{c_{i:j,k-1}^t - 2c_{i:j,k}^t + c_{i:j,k+1}^t}{\Delta z^2} \right) + s_i - \gamma_i c_{i:j,k}^t $$

where $D_i$, $s_i$, $\gamma_i$ are the diffusion coefficient, production rate, and decay rate of i. Since we use the cube as finite element, $\Delta x$, $\Delta y$, $\Delta z$ are the same as the width of the elementary cube. To verify our cubic system, we test one simple example of placing an injured hepatocyte which produces DAMP in the center of the system. Then we solve Equation 6 by using our method and the software deal.II⁹² for the numerical solution. As shown in Figure S1B, our method generates a good approximation of the solution generated by the deal.II library. (We do not use deal.II directly in our model as it turns out to be too slow.)

### Diffusion rate

In our model, the diffusion rate of a molecule is scaled approximately as the inverse of the cubic root of the molecular weight following Good-hill.⁹³ Taking the diffusion rate of TGFβ, $D_{TGF\beta}$ as reference, the diffusion rate $D_i$ of molecule i is then approximated as $D_i = D_{TGF\beta} (\frac{w_{TGF\beta}}{w_i})^{1/3}$, where $w_{TGF\beta}$ is the molecular weight of TGFβ, $w_i$ is the molecular weight of i. The diffusion rates of all molecular signals in our digital twin candidates are listed in Table S1.

### Production rate

The concentrations of above molecules in the in vitro studies in serum or cell are around 1–20 ng/mL (Al-Alwan et al.,⁹⁴ De Donatis et al;⁹⁵ Gouwy et al;⁹⁶ Jube et al;⁹⁷ Sieber et al;⁹⁸). In our model, the concentration unit of any molecule type is assumed to be 1 which equals to 5 ng/mL, which we use as reference scale unit to non-dimensionalize our reaction-diffusion equation. The production rate $s_i$ of molecular i is arbitrarily fixed as 0.1.

### Decay rate

In our model, the decay rate $\gamma_i$ of molecule i is approximated according to its half-life time: $\gamma_i = \ln 2/T_{i,1/2}$, where $T_{i,1/2}$ is the half-life time of molecule i. The decay rate of all molecular signals in our digital twin candidates is listed in Table S1.
## Scenario of lobule regeneration

In our model, the cell behaviors are regulated by certain type of signals. Different types of cells communicate with each other by responding to specific signals and collaborate to achieve the clearance of the necrotic region and recovery from the toxin-induced injury.

### Behaviors of hepatocytes

We assume that the hepatocytes within a circle of 164 µm radius (the lobule lesion size due to CCl₄-induced injury is taken from Hoehme et al.²⁶) around the central vein of the lobule are CCl₄-induced injured cells. In our model, they are marked as *injured* hepatocytes, which can synthesize DAMPs⁴,⁵ to activate, among others, Kupffer cells and infiltrating macrophages.³⁵,⁹⁹ Activated Kupffer cells and activated HSCs can produce CXCL1 to attract neutrophils to migrate into the lesion area,⁶ and CCL2 to attract infiltrating macrophages.

We have observed that the majority of the hepatocytes are dead 2 h after the administration of APAP *in vitro* (experimental data not shown here). A previous study also showed that neutrophils can induce necrosis of hepatocytes upon administration of APAP.⁷⁵ In our model, we assign random waiting times for each *injured* hepatocyte *i* according to a Gamma distribution Γ(α = 5, β = 5). After that waiting time, *i* is marked as *dead* (necrotic) hepatocyte. The choice of α is because the necrotic pathway triggered for hepatocytes due to AIL1 (acetaminophen-induced liver injury) involves about 5 reactions¹¹ and the choice of β is to fit the data that almost all hepatocytes are dead 2 h after the administration of APAP i.e., the waiting time distribution until the death event is peaked and not monotonically decreasing as for example in case of a Poissonian distribution. In addition to APAP, neutrophils can also induce necrosis of hepatocytes. For any *injured* hepatocyte *i*, if its distance from a neutrophil *j* satisfies $d_{ij} < r_i + r_j$, it is marked as *dead* (necrotic) hepatocyte (neutrophil can evoke the necrosis of a hepatocyte⁵). We assume that any dead hepatocyte *i* is eliminated by macrophages at least 24 h after it is killed (it takes roughly 24 h for a dying cell to loose its membrane integrity and collapse into fragments¹⁰). Since it usually takes 2–4 h for a macrophage to engulf and degrade foreign objects such as dead cell bodies,⁷⁹ we simulate the degradation process of the dead body of a hepatocyte *i*, if it gets contact with a Kupffer cell or infiltrating macrophage *j* with Ly6C-high phenotype, as follows: if the distance between *i* and *j* satisfies $d_{ij} < r_i + r_j$, *i* is removed from the system 3 h after the contact.

### Behaviors of Kupffer cells

We assume that for a Kupffer cell *i*, if its local concentration of DAMP, $\varphi_{DAMP,i}$ is higher than a threshold $\varphi_{DAMP,activate}$, *i* is activated and can synthesize TGFβ, CXCL1, CCL2, and PDGF.⁶,³⁹,⁴²,⁴⁸ A previous study has reported that due to their highly stationary behavior, Kupffer cells alone "are not suited" to migrate to the injury site,⁸⁹ but the presence of HSCs does promote the migration of Kupffer cells toward HSCs.⁸⁰ We assume that an activated Kupffer cell *i* can migrate toward an activated HSC *j* if the distance between them satisfies $d_{ij} < r_i + r_j + I_{HSC,branch}$, where $I_{HSC,branch}$ is the length of the HSC branch. A migration force $\vec{F}_{mig,i} = f_{mig,i} \frac{\vec{e}_i}{||\vec{e}_i||}$, $f_{mig,i} \sim \mathcal{N}(F_{mig,KC,mean}, F_{mig,KC,sd})$ on *i* is then added to the equation of motion of cell *i*, where $\vec{e}_i$ is the orientation unit vector from *i* to *j*, $F_{mig,KC,mean}$ and $F_{mig,KC,sd}$ are the mean and standard deviation of migration force magnitude approximated from the mean and standard deviation of migration speed of Kupffer cells (due to the highly stationary behavior of Kupffer cells,⁸⁹ we arbitrarily set a low migration speed for them. We also test higher and even lower migration speed of Kupffer cells in Figure S3) as $F_{mig,KC,mean} = \Gamma_{ECM,KC}v_{KC,mean}$, $F_{mig,KC,sd} = \Gamma_{ECM,KC}v_{KC,sd}$ (see value of $v_{KC,mean}$ and $v_{KC,sd}$ in Table S1). Since dead hepatocytes can be engulfed by Kupffer cells⁴⁹ and the phagocytosis ability of macrophages is decreased by 50% after it uptakes dead cells,¹⁰¹ we assign a phagocytosis probability $p_{phag,i}$ for each activated Kupffer cell *i*, where $p_{phag,i}$ is initiated as 1. If the distance between *i* and a dead hepatocyte *j* satisfies $d_{ij} < r_i + r_j$, and a random number sampled from uniform distribution U(0, 1) is less than $p_{phag,i}$, then *i* stops moving and remains with hepatocyte *j* for 3 h to engulf and degrade the dead hepatocyte. $p_{phag,i}$ is then divided by two to mimic the decreased phagocytosis ability of *i*. If there are more than one dead hepatocyte in contact with *i*, only one of them is randomly selected as the one to be engulfed by *i*. After 3 h, *j* is removed from the system and *i* can move again. In our model, we consider two alternative time courses for Kupffer cells. Either their population size drops upon activation, recovering from day 2 on, or the Kupffer cell population remains constant. Since in the latter case the life span of a Kupffer cell is no longer than about 4 days, as reported in Naito et al.,¹⁰² we assume that each Kupffer cell would revert to a quiescent mode after it is activated.

### Behaviors of HSCs

PDGF has been reported to attract HSCs to migrate.⁴⁰,⁴¹ We assume that for any HSC *i*, if the concentration of PDGF at the location of *i*, $\varphi_{PDGF,i}$ becomes higher than a threshold $\varphi_{PDGF,migrate}$, a migration force $\vec{F}_{mig,i} = f_{mig,i} \frac{\nabla\varphi_{PDGF,i}}{||\nabla\varphi_{PDGF,i}||}$, $f_{mig,i} \sim \mathcal{N}(F_{mig,HSC,mean}, F_{mig,HSC,sd})$ is applied on *i*, where $F_{mig,HSC,mean}$ and $F_{mig,HSC,sd}$ are the mean and standard deviation of migration force magnitude approximated from the mean and standard deviation of migration speed of HSC¹⁰³ as $F_{mig,HSC,mean} = \Gamma_{ECM,HSC}v_{HSC,mean}$, $F_{mig,HSC,sd} = \Gamma_{ECM,HSC}v_{HSC,sd}$ (see value of $v_{HSC,mean}$, $v_{HSC,sd}$ in Table S1). Every HSC is initiated as in *quiescent* mode. TGFβ has been reported to activate HSCs, but it was also suggested that it is not TGFβ but a factor produced by infiltrating macrophages that activates HSCs.⁴⁵ Here we assume two alternative mechanisms for HSC activation. One is that HSCs are activated by TGFβ produced by Kupffer cells. If the local concentration of TGFβ, $\varphi_{TGF\beta,i}$ at HSC *i* is higher than a threshold $\varphi_{TGF\beta,activate}$, *i* is activated into the *activated* mode. The other mechanism is that HSCs are activated by a factor produced by infiltrating macrophages (the type of this factor is not known, so we assume it has the same diffusion constant and decay rate of TGFβ). If the local
concentration of this factor, $\phi_{factor,i}$ at HSC $i$ is higher than the threshold $\phi_{TGF\beta,activate,i}$, $i$ is activated into activated mode. Upon activation, HSCs can synthesize CXCL1 and CCL2.<sup>46,47</sup> HSCs remain activated even if TGFβ falls below the threshold again. As reported before, activated HSCs are phagocytosed by Ly6C-low phenotype infiltrating macrophages in the later stage of liver regeneration<sup>8</sup> (mechanism 1). In addition to phagocytosis, we assume another mechanism (mechanism 2) for the fate of activated HSCs: Ly6C-low phenotype infiltrating macrophages can revert activated HSCs to the quiescent mode. In our model, these two mechanisms are simulated as follows: if the distance between an aHSC $i$ and an infiltrating macrohage $j$ satisfies $d_{ij} < r_i + r_j$, $i$ is eliminated under mechanism (1) or reverted to quiescent mode again under mechanism (2). To keep the same density of HSCs after liver regeneration, we assume that under mechanism (1), once $i$ is activated, it grows and divides after 24 hours.

## Behaviors of infiltrating macrophages

DAMPs are suggested to activate macrophages.<sup>99</sup> We assume that for any infiltrating macrophage $i$, if the local concentration of DAMP at the location of $i$, $\phi_{DAMP,i}$ is higher than a threshold $\phi_{DAMP,activate,i}$, $i$ is activated and can synthesize CCL2, and PDGF.<sup>39,47</sup> During liver regeneration, macrophages infiltrate into the liver following the gradient of the chemoattractant CCL2.<sup>47</sup> We assume that for any infiltrating macrophage $i$, if its local concentration of CCL2, $\phi_{CCL2,i}$ is higher than a threshold $\phi_{CCL2,migrate}$, a migration force $\vec{F}_{mig,i} = f_{mig,i} \frac{\nabla \phi_{CCL2,i}}{|\nabla \phi_{CCL2,i}|}$, $f_{mig,i} \sim \mathcal{N}(F_{mig,IM,mean}, F_{mig,IM,sd})$ is applied on $i$, where $F_{mig,IM,mean}$ and $F_{mig,IM,sd}$ are the mean and standard deviation of the migration force magnitude, approximated from the mean and standard deviation of migration speed of infiltrating macrophages<sup>104</sup> as $F_{mig,IM,mean} = \Gamma_{ECM,IM} v_{IM,mean}$, $F_{mig,IM,sd} = \Gamma_{ECM,IM} v_{IM,sd}$ (see value of $v_{IM,mean}$, $v_{IM,sd}$ in Table S1). Each infiltrating macrophage is initiated as Ly6C-high phenotype, which is responsible of degrading necrotic cells.<sup>8</sup> After a certain duration it transforms into a Ly6C-low phenotype, which is responsible of anti-fibrosis.<sup>8,105</sup> If the distance between an infiltrating macrophage $i$ of Ly6C-high phenotype and a dead hepatocyte $j$ satisfies $d_{ij} < r_i + r_j$, $j$ is removed 3 hours later to mimic the process of engulfment and degradation of apoptotic hepatocytes.<sup>50</sup> The removal process of macrophages is as described for Kupffer cells in section "behaviors of Kupffer cells" above: A phagocytose probability $p_{phag,i}$ is assigned to an infiltrating macrophage $i$ to mimic the decreased phagocytosis ability.<sup>101</sup> We assume two alternative mechanisms of how the infiltrating macrophages could change the fate of an aHSC. Mechanism (1): infiltrating macrophage removes the aHSC. Mechanism (2): infiltrating macrophage reverts the aHSC to quiescence. In our model, if the distance between an infiltrating macrophage $i$ of Ly6C-low phenotype and an activated HSC $j$ satisfies $d_{ij} < r_i + r_j$, we generate a random number sampled from $U(0, 1)$. If this random number is less than $p_{phag,i}$, $j$ is removed under mechanisms (1) or reverted to quiescent mode under mechanism (2). If mechanism (1) is used, the value of $p_{phag,i}$ is then halved. The infiltrating macrophages are initiated as monocytes with Ly6C-high phenotype, and after a period of time between 2 and 3 days they transform into macrophages with Ly6C-low phenotype.<sup>60,63</sup> The transforming time from Ly6C-high phenotype to Ly6C-low phenotype is chosen as 72 h according to the observation from Zigmond et al.<sup>63</sup> (the mass peak of Ly6C-low phenotype macrophages is at 3 days after the injury). The Ly6C-high phenotype may not be stained by F4/80 (as shown in Figure 2 of Dragomir et al;<sup>59</sup>). Hence, at day 2 after the injury, F4/80 staining shows only the Ly6C-low phenotype, which are the KCs (Figure 1A). The lifetime of an infiltrating macrophage in the system is set to be 96 hours (estimated from the observation from the APAP-induced liver injury done by Zigmond et al.,<sup>63</sup> where the mass peak of infiltrating macrophages is at 1 day after the injury and they are hardly to be seen 5 days after the injury). Infiltrating macrophages are added into the lobule system with uniform rate during the first 24 hours. The initial position of an infiltrating macrophage is randomly sampled inside the lobule but outside of the lesion (164 µm away from the central vein) reflecting their approximate distribution at the time when they start to migrate on response to CCL2, whose peak is at about day 1.

## Behaviors of neutrophils

Neutrophils are recruited by CXCL1 at the early stage of liver injury.<sup>106</sup> We assume that for any neutrophil $i$, if the local concentration of CXCL1 at the location of $i$, $\phi_{CXCL1,i}$ is higher than a threshold $\phi_{CXCL1,migrate}$, a migration force $\vec{F}_{mig,i} = f_{mig,i} \frac{\nabla \phi_{CXCL1,i}}{|\nabla \phi_{CXCL1,i}|}$, $f_{mig,i} \sim \mathcal{N}(F_{mig,Neutrophil,mean}, F_{mig,Neutrophil,sd})$ is applied on $i$, $F_{mig,Neutrophil,mean}$ and $F_{mig,Neutrophil,sd}$ are the mean and standard deviation of the migration force magnitude, approximated from the mean and standard deviation of migration speed of neutrophils<sup>70</sup> as $F_{mig,Neutrophil,mean} = \Gamma_{ECM,Neutrophil} v_{Neutrophil,mean}$, $F_{mig,Neutrophil,sd} = \Gamma_{ECM,Neutrophil} v_{Neutrophil,sd}$ (see value of $v_{Neutrophil,mean}$, $v_{Neutrophil,sd}$ in Table S1). If its distance to an injured hepatocyte $j$ satisfies $d_{ij} < r_i + r_j$, $j$ is marked as a dead (necrotic) hepatocyte (neutrophils induce the necrosis of hepatocytes during acute liver injury, Ramaiah and Jaeschke<sup>57</sup>). The lifetime of a neutrophil in the system is set to be 48 hours.<sup>60</sup> Neutrophils are added into the lobule system with uniform rate during the first 24 hours. They migrate very quick. In the model, the initial position of each neutrophil is randomly sampled inside the lobule but outside of the lesion (164 µm away from the central vein).

## Behaviors of platelets/sinusoids

Following liver injury, platelets are recruited to the liver and adhere to the endothelium to generate factors such as PDGF and HGF.<sup>36,38</sup> In our model, platelets are not explicitly modeled. The sinusoid spheres within the lesion (<164 µm away from the central vein) are considered as endothelium adherent for platelets, which are an additional source of PDGF. The presence of platelet-adherent sinusoid spheres lasts until
day 2 (the count of platelets resumes to normal after day 2, as reported in Stravitz et al. ¹⁰⁾. The promotive effect of platelet released HGF on the proliferation of hepatocytes is not modeled explicitly. The proliferation of hepatocytes is modeled by applying a spatio-temporal proliferation pattern extracted from experimental data on healthy hepatocytes in the liver according to its location and time. ²⁸

Dynamic of signals

In our model, DAMPs are synthesized by injured and dead hepatocytes; TGFβ is synthesized by activated Kupffer cells; CCL2 is synthesized by activated Kupffer cells, activated HSCs, activated infiltrating macrophages, and platelet-adherent sinusoids; CXCL1 is synthesized by activated Kupffer cells and activated HSCs; PDGF is synthesized by platelet-adherent sinusoids and Kupffer cells. There is also a factor synthesized by activated infiltrating macrophages. We assume that this factor takes the same diffusion coefficient and decay rate as TGFβ. The simulated sensitivity test of the concentration of each signal to regulate the behavior of a certain cell type is shown in Figure S3. The production rate coefficient for TGFβ is multiplied by 5 when its source cell, Kupffer cell is phagocytosing dead hepatocytes, by taking into account the observation that the TGFβ expression level in Kupffer cells incubated with apoptotic cells is 5-fold higher than those incubated without apoptotic cells. ⁴⁹

QUANTIFICATION AND STATISTICAL ANALYSIS

All analyses were performed in the software R using custom code (code is available in https://zenodo.org/record/8252951). The statistics were reported as mean values and standard deviation across 4 simulation runs in each case.