A computational model of liver tissue damage and repair

Drug induced liver injury (DILI) and cell death can result from oxidative stress in hepatocytes. An initial pattern of centrilobular damage in the APAP model of DILI is amplified by communication from stressed cells and immune system activation. While hepatocyte proliferation counters cell loss, high doses are still lethal to the tissue. To understand the progression of disease from the initial damage to tissue recovery or death, we computationally model the competing biological processes of hepatocyte proliferation, necrosis and injury propagation. We parametrize timescales of proliferation (α), conversion of healthy to stressed cells (β) and further sensitization of stressed cells towards necrotic pathways (γ) and model them on a Cellular Automaton (CA) based grid of lattice sites. 1D simulations show that a small α/β (fast proliferation), combined with a large γ/β (slow death) have the lowest probabilities of tissue survival. At large α/β, tissue fate can be described by a critical γ/β* ratio alone; this value is dependent on the initial amount of damage and proportional to the tissue size N. Additionally, the 1D model predicts a minimum healthy population size below which damage is irreversible. Finally, we compare 1D and 2D phase spaces and discuss outcomes of bistability where either survival or death is possible, and of coexistence where simulated tissue never completely recovers or dies but persists as a mixture of healthy, stressed and necrotic cells. In conclusion, our model sheds light on the evolution of tissue damage or recovery and predicts potential for divergent fates given different rates of proliferation, necrosis, and injury propagation.


Introduction
The liver performs a variety of essential functions including the clearance and metabolism of toxins in the bloodstream. It is generally considered to consist of repeated units called lobules, with a human liver containing a billion lobules. Individual lobules are "plumbed" in parallel a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

APAP induced liver injury
The majority of a therapeutic dose of APAP is metabolized by glucuronidation or sulfation and rapidly excreted in urine. The remainder is converted by CYP2E1 to the reactive compound N-acetyl-p-benzoquinone imine (NAPQI) that readily reacts with sulfhydryl groups on proteins. At therapeutic APAP dosages the small amount of NAPQI produced is detoxified by binding to the cysteine thiol of glutathione (GSH) and excreted from the liver via bile. With higher doses of APAP the glucuronidation and sulfation pathways become saturated and a larger proportion of the APAP is metabolized by CYP2E1 to NAPQI. Cellular GSH then becomes depleted and the excess NAPQI is free to form harmful adducts with other proteins within hepatocytes [5,10]. The most harmful effects result from NAPQI binding to mitochondrial proteins [11][12][13].
Intracellular events due to APAP induced hepatotoxicity, as outlined by Jaeschke et al [14] can be characterized by a phase of initiation in which APAP overdose results in NAPQI formation and GSH depletion [5,15,16], followed by APAP protein adduct formation, mitochondrial stress and compromised cellular respiration [10]. This is followed by the amplification of the injury in which oxidant stress results in mitochondrial DNA damage and loss of mitochondrial membrane integrity, collapsing the membrane potential and disabling ATP production. Without cellular energy production, membrane integrity is lost and nuclear DNA fragments, resulting in necrotic cell death.
The initial APAP injury is pericentral and coincides with the pericentral localization of high CYP2E1 expression. This damage is characterized by the early appearance of stressed cells near the central vein (CV) [15,17]. At a sufficiently high APAP dose, the injured region expands towards the periportal region [18] and can result in complete liver failure.
Some researchers have attributed this progression of damage to intercellular communication mediated by gap junctions. Experiments with mice deficient in connexin 32, a key gap junction protein, given a normally toxic dose of APAP, showed reduced necrosis and that necrosis did not propagate throughout the lobule [19]. Gap junctions appear to amplify the damage by propagating oxidative stress signals between adjacent hepatocytes. Other experiments have also seen synchronized deaths in hepatocytes mediated by gap junctions [20].
In response to APAP induced injury, hepatocytes attempt to compensate for the loss in liver mass by proliferating [21]. In general, this regeneration is related to cytokines produced by the innate immune system of the liver.
Immune system components are activated as necrotic hepatocytes release components, including nuclear DNA fragments, formyl peptides and HMGB1 (High-Mobility Group Box-1 Protein), that act as damage-associated molecular patterns (DAMPs) [14,22]. HMGB1 proteins bind to the liver resident macrophages, called Kupffer cells, through toll-like receptors (TLRS) [23], inducing further cytokine secretion that recruits neutrophils and Ly6C high monocytes. These monocytes and Kupffers secrete pro-inflammatory cytokines such as TNF-α, and IL-6 that also have pro-mitotic effects by priming healthy hepatocytes to be more responsive to growth factors [21,[24][25][26][27]. At the same time these cytokines also activate death pathways in cells [28,29]. The determining factor of whether a hepatocyte proliferates or undergo necrosis or apoptosis may depend on the amount of intracellular ATP [30].
These immune system components have the potential to release additional reactive oxygen species causing additional damage [31]. Similar behavior can be attributed to neutrophils that come in to clear the dead cells, and are also capable of releasing a variety of oxidants [32,33]. While neutrophils don't normally damage healthy cells they can make mistakes and kill normal and stressed cells in the vicinity of damage [33,34]. Thus, although the sterile inflammatory response is necessary to remove cellular debris and activate liver regeneration, this response also has the potential to aggravate the injury.
These observations raise the question; what tips the balance where the same set of cell behaviors that are needed for tissue repair and survival can, in some cases, lead to widespread cell death and irreparable tissue damage? To develop an understanding of the progression of liver damage in this system we have developed a Cellular Automata (CA) model of hepatocyte injury propagation, death and proliferation using a one dimensional (1D) linear chain and two dimensional (2D) hexagonal grids of simulated hepatocytes with key parameters associated to the timescales of these three processes. CA based models have been previously used to model a wide variety of processes including stock market dynamics [35], spread of forest fires [36,37], cancer progression [38], proliferation in migratory cells [39], invasion [40], neurosphere growth [41]. These models have been shown to be capable of producing rich parameter sets of outcomes and have been additionally used to describe phase transitions of infection dynamics [42,43] and traffic flow [44].

Approach
Our CA model consists of hepatocytes as discrete sites on either a one dimensional (1D) chain of cells, or as a two dimensional (2D) grid on a hexagonal lattice. Both 1D and 2D simulations represent typical representations of a liver lobule.
Our 1D model consists of a string of hexagonal hepatocytes representing a simple hepatocyte cord running from a portal triad to the central vein (PT to CV). A 2D simulation consists of a regular hexagonal array of cells representing a hexagonal lobule, centered on a CV, consisting of hexagonal hepatocytes. The 1D and 2D arrangements allow us to explore the spatial effects on the different processes based on the positional inputs and number of neighbor cells.
A cell lattice site can be a cell in a healthy (H) state, a cell in a stressed (S) state, or an empty site where a cell has died (D). Similar to what has been implemented in [45], a healthy cell can be converted to a stressed one due to the presence of other stressed cells in the neighborhood.
A stressed cell eventually dies, leaving behind an empty space. However, recent intravital observations have also revealed that stressed hepatocytes can recover from APAP induced injury through mitochondrial repolarization [16]. Our current model does not include this process but it could be included as a transition from a stressed cell back to a normal cell. Healthy cells can proliferate and repopulate neighboring empty sites, which couples cell division to the death of a neighboring cell.
Our model has three different timescales -(A) α, associated with the cell proliferation timescale, (B) β, associated with the conversion of a healthy cell to a stressed cell, and (C) γ, associated with the process of the stressed cell dying. Varying these parameters allows us to observe the phase space of outcomes leading to either tissue recovery or complete tissue death. We measure the average number of healthy cell states as the characteristic metric at the end of each simulation.
The CA model uses stochastic transitions based on our parameters of α, β, and γ and the state transition rules. We initialize the CA (Fig 1A) to represent the experimentally observed first appearance of pericentral stressed cells at around 2 hours (see Fig 2H) and assume an initial region, defined by a distance from the CV, consisting completely of stressed cells. NAPQI protein adduct studies at different doses of APAP have resulted in different centrilobular of updating the model with α as the proliferation timescale, β as the timescale of conversion of a healthy to a stressed cell and γ as the timescale of stressed cell death. The cell to which the rule is being applied is marked by an X. The Conversion Rule converts a cell from normal to stressed and is dependent on the number of stressed neighbors around the healthy cell. Two random numbers are chosenr 1 has to be lesser than the conversion parameter and r 2 has to be lesser than the total number of stressed neighbors (N s ) normalized over the total number of the cells' neighbors (T). The Death Rule kills a stressed cell if a random number is less than the death parameter P d . The Proliferation Rule fills in an empty location (created by a dead cell) and will only occur if the space is next to a healthy cell. For each healthy cell, a random number r p is chosen and proliferation occurs at that time step if it is less than the proliferation parameter P p and one of its six immediate neighbors contains an empty site (N d >1).
https://doi.org/10.1371/journal.pone.0243451.g001 staining of mouse liver 24 hours after APAP at 0, 10, 100, 250 and 500 mg/kg. Bar = 300 μm. (F) Serum AST and ALT 24 hours after APAP at 0, 10, 100, 250 and 500 mg/kg (n = 3 animals in each group). No necrosis is evident at 10, or 100 mg/kg APAP although modest increases in AST/ALT were observed. Pericentral necrosis, visible as regions of hepatocytes with lighter pink cytoplasm and condensation and loss of hepatic nuclei, has developed at 250 and 500 mg/ kg and is accompanied by major increases in AST/ALT (mean ± s.e.). Time Course: (G-J, D, K) H&E staining of mouse liver at 0, 2, 6, 12, 24 and 48 hours after APAP at 250 mg/kg. (G-K) Bar = 300 μm. (L) Serum AST and ALT at 0, 2, 6, 12 and 48 hours after APAP at 250 mg/kg (n = 3 animals in each group). Data re-visualized as in Dunn et al [16]. At 250 mg/kg the cytoplasm of cells around the CVs begins to appear lighter by H&E staining as early as 30 minutes after administration of APAP. The lightened region becomes easily discernible around 2 hours after APAP treatment, patterns by 2 hours, with higher doses showing the adduct levels across a larger area [18]. We implement this dose effect as a shift in the initial position of the boundary between healthy and stressed (or dead) cells towards the portal region, resulting in an increasing number of stressed cells and decreasing number of healthy cells.

Rules for updating the CA
At any point in time there are N total sites made up of healthy cells, stressed cells and dead/ empty sites (Fig 1B). For each cell we keep track of how many of its neighbors sites are stressed (N S ), healthy (N H ), or empty (dead) (N D ). The automaton evolves through the interaction of these states through equally spaced time points with a time step of Δt. Each of these cells can choose to update their state according to the three rules described below at every time point.
(A) Proliferation. Hepatocytes in the liver try to proliferate to maintain a constant mass. In our model only healthy cells can proliferate. We assume a characteristic time scale of α associated with the proliferation process, such that the probability of a hepatocyte dividing in a computation time step becomes P p = Δt/α. Since this is a probability, we require that 0<Δt�α. All the other parameters associated with timescales will also obey this relation.
Additionally, we assume that this process is limited by space and that the ability of a cell to divide is dependent on its immediate cell density. This is based on the observation that division occurred adjacent to necrotic regions and similar assumptions have been used in [40,41], based on cells' ability to mechanically sense its neighborhood [46]. Contact-inhibited growth has also been seen in cultured hepatocytes in [47].
Computationally to achieve both of these criteria, we begin by picking a random number r p from a uniform distribution for each healthy hepatocyte and checking if r p �P p . If yes, the healthy cell is added to a queue. Next, for all the cells in the proliferation queue we check if each hepatocyte has free space to move into by checking for the presence of empty spaces left behind by dead cells. In the 1D case, each healthy hepatocyte can look at its two nearest neighbors (to its left and right) while in the 2D case hepatocytes can look at six of its nearest neighbors. We make a binary choice here, if any of its neighboring sites are empty the cell will choose to divide and if not, the cell won't. Upon division the hepatocyte will take the place of the empty site. If multiple cells pick the same empty site for division, the cells are randomly shuffled and one cell is picked.
(B) Cell death. Stressed cells undergo necrosis by losing plasma membrane integrity. We assume that there is a characteristic time scale associated with this process of stressed cells dying, denoted by γ. We assume that this process is independent of neighbor density and any stressed cell can die with a probability of P d = Δt/γ. Similar to choosing cells for proliferation, we pick a random number r d for each dead cell and check if r d �P d .
(C) Conversion of a healthy to a stressed cell. We assume that stressed cells can communicate and trigger the conversion of a healthy to a stressed cell. We assume that there is a characteristic time scale β associated with this process and the process is dependent on the presence of other stressed cells in the healthy cell's neighborhood.
We pick two different random numbers for this process -(a) A random number r 1 is first chosen for each healthy hepatocyte to see if its lesser than the conversion parameter P c = Δt/β. If the condition is satisfied, the cell is added to a queue, (b) A second number r 2 is chosen and compared to the total number of stressed neighbors (N s ) in its neighborhood. A cell is but hepatic nuclei remain intact and AST/ALT is unchanged from control. Pericentral necrosis has begun by 4 hours after APAP treatment and is accompanied by increased AST/ALT. Pericentral necrosis is prominent at 6, 12 and 24 hours and is accompanied by further increased AST/ALT. (K) By 48 hours pericentral region shows signs of recovery. https://doi.org/10.1371/journal.pone.0243451.g002 converted to a stressed type if r 2 �N s /T, where T is the total number of neighbors. Previously, majority rules, where the fate of each cell depends on the state of the majority of the neighbors have been also used to model opinion dynamics [48,49]. For 1st order neighbors, we consider 2 nearest cells in 1D and 6 nearest cells in 2D hexagonal lattice.
A 2D schematic of the rules is listed in Fig 1D-1F. At the end of our simulations, we classify any case that recovers all it's healthy population (H = N, S = 0, D = 0) as recovery of the tissue and any simulation that results in all dead states as one of tissue death (H = 0, S = 0, D = N).

(1) 1D pericentral necrotic damage and recovery
We first examined hematoxylin-eosin (H&E) stained liver tissue from mice sacrificed 24 hours after IP administration of APAP at 0, 10, 100, 250 and 500 mg/kg for evidence of hepatotoxicity. Necrotic zones are characterized by cytoplasmic vacuolization, loss of hepatic cytoplasm and nuclei, and congestion. In control sections and at low doses of APAP no discernable necrosis was present (Fig 2A-2C). At higher doses (Fig 2D and 2E) widespread centrilobular necrosis was evident. We also measured serum levels of the liver enzymes, alanine aminotransferase (ALT) and aspartate aminotransferase (AST), which are common clinical biomarkers of liver injury. Leakage of ALT/AST into the bloodstream occurs following hepatocyte injury and high levels can be indicative of significant necrosis. At 10 mg/kg and 100 mg/kg we saw modest increases in both ALT and AST ( Fig 2F) even though no necrosis was evident by H&E (Fig 2B  and 2C). At 250 and 500 mg/kg APAP we saw striking increases in both ALT and AST ( Fig  2F), accompanying the widespread centrilobular necrosis seen by H&E (Fig 2D and 2E). This is in agreement with a detailed study by Roberts et al. (18) where they documented a doserelated increase in extent and area of centrilobular necrosis, accompanied by increases in centrilobular APAP adduct localization, and of GSH and ALT in serum. Further, by 48 hours mice recover from 250 mg/kg APAP, while 500 mg/kg is lethal. The questions that this suggests are: Given this damage pattern, what cell behaviors might underlie recovery in one case and death in the other; and what might be the outcomes of changes in extent of necrosis?
We then examined the time course of hepatocyte damage and recovery at 250 mg/kg APAP by H&E (Fig 2G-2K). Signs of stressed cells are seen by 2 hours after APAP, characterized by the appearance of hepatocytes with lightened cytoplasm around the CV ( Fig 2H). Nuclei were still present in these stressed cells and appeared normal in size and morphology. Analysis of the initial damage pattern at 2 hours showed that the stressed area extended 32.6% of the distance along the CV to PT axis with a standard deviation of 16%. ALT/AST levels unchanged from control at this point indicating that hepatocytes were still intact ( Fig 2L). Necrosis was progressive over time. By 4 hours necrosis was visible as lightened cytoplasm and condensation and loss of hepatic nuclei and was accompanied by congestion of the sinusoids with red blood cells (RBCs) around the periphery of the necrotic region. Necrosis was more prominent by 6 hours, with extensive loss of hepatic nuclei and continued RBC congestion. At both 4 and 6 hours after APAP treatment, necrosis was accompanied by an increase in ALT/AST, with the peak value occurring by 6 hours. The ALT/AST levels remained elevated and hepatocytes around the CV still appeared damaged at 24 hours. By 48 hours ALT/AST levels were near normal and the tissue had recovered with centrilobular hepatocytes appearing similar to the unaffected periportal cells. Taken together, the H&E data (Fig 2) gives us spatiality and a time scale for when cells become stressed (2 hours after APAP), have undergone necrosis (6 hours after APAP) and when the damaged centrilobular region is on the path to recovery (48 hours after APAP). [31,50]. Kupffer cell activation is accompanied by their release of pro-inflammatory cytokines that are largely cytoprotective but can also mediate cell damage [31]. We quantified the levels and localization of resident and recruited Kupffer cells after APAP treatment, using F4/80 labeling (Fig 3A-3E). We analyzed multiple lobular regions (n) around the CVs for a single animal; n = 9 (0 hours), 14 (1 hours), 7 (2 hours), 7 (4 hours), 6 (6 hours), 5 (12 hours), 4 (48 hours) and 12 (72 hours). In untreated mice there were very few Kupffers and those that were present were rarely found in the centrilobular region ( Fig 3A, 3D and 3E). Numbers and distribution of Kupffer cells remained unchanged from control through 4 hours after APAP treatment (Fig 3D and 3E). By 6 hours after APAP, Kupffer cell numbers had increased and they were localized primarily within the fourth cell layer out from the CV, at the periphery of the damaged centrilobular area, at the boundary between healthy and necrotic tissue (Fig 3B, 3D and 3E). This coincides with the worsened necrosis and peak ALT/AST levels also seen at 6 hours after APAP treatment (Fig 2). Proximity of the Kupffer cells to the edge between healthy and necrotic tissue suggests that any healthy cell at this location is at a greater risk of being affected by the damage. Kupffer cell numbers increase over the next 42 hours and their localization progresses pericentrally as the necrotic tissue is repaired (Fig 3C-3E), evidenced by the reduced ALT/AST levels at 48 hours after treatment. By 72 hours after APAP the Kupffer cells have returned to control levels and distribution ( Fig 3E).

Resident and recruited Kupffer cells are activated by 2 hours after an APAP overdose
We examined changes in cell proliferation levels and localization using phosphohistone H3 labeling (Fig 3F-3J). We performed our analysis for a single animal at each of 6 time points with multiple (n) regions sampled; n = 10 (0 hours), 9 (1 hours), 9 (2 hours), 8 (4 hours), 8 (6 hours) and 8 (48 hours). In control liver sections dividing cells are sparse and not concentrated in either periportal or pericentral regions (Fig 3F, 3I and 3J). At 1, 2, and 4 hours after APAP proliferation levels and localization are unchanged from control ( Fig 3I and 3J). By 6 hours after APAP increased proliferation is evident and is concentrated at the periphery of the centrilobular region of necrotic cells (Fig 3G, 3I and 3J). The localization of dividing cells at 6 and 48 hours after APAP correlates with the leading edge of tissue recovery as it progresses toward the CV (Fig 3G-3J). This localization of hepatocyte proliferation during recovery of centrilobular necrosis being predominantly at the margin of the necrotic region supports our model assumption that loss of contact inhibition contributes to the increase in hepatocyte proliferation (Figs 2 and 3).
Based on our experimental results with APAP, we next show results for a CA model with initial conditions around 2 hours, characterized by the first appearance of stressed cells (Fig 2). On assuming a hepatocyte width of 30 μm and dividing the average CV to PT distance, we estimate 10 hepatocytes to be present in this region. We run our base model by setting N = 10 with 40% of these hepatocytes to be stressed, i.e, H i = 6, S i = 4, D i = 0, since the H&E slides show approximately 30-40% of the CV to PT region as stressed. The number of healthy or stressed cells and the size of the tissue are all parameters in our model and we also make observations on varying them. We use a timestep of Δt = 0.25 (hours) for the simulations and unless otherwise stated, all of our results are performed over 100 trial replicas. Our results explore the parameter space in terms of ratios of the timescales, by fixing any one timescale to a value and incrementing the other two timescales in steps of 0.50 (hours). This allows us to make general comments based on the output of the simulations in terms of biology of the system.

(2) Small α combined with a large γ leads to low recovery probabilities
In the 1D model, we identify two steady states for our simulations. H = N, S = 0, D = 0 is reached when the healthy states cells grow and the cells all completely fill up the simulated . We find that the probability of the tissue surviving decreases with increasing ratios of γ/β. This can be understood as follows -if γ is small compared to β, stressed cells die rapidly before influencing the adjacent healthy cells. As γ increases, the long life of the stressed cells leads to greater damage to healthy cells. However this increase in γ does not lead to a sharp transition to the steady state H = 0, S = 0, D = N, but the stochastic model instead shows regions of bistability in which the tissue could either survive or die. These bistable regions also show a gradient of output with the tissue recovery probability steadily decreasing as γ increases (Fig 4A and 4B).
We define a criteria to mark the transition into the bistable region by tracking γ/β and α/β values in which the probabilities of survival fall from 1 to 0.95 (points colored in red in Fig 4A). We find these points by finding the first occurrence of γ/β for every α/β for which the survival probability is close to 0.95±0.01. The trend followed by the points indicate that this transition in 1D is almost independent of the proliferation timescale α, at large values of α, and can be approximately described by a critical ratio of γ/β � = 2.25±0.37 (red line in Fig 4A). Sample outputs showing some of the trajectories along α/β = 1 for different γ/β ratios are shown in S1 Fig. For smaller timescales of proliferation, the parameter space shows variation in the survival probabilities as a function of γ/β and this behavior is more prominent beyond our critical ratio of γ/β � . Fig 4B shows the survival probabilities as a function of γ/β for α/β = 0.25, 3.25, 5.25 shown in blue, black and magenta respectively. Interestingly, we find that that for the same value of γ/β, the probability of the tissue surviving is actually lowest for a small (fast) timescale the CV (within 1 cell diameter from the CV). (J) Dividing cells are concentrated at the periphery of the necrotic region at 6 hours after APAP. By 48 hours after APAP dividing cells are concentrated near the CV.
https://doi.org/10.1371/journal.pone.0243451.g003 of proliferation. We interpret this as follows -in this region (γ/β>γ/β � ) the conversion from a healthy to a stressed cell is more rapid than the clearance of the stressed cells. Cells will proliferate till they are limited by the number of spaces and the slow clearance of stressed cells. If proliferation is fast, healthy cells rapidly provide a domain for the stressed pattern to propagate into, as opposed to the inability of stressed cells to propagate into stressed or dead regions. This can also be seen in the S2 Fig that shows two different trajectory trials for a small proliferation timescale α/β = 0.05 (S2A and S2B Fig) and a larger proliferation timescale of α/β = 1.05 (S2C and S2D Fig). Both are for a large (slow) death timescale with γ/β = 8.05. For the smaller proliferation timescale, we see the healthy population, shown in blue, peaking first, indicating rapid proliferation, that then begins to decay as the stressed cells take over and kill the healthy population. On the other hand, the few trajectories that survive for α/β = 1.05 don't usually show signs of increased proliferation till after the stressed population has decayed. Fig 4C also shows this behavior for three different values of γ/β = 0.25, 3.25, 5.25 (blue, black, and magenta respectively) indicating that this behavior is only amplified for larger γ; the magenta curve shows a lower survival probability for small α before saturating to a probability of about 0.6.
We also plot the parameter space for β vs γ by fixing α = 5 (hours) ( Fig 4D). As expected from Fig 4A, the transition to the bistable region is approximately described by a linear function with a constant ratio of γ/β � . By averaging over γ/β ratios where the survival probability of the tissue falls to 0.95, we obtain γ/β � = 1.69±0.24 (close to the ratio obtained in the previous case), indicated by the red line. For very low γ, we see that the probability of tissue survival is very high for all values of β as the stressed cells die rapidly before being able to affect other cells. As γ increases the bistability region expands leading to fewer survival chances for the tissue as a function of β. Fig 4E shows  In between these two cases of small and large γ, we see a peak in the number of new hepatocytes created. This peak is more prominent at smaller proliferation timescales. Additionally, since smaller proliferation timescales also result in a greater probability of tissue death, we see a decrease in the number of new hepatocytes created for small α and large γ (Fig 5A and 5C). Thus, smaller proliferation timescales produce a larger peak indicating a large regenerative capacity for moderate values of γ, before declining more for large values of γ, as compared to large values of α. Additionally, the line in Fig 5A indicates γ/β � = 2.25±0.37, denoting the ratio of γ/β above which the survival probability falls below 0.95. We note that the peak of maximum divisions is slightly to the right of this line, suggesting that in some cases even having the ability to regenerate many times isn't still sufficient to stop the tissue from dying. Simulations with larger N show the same behavior with now the γ/β � and the peak of maximum proliferation shifted to the right (S4 Fig). Fixing α = 5 (hours) shows a region characterized by γ/β ratios for which maximum proliferation can be now seen (Fig 5B). Horizontal lines across this plane are shown for fixed γ/α = 0.15 (blue), 3.25 (black) and 5.25 (magenta) (Fig 5D).

(4) 1D critical threshold of damage pattern propagation
Our in vivo mouse dose response studies indicated that APAP at high doses (500 mg/kg) is often fatal by 48 hours, presumably due to massive liver necrosis, while animals treated with a lower dose of 250 mg/kg recover. Here we ask the question -Is there an extent of tissue damage beyond which damage is fatal and unrecoverable? To do this, we look at the minimum number of healthy cells min(H) that simulations can reach from an initial number of healthy cells (H i ) and still eventually survive. We quantify this loss of healthy cells as a function of the parameters with We make our observations with a fixed α and quantify our results in terms of γ/β ratios. Fig 6A shows how the tolerable extent of damage is dependent on γ/β (colorbar shows the loss). For tissues that have a survival probability of 0.95 and above, the maximum damage is seen to be around 20% as γ/β approaches 2.0. As γ/β is increased even further, much larger extents of damage can cause lower survival probabilities. However, parameters set with a very high ratio of γ/β showed reduced signs of damage for simulations that progressed to survival of the tissue. This is because under these parameter combinations, the simulations that do survive must only sustain minimal losses from the initial number of healthy cells.
Alternatively  across the entire studied parameter range is below 10%. The gain for high γ/β is low, as damage is easily propagated due to the longer presence of stressed cells. The gain in cells increases as γ/ β decreases, before falling again as γ/β � is reached where stressed cells are now more easily eliminated. There are very few simulations that progress to complete cell death in this region. A few outliers are seen, characterizing lone trials where this gain is as large as 50%.
Next, while this loss or gain percentages are parameter dependent, we also checked if the tissue shows a critical size of healthy populations. We do this by plotting hmin(H)i as a function of the binned final average tissue length, hHi/N. Since this final length also represents the probability of the tissue surviving, this will give us a correlation on how likely it is for the simulation to be successful after reaching a certain minimum healthy cell population. We find that with N = 10, and H i = 6, S i = 4 (initial damage of 40%), the average minimum length in simulations that survive saturate to around 4 healthy cells, indicating a lower threshold for parameter combination with large γ/β ratios (Fig 6B). Additionally, we find that this lower threshold is dependent on the tissue size (N). Varying the size of the tissue with different N and the same damage fraction (around 40%) shows different critical lengths for the minimum number of healthy hepatocytes. However, scaling these trends by N shows that the critical length is always the same for the same initial amount of damage fraction (Fig 6C). In this case, this fraction is 0.37±0.04, indicating that if the healthy population drops below about 37% of the total population, the damage is likely to progress to complete cell death. However, this critical length is dependent on the initial amount of damage as shown in Fig 6D. Here the different lines are

(5) Critical ratio γ/β � is proportional to the size of the tissue
The previous simulations indicated that a minimum threshold is dependent on the size of tissue. We additionally quantify how the phase spaces behave as we vary the number of initial stressed cells and the total size of the tissue. On changing H i = 2, S i = 8 (initial damage at 80%) we find that the region where the tissue survives is reduced drastically (Fig 7A). Approximating the boundary where the tissue survival (and the average length) falls to 0.95±0.01, we measure a slope (γ/β � ) of 0.38±0.11. This slope value increases to 1.07±0.18 for H i = 4, S i = 6 (initial damage at 60%) ( Fig  7B) and 2.53±0.35 for H i = 8, S i = 2 (initial damage at 20%) (Fig 7C). These results again show how the probability of recovery is dependent on the initial length of damage in the 1D model.
Similarly, we can vary the total length of the tissue (N), simulating the highly variable CV-PT distance observed in liver lobules, and evaluate how the initial damage pattern now affects the survival probability of the tissue. We found that for a given fraction of initial damage, a larger tissue size expands the parameter region where the tissue survives (Fig 7D-7F).
We can see the combined trends in Fig 7G, with the critical ratio γ/β � plotted on the y axis. We find that rescaling all the curves by the length of the tissue essentially collapses them ( Fig  7H), showing that γ/β � is dependent on the tissue size. For 40% initial damage, the approximate value of γ/β � = (0.17±0.02)N.  of stressed cells centered at the CV to represent 40% damage along the CV-PT axis. We run our parameter scans for a fixed β = 5 (hours). Additional examples of exploring the phase space for different values of α is shown in S7 and S8 Figs. We find that now we have a richer phase space comprising solutions in which the tissue fully recovers (H = N, S = 0, D = 0), where the tissue completely dies (H = 0, S = 0, D = N) and a "coexistence state" where neither the healthy nor the dead population take over but instead the three cell types coexist indefinitely (H � , S � , D � ). This third state appears to be caused partly by cells having a larger number of contact neighbors (six in 2D versus two in 1D) which provide additional information to the cell, and partly from parameters being in regions of maximum regenerative capacity (Fig 5). Fig 8A shows the phase space of all the outcomes from the trials. In Fig 8A, any parameter combination over all the trials that leads to the complete recovery of healthy cells is colored yellow (marked as S), simulations where all the cells die are colored green (marked as D) and simulations where coexisting populations are obtained are shown in red (marked as C). We note that the transitions between these different phase regions are not sharp, but rather gradual as indicated in Fig 8A, apparently because of the stochasticity of the model. Unlike the 1D case, the 2D case does not predict a single γ/β � ratio for the transition from recovery to the different phase outcomes. However, we note that signs of tissue death in Fig 8A begin   cells die rapidly and healthy cells quickly take their place. As γ/β is increased to 1.45 (Fig 8D  and S6B Fig), stressed cells live longer affecting the tissue to a greater degree, but the tissue can still recover. On increasing γ/β = 2.45, we reach a domain where stressed cells and healthy cells coexist without any of the populations ever winning out (Fig 8E and S6C Fig). The average number of healthy or stressed cells depends on the value of γ/β ratio and as this value is increased to 4.45 (Fig 8F and S6D Fig) As γ/β is increased even further to 9.05 ( Fig 8G and S6E Fig), the stressed cells persist longer and are able to affect other cells. This eventually pushes the tissue towards death. Fig 8H focuses on the healthy steady state populations across different values of γ/β for fixed α/β = 1.65(blue), 5.95 (black) and 10.25 (magenta). The first two lines cross the coexistence regions and a dip in the population numbers can be seen in the black line. This is interesting as it shows that coexisting states with a lower healthy population exist next to parameter regions showing a high degree of survival, before transitioning into bistable regions again. Additionally, as in the 1D case, fast proliferation appears to lead to a decreased healthy state population.
We also find that any initial asymmetries in the distribution of stressed cells doesn't change the nature of the phase space outcome (S13 Fig). However, simulations with a larger number of cells (N L = 50, N = 7650) showed expansion in the coexisting state space regime (S14A Fig). Varying the initial conditions in this case showed shifts in the SD boundary as in the 1D case with little changes in the localization of the coexisting states (S14B and S14C Fig).

(7) 2D Space with second order neighbor effects
We explored the parameter and phase space of the 2D model in the presence of second order neighbor effects. In this case, a given healthy cell responds based on the state of both its six first order neighbors and on its 12 second order neighbors. This simulates the effect of a short range diffusible signal that increases the range of cell-cell communication from contact neighbors to neighbors separated by one intervening cell. We find that while the nature of the 2D phase space doesn't change greatly (Fig 9A), a shift in the boundary between the survival (S) and the bistable (SD) regions can be seen (compare Figs 8A and 9A). Additionally, at large γ/β, for the same parameter values as in Fig 8A, the tissue shows a small increase in the probability of tissue death. Fig 9B plots

(8) 2D Coexisting states depend on the number of divisions
In the previous 1D and 2D simulations, each hepatocyte can divide as many times as needed. Actual cells may have a limited division potential. Hepatocytes are known to divide throughout the life of an individual, although some experiments suggest that the number of divisions per hepatocyte may be limited [51,52]. We find that fixing the number of maximum allowed divisions per hepatocyte to 1 removes the coexistence states ( Fig 9D). However, including this limitation also limits the hepatocytes' ability to completely repopulate the liver lobule space. Therefore, we relax our criteria to denote our survival 'S' states as any simulation outcome that has healthy cells at the end without any stressed cells. Fig 9D indicates that we begin seeing signs of complete tissue death near γ/β = 5.25±0.29, this is near the same ratio of γ/β = 5.04 ±0.30 where we saw the first appearance of complete tissue death in simulations with no limit on the maximum number of divisions. Fig 9E shows the average healthy population as a function of γ/β and α/β. Fig 9F shows

Discussion
In this paper we have explored a CA based model of hepatocyte behavior in the liver following a drug induced centrilobular injury. Based on observations from our experimental data in mice, we have represented hepatocytes in three different states -healthy, stressed and necrotic/dead. The simulated cells change states based on their current and neighboring cells' states. Healthy cells can divide with a proliferation time scale of α to regain liver mass by replacing adjacent empty spaces left by dead cells. Healthy cells can also be converted into stressed cells due to interaction with neighboring stressed cells with a timescale of β. The stressed cells disappear (die) after becoming necrotic with a timescale of γ, leaving behind empty spaces. We model the process of response to an initial large scale tissue damage as a competition between these different processes.
To facilitate exploring the parameter space, we vary two of the parameters while keeping the third fixed and expressing the varied parameters as a ratio with the fixed parameter. By varying the ratios of any two timescales, we discover parameter combinations that lead to either tissue recovery or death. We find that the transition from recovery to death is not sharp but characterized by parameter sets where either outcome is possible thus showing probability gradients of tissue recovery. For larger values of the proliferation timescale α, the probability of the tissue surviving can be described by the γ/β ratios alone ( Fig 4A). However, at smaller values of α, variations in α become important with the probability of tissue death being the highest for smaller α/β ratios and large γ/β ratios (Fig 4A and S2 Fig). We can characterize the parameters that show tissue death by approximating critical ratios of γ/β � beyond which the tissue has less than a 0.95 probability of survival. This γ/β � is dependent on the size of the tissue (N) and initial amount of damage (number of stressed cells). For an initial 40% damage, γ/β � � (0.17±0.02)N (Fig 7). Our CA model also predicts a minimum healthy population size of 0.37N for the initial 40% damage, below which the liver damage is likely irreversible (Fig 6).
Applying the same set of rules to a 2D model produces a richer parameter and phase space where additional coexisting states can be identified (Fig 8A). In the coexisting state, the cells persist in each of the three possible states without the system ever going to complete recovery or complete death. The coexisting states are stable to external perturbations that give rapid changes in the size of one of the three populations (S9-S11 Figs). These coexisting states are localized around regions of maximum proliferation ( Fig 5) and limiting the number of divisions per hepatocyte to 1 removes these coexisting states (Fig 9). However, doing so also limits the hepatocyte's ability to completely fill the liver tissue. This suggests that other requirements need to be put in place, like an accompanying increase in the healthy hepatocyte sizes during recovery as seen in experiments [52]. This simple CA model makes interesting predictions on the parameter and phase spaces of hepatocyte behavior in a spatially defined context. Our model can be compared to biological parameter regimes by observing that the first signs of stressed cells appear around 2 hours and prominent necrosis appears over the next 4 hours. Previously reported literature values have also suggested that the hepatocyte doubling time is approximately 24 hours [53]. Based on these observations we can approximate our parameters as a�20 hours, with β�1−2 hours and 3 stressed layers damage along the CV-PT axis). The high probability values at these ratios indicate that the tissue will recover under these conditions. This is also seen in our histology slides. Additionally, S16 Fig showing the normalized damaged fraction at the different time points from H&E slides indicate that the fraction remains relatively unchanged, with no large losses to the healthy cells, as is also seen in these simulations. However, this analogy is limited somewhat by the limited number of time points (four) in the in vivo studies. Future mouse models with a finer time resolution of spatial necrosis patterns could help test these predictions based on the three suggested parameters. Additionally, we can also compare the time estimate of tissue recovery after the APAP treatment by adding an initial 2 hours for time taken to reach the steady state in our simulations. For our 1D simulations, for the chosen sets as above, this estimate is around 25-30 hours while for the 2D simulations it is 22-25 hours. This average value of around a day is lower than recovery period of >48 hours as is seen in the H&E slides. However, since this time is heavily dependent on the value of α, setting the time for how fast the healthy cells proliferate, a better estimate of this value will be able to approximate the model behaviour better. Finally, while we don't expect an additional recovery timescale to change our results qualitatively, we expect that the boundaries of the different domains would shift, which could also lead to finer predictions for the experiments.
While these results have been based on a pipe or a planar layout, these simulations can also be extended to a 3D framework where cells are in contact with additional neighbors in planes above and below making the total number of neighboring sites 8. Compared to the 2D case, these simulations show that the coexistence regions expand indicating that the neighbor interaction range could be an important parameter (S17 Fig). In conclusion, we have developed stochastic CA models of progression of injury in liver lobules following a pericentral toxic challenge. The modeling delineates the parameter space and resultant tissue response phase space including complete tissue recovery, complete tissue death and situations where the tissue appears to settle into a partially damaged state that persists indefinitely, as long as there is sufficient proliferative capacity from existing hepatocytes. This partially damaged persistent state, though probably not relevant to APAP toxicity, may be relevant to chronic liver disease such as hepatitis and alcoholic cirrhosis. Our results also show that rapid proliferation does not always ensure tissue recovery, indicating possible situations where proliferation might not be able to limit liver damage. Biological perturbations to our timescales, aided by a hepatocyte's interaction range, will be able to test out these predictions, including the presence of critical ratios of tissue recovery and the dependence on a minimum healthy hepatocyte population beyond which damage could be irreversible. Regions in the parameter space where it is not possible to predict the outcome of the simulation would be of particular interest as they represent situations where the outcome might be critically dependent on small changes in other factors.

Experimental animal model
All animal experiments were approved and conducted according to Institutional Animal Care and Use Committee guidelines of Indiana University, and adhered to the NRC guide for care and use of animals. Twelve hours before APAP administration, 7-10 week old male C57BL/6 mice were moved to clean cages without food. The day of the study mice received an IP injection of APAP (dose of 250 mg/kg) or saline vehicle (at a volume of 0.2 mL/10 g body weight), and were then returned to cages with food. Mice were then analyzed over time.
Liver enzymes. Liver enzymes, alanine aminotransferase (ALT) and aspartate aminotransferase (AST), were measured immediately following serum separation from whole blood using Infinity™ ALT and AST kits (Thermo scientific, Middletown, VA) with a UV microplate reader (Tecan Infinite M200Pro, San Jose, CA).
Histopathology. Following euthanasia, mouse livers were perfused with warm saline and removed and fixed in 10% neutral-buffered formalin. Fixed liver tissues were paraffin embedded, sectioned (4−5 μm), and stained by hematoxylin and eosin (H&E) for the evaluation of pathological changes.

Fluorescent labeling
Following euthanasia, mouse livers were perfused with warm saline, then with warm 4% paraformaldehyde in PBS. Livers were then removed and incubated overnight in 4% paraformaldehyde in PBS. Fixed livers were vibratome sectioned (200 mm). Liver sections were blocked and permeabilized overnight using 0.5% Tween in Carbo-Free Blocking Solution (Vector Laboratories). Sections were then either incubated in anti-F4/80 antibody (Serotech) to label resident (Kupffer) and infiltrating macrophages [54] or anti-histone-H3 antibody (Invitrogen) to detect proliferation [55], followed by fluorescent secondary antibodies for visualization (Invitrogen). Morphological staining was applied to all sections using rhodamine-conjugated lens culinaris agglutinin (LCA, Vector Laboratories) to highlight cells and tissues and Sytox Green (Invitrogen) to label nuclei. In some sections cell membranes were also labeled using 633-conjugated Tomato Lectin (TL, Vector laboratories). Labeled sections were either mounted in PBS or optically cleared and mounted [56].

Microscopy and image analysis
We acquired confocal images of the fluorescently labeled tissues using Leica SP2 confocal/MP and Leica SP8 confocal microscope systems. We conducted quantitative image analysis using FIJI [57,58]. Using the Cell Counter plugin, we counted F4/80 and phosphohistone H3 positive cells and classified them by distance in number of cell diameters from a central vein. To quantify the proportion of damaged cells in the H&E images, we used Color Deconvolution [59] to separate color channels followed by segmentation using Weka [60] using the intensity of the cytoplasmic staining to identify the transition from healthy to damage cell regions. Along a line from each CV to the nearest PT, we measured the proportion of damaged cells.

Simulations
CA model was implemented using Python. Coexisting states shown in red (marked with C), parameters for which the tissue always survives shown in yellow (marked with S) and regions where the tissue always dies shown in green (marked with D). Similar to the parameter space seen in Fig 8A, for ratios of β/ α below 0.15 (α/β>6.7), the transition from recovery to death goes through regions of bistable outcomes of both tissue survival and death. By increasing that ratio to around β/α = 0.25 (α/β = 4), a mixture of coexistence and survival states is seen which again transitions into regions characterized by death, recovery and coexistence as γ/α is increased. As β/α is further increased (α/β is decreased) pure coexisting states can be seen wedged between different recovery and death regions. in the left column are when 95% of stressed cells are converted into the healthy or dead states while the right column shows output when all but one stressed cell remains. Coexisting states are stable and the perturbations die down even after 95% removal of the stressed cells. (S1 Table).