article: Acetaminophen-induced liver injury in mice is a model for drug-induced liver injury in humans. A precondition for improved strategies to disrupt and/or reverse the damage is a credible explanatory mechanism for how toxicity phenomena emerge and converge to cause hepatic necrosis. The Target Phenomenon in mice is that necrosis begins adjacent to the lobule’s central vein (CV) and progresses outward. An explanatory mechanism remains elusive. Evidence supports that location dependent differences in NAPQI (the reactive metabolite) formation within hepatic lobules (NAPQI zonation) are necessary and sufficient prerequisites to account for that phenomenon. We call that the NZ-mechanism hypothesis. Challenging that hypothesis in mice is infeasible because 1) influential variables cannot be controlled, and 2) it would require sequential intracellular measurements at different lobular locations within the same mouse. Virtual hepatocytes use independently configured periportal-to-CV gradients to exhibit lobule-location dependent behaviors. Employing NZ-mechanism achieved quantitative validation targets for acetaminophen clearance and metabolism but failed to achieve the Target Phenomenon. We posited that, in order to do so, at least one additional feature must exhibit zonation by decreasing in the CV direction. We instantiated and explored two alternatives: 1) a glutathione depletion threshold diminishes in the CV direction; and 2) ability to repair mitochondrial damage diminishes in the CV direction. Inclusion of one or the other feature into NZ-mechanism failed to achieve the Target Phenomenon. However, inclusion of both features enabled successfully achieving the Target Phenomenon. The merged mechanism provides a multilevel, multiscale causal explanation of key temporal features of acetaminophen hepatotoxicity in mice. We discovered that variants of the merged mechanism provide plausible quantitative explanations for the considerable variation in 24-hour necrosis scores among 37 genetically diverse mouse strains following a single toxic acetaminophen dose.
  Acetaminophen (APAP) -induced liver injury (AILI) in mice is the widely used model for drug-induced liver injury in humans. There has been dramatic progress in identifying involvement of a variety of molecular level events and pathways [1–3]. To use that knowledge effectively in predicting injury and developing new treatment strategies, we need more credible explanations of how, when, and where key injury features emerge within hepatic lobules, and why humans exhibit a wide diversity of responses. A characteristic early feature is that necrosis, proceeded by covalent adduct formation, begins perivenous, adjacent to the central vein (CV) of hepatic lobules, and then progresses (radially) outward [4,5]. Given that the fraction of APAP converted to the reactive metabolite (NAPQI: N-acetyl-p-benzoquinone imine) within hepatocytes increases from lobule’s periportal (PP) space to CV (called zonation), the prevailing mechanistic explanations assume that zonation of NAPQI formation is a necessary and sufficient prerequisite to explain early perivenous necrosis. That NAPQI zonation hypothesis provides an explanatory foundation for current efforts to better understand downstream complexities leading to tissue regeneration or liver failure. However, challenging that explanation directly in mice is currently impracticable because it requires exercising precise control over zonation of many mechanism features along with methods to measure intralobular features sequentially within the same mouse. Those barriers have impeded progress. Here, we successfully circumvent those impediments using multi-attribute software-based experiments. We challenge the NAPQI zonation explanation by experimenting on a strongly analogous software mechanism (as distinct from being described mathematically), called NZ-Mechanism, that is instantiated within a virtual mouse. This software analog is engineered to be quantitatively and qualitatively biomimetic across all anatomical, hepatic zonation, and cell biology features currently believed relevant to challenging the NAPQI zonation hypothesis. Results from virtual experiments provide strong quantitative evidence that zonation of NAPQI formation alone is insufficient to explain early perivenous necrosis. A parsimoniously more complex explanatory mechanism is needed, one involving additional feature zonation. After testing several, we focused on two equally plausible mechanism zonation features, but they too failed to meet stringent tests of sufficiency. However, when the two zonation features that had failed were combined into a unified mechanism, sufficiency was achieved, thus demonstrating the potential power of the virtual experiment approach. Physiologically reasonable perturbations of the merged mechanism caused significant changes in both the numbers of hepatocytes killed and the intralobular patterns of dead hepatocytes for a fixed toxic APAP dose of 300 mg/kg. Injury diversity was similar to that reported for humans [6,7] and to differences in 24-hour necrosis scores among 37 genetically different mouse strains [7]. For the latter, we hypothesized that variation of lobule location-dependent features within the combined mechanism, resulting from changing selected mechanism configurations, can be sufficient to account for such necrosis score diversity. We present results of virtual experiments supporting that operating hypothesis with the intent that these results will help formulate new wet-lab experiments to further improve explanatory mechanism insight (Fig 1).
  Mouse Analog is the multiscale biomimetic software system (Fig 2) designed and built intentionally to be scientifically useful [8], especially, as in this work, for testing mechanism hypotheses of AILI, starting with the NAPQI zonation hypothesis: a configuration of NZ-Mechanism will be sufficient to explain early necrosis adjacent to CV. To first establish credibility that Mouse Analog can be used to challenge the NAPQI zonation hypothesis, we needed to achieve a variety of prespecified Target Attributes (performance requirements). A Target Attribute is a characteristic feature of APAP induced liver injury to which a prespecified Similarity Criterion is assigned. Each wet-lab measurement that we seek to mimic, such as hepatic extraction ratio, becomes a Target Attribute. The Similarity Criterion specifies the degree of similarity that must be achieved. An example is the mean virtual experiment measurement falling within ± 1 standard deviation of the mean wet-lab experiment measurement. The following are three of the Target Attributes achieved: measurements of APAP 1) hepatic extraction ratio, 2) intrinsic clearance, and 3) dose-dependent pharmacokinetic profiles (Supporting S1 Fig). By increasing the strength and variety of analogies between measurements made during experiments on Mouse Analog and corresponding measurements made during experiments on mice (Fig 3A), the credibility of explanatory inferences drawn from virtual experiments increases. Hereafter, Mouse Analog components and characteristics are capitalized to distinguish them from mouse counterparts; feature and parameter names are italicized. For a particular Mechanism variant, Target Attributes are achieved by changing feature configuration values (Supporting S1 Table). Regardless of zonation, the probability of an event occurring within each Hepatocyte is fixed for the duration of the experiment. Probability of occurrence is the same within all Hepatocytes for events not subject to zonation. Note that to keep mechanisms parsimonious, some Fig 3B features, such as creation of Mitochondrial Damage (mitoD) objects, do not have direct mouse counterparts but instead subsume many. For NZ-Mechanism, the zonation of APAP Metabolism events (Fig 4A) causes the largest amount of NAPQI to be generated adjacent to CV. As depicted in Fig 3B, at each time step, each bound APAP may be Metabolized (maps to rate of metabolism); the probability of a Metabolic event is Lobule-location dependent, increasing from 0. 35 to 0. 95 PP to CV, which maps to increasing expression of participating metabolic enzymes. The probability that the Metabolite is NAPQI increases from 0. 33 to 0. 9 PP to CV (otherwise, G&S Metabolites are produced). That increasing probability maps to increasing expression of enzymes responsible for NAPQI formation, primarily CYP2E [1–4]. Glutathione (GSH) Depletion Threshold and probability of a mitoD Mitigation event are independent of Lobule location for NZ-Mechanism. Changing one or both alters total Hepatocyte Death but does not alter their Lobular locations. However, making one or both location-dependent (Fig 4B) changes the causal cascade and consequently alters Hepatocyte Death locations. Because hepatic lobules approximate regular polyhedra with CV located centrally, the percent of a given lobule’s hepatocytes within different regions varies dramatically. The same is true within Mouse Analog Lobules (Fig 4C). Note that Lobule does not have a direct mouse counterpart (Methods). Rather, it maps to a small random sample of all hepatic periportal to perivenous flow paths. Because of periportal interconnections among sinusoids, the number of hepatocytes encountered by a compound such as sucrose moving through a lobule is typically greater than the number of hepatocytes along straight-through sinusoids. Flow through a Lobule enables APAP and other mobile objects to mimic that phenomenon. The temporal patterns for average location of Death events and Death trigger events differ only by the value of Death Delay intervals, which range from 1. 2 to 12 hours (see Death Delay subsection in Methods). We explain in the Death Delay subsection that measuring Death events while keeping measurement of events separate from the Mouse Analog mechanism is computationally expensive. Thus, we elected to use measurements of Death trigger events for Fig 5 to challenge the NAPQI zonation hypothesis. Even though the highest probability of NAPQI formation occurs adjacent to CV, NZ-Mechanism is falsified because early average trigger events failed to fall within the Zone 3 target range close to CV (Fig 5). Falsification is discussed further in the Validation subsection under Methods. Prior to about 6 minutes, the location of early average trigger events shifts toward CV, but at about 10 minutes the trend shifts in the PP direction, away from CV. Such an unusual pattern suggests that NZ-Mechanism may not be biomimetic. Seeking evidence that challenges that prediction would be problematic in mice at such early times. Nevertheless, what is the explanation for this Mouse Analog phenomenon? Because of event stochasticity, some Hepatocytes at similar distances from the PP space will experience trigger events before their neighbors. Consider two regions: RCV (nearer to CV) and RPP (near PP space). The fraction of Hepatocytes experiencing an early trigger event will be larger in RCV. However, the total number of Hepatocytes in RPP is larger (Fig 4C). Consequently, the number of early trigger events in RCV and RPP can be similar, sometimes even larger in RPP. Thus, the average location of all early trigger events (circled trend 1, Fig 5) can appear skewed toward RPP. Later, remaining RPP Hepatocytes will be less vulnerable as the amount of APAP decreases, which causes the average trigger event location to shift toward CV. Concurrently, the fraction of Hepatocytes that have already experienced trigger events increased faster in RCV than in RPP, explaining the direction change around 10 minutes. The hypothesis that NZ-Mechanism alone is sufficient to explain early periportal necrosis is clearly falsified. We found no configurations that achieved that Target Phenomenon and shifted Hepatocyte Death patterns toward the CV (Fig 5). The results of these virtual experiments provide clear, strong evidence that an alternative, somewhat more complicated explanation is needed. Having falsified NZ-Mechanism, we posited that at least one additional feature of the Mechanism must exhibit zonation. The following two (Fig 4B) seemed equally plausible. 1) GNZ-Mechanism specifies that GSH Depletion Threshold values decrease PP to CV and NAPQI formation values increase PP to CV. 2) MNZ-Mechanism specifies that each Hepatocyte’s ability to mitigate mitoD diminishes sigmoidally PP to CV and NAPQI formation increases PP to CV. Predecessors of MNZ-Mechanism employing a parsimonious linear mitoD Mitigation gradient failed to shift trigger events sufficiently in the CV direction. Using a sigmoidal distribution improved similarity to the Target Phenomenon, but not sufficiently. The results require that both MNZ- and GNZ-Mechanisms be falsified. We inferred that combining the two Mechanisms into MGNZ-Mechanism might be sufficient to achieve the Target Phenomenon, and identified values of MGNZ-Mechanism configurations that did so. Having achieved the Target Phenomenon, MGNZ-Mechanism stands as a plausible explanation of the target APAP-induced necrosis pattern. Contents of each Hepatocyte were measured within three 5-grid-space-wide bands (Fig 6A) during execution of MGNZ-Mechanism. The APAP pharmacokinetic profile (Fig 6B) is sufficiently similar to a plasma profile in mice [16] (Supporting S1 Fig). To mimic hepatic blood flow, a particular Mouse Analog feature configuration determines the fraction of APAP in Mouse Body that is transferred to Liver each time step. G&S are transported out of Hepatocytes and allowed to accumulate in Mouse Body. NAPQI peaking first in CV region (Fig 6C) is a consequence of differences in APAP exposure, feature zonation, and the fact that Hepatocytes that have experienced a Death trigger event stop Metabolizing APAP. Mean amounts of G&S per Hepatocyte (Fig 6E) also reflect NAPQI patterns. G&S amounts in CV region prior to 20 minutes are greater than PP region amounts, which may seem counter-intuitive. Zonation of Metabolism is one explanatory factor, but zonal differences in Hepatocyte numbers (Fig 4C) are more important: 742 in Zone 3,3, 948 in Zone 2, and 9,310 in Zone 1. Prior to 80 minutes, Death has been triggered and APAP metabolism has been halted in most Hepatocytes in CV region but not Midzonal region. Thus, the value of cumulative Midzonal GSH Depletion events continues to increase after 80 minutes (Fig 6H). Having early mitoD values in CV region that are much larger (10x) than those in Midzonal region (Fig 6I) resulted in Hepatocyte Death occurring first near CV. The probability of a mitoD Mitigation event is smallest adjacent to CV, yet that is where the number of early Mitigation events is largest (Fig 6J). Therefore, although the probability of any one Mitigation Event is small, the cumulative number of such events is large. The space of MGNZ-Mechanism configurations contains both biomimetic and non-biomimetic variants. A variant that, for example, eliminates zonation of APAP metabolism (Fig 4A) is non- biomimetic. To test sensitivity and robustness of the MGNZ-Mechanism configuration, we conducted experiments using 64 plausibly biomimetic Mouse Analog variants of the MGNZ-Mechanism configuration in which values assigned to one-to-seven of ten influential configurations (shown in Supporting S2 Table) were changed. Each variant represents an arbitrary virtual mouse strain. Most changes included adjusting one or more zonation configuration. The Mechanism configurations included the five in Fig 4A and 4B plus the following: the probability that a NAPQI object will react, Death Delay value, the probability that a NAPQI reaction product will be mitoD or nonMD (non-mitochondrial damage products, discussed below), the probability of a nonMD Mitigation event, and/or Death Trigger Threshold value. Total Dead Hepatocytes 24 hours after dosing for all 64 variants were ordered in ascending order and plotted (Fig 7A). Comparing Deaths per zone rather than total Deaths (inserts, Fig 7A) provides concrete evidence of differences in how events within each Mechanism variant unfolded. Supporting S2 Fig provides additional examples. We compared those results to data from Harrill et al. [7] (red bars, Fig 7B) in which mean necrosis scores of 37 different mouse strains were recorded 24 hour after dosing with 300 mg/kg of APAP administered intragastrically. We assumed a direct correlation between necrosis score and both numbers of necrotic cells (in mice) as well as total Dead Hepatocytes (Fig 7B, right axis). All Mouse Analog data in Fig 7A were scaled the same so that at least one (of the 64) Mouse Analog bar heights was similar to the smallest and largest values in Fig 7B. We arbitrarily matched the 24-hour necrosis score for strain FVB/J to total Dead Hepatocytes to Mechanism variant #37, and then selected the one variant that provided the closest match to the necrosis score for each of the other mouse strains (white bar).
  Our intent is that parallel virtual and wet-lab experimentation will integrate into a faster and more effective scientific approach [17] (Fig 1). Most researchers seem to prefer including ample realistic detail in explanatory models. However, mechanisms like those employed in Mouse Analogs are scientifically preferable and more informative because they are only as detailed as necessary for explaining a particular phenomenon or challenging a particular hypothesis. A validated model Mechanism with limited detail subsumes many different yet equally plausible more detailed variants. When we falsify such a Mechanism (show that it is inadequate for its Target Attributes), we significantly shrink the space of equally plausible variants by eliminating whole classes of more detailed variants. The results in Fig 5 document the first use of virtual experiments to challenge competing biomimetic Mechanism hypotheses. The fact that three were falsified demonstrates how virtual experiments can shrink the space of plausible explanatory Mechanisms and facilitate the more fundamental work of wet-lab experimentation. When we started, we expected Mouse Analogs using NZ-Mechanism to show Hepatocyte Deaths occurring first near CV. We acquired new knowledge by examining why it, as well as MNZ- and GNZ-Mechanisms, failed to mimic that phenomenon. Explanatory insight was further improved by examining how MGNZ-Mechanism successfully achieved the Target Phenomenon. An example is the importance of zonated repair pathways in mitigating damage (Fig 6). We anticipate that other comparably parsimonious yet somewhat different Mouse Analog Mechanisms can do the same. The Iterative Refinement Protocol (described under Materials and Methods) outlines how to identify and challenge comparably parsimonious explanations of target phenomena. Mouse Analog concreteness, coupled with the ability to observe and measure propagating causal events, makes it easy for experts and non-experts to form and discuss opinions about the credibility of similarities between virtual experiment results and those from real or envisioned wet-lab counterparts. NZ-Mechanism was falsified because a small fraction of Hepatocytes in Zones 1 and 2 experienced an early Death trigger event (Fig 5) and consequently will experience an early Death event. The same phenomenon may occur in mice, but it has escaped detection because experiments were focused elsewhere. If deemed important, that phenomenon—that prediction—could be challenged directly using wet-lab experiments. During execution, MGNZ-Mechanism is a concrete hypothesis for how key features of APAP hepatotoxicity in mice are generated (Fig 3A). Suppose that strong analogies do exist between virtual and real mice. It is clear from concurrent differences within the three regions (Fig 6C–6H) that explanatory inferences drawn from whole liver measurements or hepatic biopsies will at best be flawed and likely misleading. Differences in zonal predictions can guide design of wet-lab experiments intended to challenge and possibly falsify the MGNZ-Mechanism hypothesis. For example, 1) at 30 min post-dose, GSH Depletion in the Midzonal region is significantly greater than in either CV or PP regions; and 2) at 10 min post-dose, counts of mitoD objects (maps to mitochondrial damage) in CV Hepatocytes are an order of magnitude greater than in PP Hepatocytes. Such experiments will be win-win: no matter the outcome, we will have useful new knowledge, and the new results are expected to motivate additional rounds of virtual experiments [17] (Fig 1). Note that, absent the data in Fig 6 there would be no basis to undertake such narrowly focused experiments. Within individuals, expressions of hepatic proteins are known to change significantly in response to a variety of environmental and health factors. We simulated possible examples (Fig 7A). Genetic differences among mouse strains that influence hepatic zonation more than they influence particular pathways or molecular level events may account for the diversity of differences in necrosis scores (Fig 7B). However, considerable uncertainty remains. As illustrated by the green bars in Fig 7A, for each Mouse strain explanation, there will exist a sizable set of Mechanism variants that would be equally explanatory. Further, MGNZ-Mechanism is just one realization drawn from a plausible yet finite Mechanism space. Automated methods, as they become available, can be used to more systematically explore that space, prior to or in parallel with new wet-lab experiments to challenge MGNZ-Mechanism. The data in Fig 7 demonstrate that, when targeting just one attribute, such as total Hepatocyte Death 24 hours post-dose, modest configuration changes for two or more Mouse Analog features can have counterbalancing consequences resulting in altered event details yet essentially no significant change in the measured feature. On the other hand, modest changes can have synergistic consequences resulting in a 3-fold change or more. Toxicity enhancing synergistic changes in cascade zonation, possibly initiated by changes in environmental and health factors, may be a contributor to idiosyncratic drug induced liver injury. In combination with the adaptation hypothesis [18], such change may help explain how and why only a small percentage of susceptible individuals develop overt liver injury.
  Methods stem from four requirements [19], which are needed to implement the use cases illustrated in Fig 3B. Meeting them is outside the scope of systems biology models currently in use to predict APAP hepatotoxicity [20,21]. We use the virtual experiment approach described by Kirschner et al. [17] along with the enhanced strategies detailed by Petersen et al. [22,23]. Feature names (Supporting S1 Table) are italicized. To achieve Requirement 2, Mouse Analogs are written in Java, utilizing the MASON multi-agent simulation toolkit [28]. In silico experiments are run using virtual machines [29] on Google Compute Engine, running 64-bit Debian 7. For longer simulations (e. g. , Fig 6) Monte Carlo trials are run in parallel. The data presented herein along with Mouse Analog code are available [30]. We started with a previously validated version of Liver [9]. The goal of an Iterative Refinement Protocol cycle is to refine a formulated explanatory Mechanism by achieving Target Attributes, thus completing one right-side cycle in Fig 1. A concrete software mechanism can be falsified—shown inadequate for its Target Attributes—in two ways: 1) it cannot exhibit a Target Attribute and/or 2) we fail to discover a Mouse Analog configuration that achieves all Target Attributes. Before starting, specify Target Attributes to be explained and rank-order them in terms of expected difficulty to achieve. For each Target Attribute, specify initial and possible subsequent Similarity Criteria. Liver composition (Fig 2) is now stable and robust, having already achieved several Target Attributes [9–12,31]. Bile Network does not influence the four Mechanisms. Because rat and mouse lobule structure and organization are similar, Liver can be used in simulations of phenomena measured in rats and mice. By altering analog-to-referent mapping functions [25], we can change experimental context. A Lobule maps to a small random sample of all lobular flow paths. The 3-zone network has 68 nodes, 45,20, and 3 in Zones 1,2, and 3, respectively. Nodes are connected using 99 edges: 55 from Zone 1-to-Zone 2; 10 within Zone 1; 24 from Zone 2-to-Zone 3; 10 within Zone 2; and 0 within Zone 3. To represent both uncertainty and variability, edges between Sinusoid Segments are randomly assigned at the start of each Monte Carlo execution. Having edges and Sinusoid Segment sizes assigned randomly at the start of each experiment simulates lobular variability within and between livers. See Validation subsection below for additional detail. Hepatocyte objects, excluding Fig 3B capabilities, are identical to those described in Petersen et al. [13,22]. Hepatocytes contain four physiomimetic modules [13]: InductionHandler, EliminationHandler, MetabolismHandler, and BindingHandler. The events and features illustrated in Fig 3B were added incrementally. Liver sinusoidal endothelial cells have been implicated in AILI. They contain CYP2E1 and GSH, although the relative amounts are very small compared to hepatocytes. An unanswered question is whether lobular differences in CYP2E1 and GSH within liver sinusoidal endothelial cells contribute to zonation of injury [32]. Addressing that question is feasible using Mouse Analog, given Target Attributes specific to liver sinusoidal endothelial cells, but it is outside the scope of this work. Adhering to a strong parsimony guideline, and given that the vast majority APAP metabolism occurs within hepatocytes, we specified that all relevant APAP metabolism and GSH Depletion occur within Hepatocytes. Within Mouse Analog Endothelial Cells, we specified that only nonspecific APAP binding occur. It is straightforward to add modules to Endothelial Cells and configure them appropriately should evidence become available that falsifies those specifications. The first performance goal was to discover Metabolism phase (Fig 3B) feature configurations that, upon achieving Target Attributes, would remain unchanged during experiments to challenge Mechanisms intended to explain necrosis patterns. There is quantitative variability (within and between experiments) in zonal measurements of CYP2E1 (primarily responsible for formation of NAPQI) [33], the relative proportion of the three primary metabolites (glucuronide, sulfate, and NAPQI) [34], hepatic clearance, and hepatic extraction ratio. However, there is qualitative agreement on their trends. An APAP object maps to a tiny fraction of the 300 mg/kg dose. There is a direct mapping between the probability of an unbound APAP object being metabolized each time step (1 second) and amounts of metabolic enzymes [13,22]. In mice, CYP2E1 levels per hepatocyte increase PP to CV by 2 to >10 fold. We specified this Target Attribute: probability of an APAP metabolic event generating NAPQI increases at least three-fold PP to CV. All other metabolites are lumped together and called G&S (maps mostly to the glucuronide and sulfate metabolites). In most reports, inactive metabolites are estimated to account for up to 85% of a dose, with NAPQI accounting for the balance. We specified that total NAPQI, as fraction of dose, be within 0. 15–0. 4; a large range was needed to ensure sufficient numbers of NAPQI are generated when dose is reduced. A Target Attribute is that toxicity be reduced by about 50% when the APAP dose is reduced by 50%. This Target Attribute was the most demanding: APAP hepatic extraction ratio = 0. 6 ± 0. 06. To simplify achieving that Target Attribute, we used a virtual single pass Liver perfusion protocol with a constant rate of APAP input. When APAP outflow stabilized within 34–46% of input values, the Target Attribute was reached. For convenience, we began the initial Iterative Refinement Protocol cycle by specifying that APAP Hepatic Clearance be similar to that of prazosin, for which the Liver had already been validated [10–12]. We cycled through the Iterative Refinement Protocol seeking changes to configurations that would enable achieving additional Target Attributes specified below under Mouse Body. The resulting configuration values are illustrated in Fig 4A and 4B. Lobules have a biomimetic PP to CV gradient [13] (Supporting S4 Fig), which Hepatocytes can use to create feature zonation. It maps to measures of one or more common blood attributes, such as pO2 [35]. However, because we needed the flexibility to explore multiple plausible, feature-specific PP to CV gradients for several different events, those in Fig 4A and 4B were implemented explicitly as functions of path length from analog Lobule PP entrance to the current Hepatocyte’s position. Path length (distance from PP entrance in grid spaces, dPP) is specified by this Hepatocyte’s distance from the inlet of its Sinusoid Segment (SS) plus the averages of the lengths of all Sinusoidal Segments in each of the previous zones. I. e. , dPP (SSz) = ∑i=z1⟨|SSji−1|⟩j, where SSz is a given Sinusoid Segment in Zone z and |SS| is its length and j indicates each SS in zone i; j indicates each SS in zone i; and dPP (Hss) = xss + dPP (SSz), where Hss is a Hepatocyte in a given Sinusoid Segment, and xss < |SS| is a position along the length of a Sinusoid Segment. The angle brackets indicate the average of its contents. Each gradient has a shape (linear or sigmoidal) and its own start and end values, which are then scaled over the common path length computation. Gradients are turned on or off independently by specifying their start and end values. Once these feature-specific gradients have stabilized, they can be replaced by rules that use the common biomimetic gradient. We started with a simple parsimonious operating hypothesis: when specific damage products accumulate within a Liver above some threshold, Death will be triggered irreversibly. The evidence shows that accumulation of APAP-induced damage products in mice occurs after available GSH has been sufficiently depleted. By analogy, we implemented a new feature: at each time step, with p = 0. 5 (specified arbitrarily because we had no wet-lab data to guide doing otherwise), each NAPQI object may be destroyed, which maps to in vivo depletion of a fraction of hepatocyte’s available GSH. A NAPQI destruction event is also a GSH depletion event. A small GSH Depletion Threshold value maps to mice that are very sensitive to APAP hepatotoxicity; a larger GSH Depletion Threshold maps to increased resistance. The Threshold value needed to be small enough to allow sufficient accumulation of Damage products (below), but large enough to achieve this Target Attribute: toxicity is reduced by about 50% when the APAP dose is reduced by 50%. Setting GSH Depletion Threshold = 5 for NZ-Mechanisms proved adequate. In wet-lab studies, GSH is often measured in whole liver homogenates. To make NZ- and MNZ-Mechanisms more directly comparable to GNZ- and MGNZ-Mechanisms at the level of whole Lobule measurements, we specified that Depletion Threshold = 3. 5 for NZ- and MNZ-Mechanisms, which is the average value for all Hepatocytes using GNZ- and MGNZ-Mechanisms. For some cases in Fig 7A, GSH Depletion Threshold values at the Lobule’s entrance of 8 or 3 are used (Supporting S2 Table). Each time step after Depletion Threshold is reached, with probability to react = 0. 5, each NAPQI object may be destroyed and replaced by a Damage product. Evidence implicates necrosis being triggered by mitochondrial damage. We parsimoniously specified that there be two classes of Damage product: “mitochondrial damage products, ” called mitoD, and “non-mitochondrial damage products, ” called nonMD. When a NAPQI object is destroyed, it is replaced with either a nonMD or mitoD object selected randomly (probability = 0. 5); we had no wet-lab data to guide specifying differently. To trigger Death, we implemented an analog counterpart to this simple yet widely accepted mechanism: upon accumulation of sufficient mitochondrial damage, a tipping point is reached and necrosis is triggered irreversibly. Triggering necrosis within one mouse hepatocyte requires reaction of possibly hundreds (or more) of NAPQI molecules, which is a very tiny fraction, f, of the administered APAP. One analog NAPQI object maps to a very small fraction, fA, of that same APAP dose. For mice resistant to APAP hepatotoxicity (case 1), it is possible that f > fA. When referent mice are much more sensitive to APAP hepatotoxicity (case 2), it is likely that fA > f. For case 2, replacing one NAPQI with one mitoD will be more than enough to trigger a Death event. Thus, the resolution [17] of one NAPQI → one mitoD is inadequate to simulate toxicity phase events. We solved that problem by specifying that each mitoD be amplified: one NAPQI → (1 + n) mitoD. Specifying that n be a random draw from either a uniform (1,6) distribution enabled achieving Target Attributes. The mitoD amplification step is an example of making an analog Mechanism locally finer grain (increasing resolution) without changing granularity elsewhere. As with Depletion Threshold, the Death Trigger Threshold value needed to be small enough to allow sufficient accumulation of mitoD without a trigger event, but also large enough so that toxicity is reduced by about 50% when APAP dose is reduced by 50%. Setting Death Trigger Threshold = 6 proved adequate. For some cases in Fig 7A, larger Death Trigger Threshold values are used. Hepatocytes utilize multiple mechanisms to mitigate or reverse different types of damage, ranging from up-regulation of GSH synthesis to mitochondrial autophagy [36]. Consistent with our strong parsimony guideline, we started with a single mitigation Mechanism, which maps to a conflation of all actual mitigation/recovery mechanisms. Each time step, with p = 0. 5 (specified arbitrarily), each nonMD and mitoD object may be destroyed. For NZ-Mechanism, lacking wet-lab data for guidance, we made probability-of-Mitigation event for each nonMD and mitoD at each time step independent of Lobule location. Because hepatocytes closer to CV are known to experience increasing oxidative stress, we inferred that mitigation mechanisms would be more robust closer to CV. However, when we increased probability-of-mitoD-Mitigation event (simply p-mitoD-Mitigation below) PP to CV, the locations of Dead Hepatocytes were shifted further in the PP direction, relative to the NZ-Mechanism data in Fig 5. To shift Death trigger events in the CV direction, it was necessary to decrease p-mitoD-Mitigation PP to CV. Instantiating a linear 0. 9-to-0 decrease in p-mitoD-Mitigation failed to produce a sufficient shift. We considered several options for moving mean trigger events toward CV. E. g. , keep p-mitoD-Mitigation constant at 0. 9 through Zone 2 and then decrease it linearly to 0 at CV. We rejected that feature because, absent any supporting observation, the risk seemed too great that it would be non-biomimetic. We considered adding a new biomimetic feature or event, but doing so conflicted with our strong parsimony guideline. Following several exploratory Iterative Refinement Protocol cycles, we opted for specifying that p-mitoD-Mitigation for MNZ- and MGNZ-Mechanisms decreases PP to CV following a reverse sigmoid with the inflection centered approximately midway in Zone 2. However, there is neither supportive nor refuting evidence for that configuration; therefore, this explanatory Mechanism hypothesis needs to be challenged. We could posit that the efficiency of one or more processes for mitigating normal mitochondrial dysfunction decreases PP to CV [36,37] and that one or more of the processes that maintain necrosis-triggering pathways is preferentially sensitive to NAPQI damage. Consequently, their combined effect maps to p-mitoD-Mitigation decreasing sigmoidally PP to CV in Mouse Analog. As with the GSH Depletion Threshold, to make NZ- and GNZ-Mechanisms directly comparable to MNZ- and MGNZ-Mechanisms at the whole Lobule level, we determined that the average p value for all Hepatocytes using is 0. 67. Therefore, we specified that p-mitoD-Mitigation be 0. 67 for NZ- and MNZ-Mechanisms in Fig 5 rather than 0. 5. Connecting Liver to Mouse Body produces Mouse Analog. Mouse Body contains a space that maps to all extrahepatic tissues including blood along with a space to contain dose, which enables simulating intravenous, intraperitoneal, and intragastric dosing. When additional details are required, new objects can be added without influencing preexisting Liver Mechanisms. To achieve APAP pharmacokinetic attributes, we adjusted values of two configuration features: one controls Absorption into Mouse Body and the other controls metering of APAP from Mouse Body to the Lobule’s entrance. The latter maps to hepatic blood flow. We specified that Absorption be first order and rapid. We prespecified two Target Attributes: 1) mimic single dose APAP pharmacokinetics in mice [16] and 2) observe characteristic dose-dependent pharmacokinetics for APAP in Mouse Body. For (1), we prespecified that Mouse Body measurements be within one standard deviation of mean wet-lab values. Supporting S1A Fig shows that the first Target Attribute was achieved. For (2), we found no fine grain dose-dependent pharmacokinetic data in mice to use as direct validation targets. However, there is pharmacokinetic data in rats [38,39] demonstrating characteristic dose-dependent pharmacokinetics. We reasoned that, by using a simple mouse-to-rat scaling, reported dose-dependent pharmacokinetic data in rats [38,39] could serve as alternative targets, despite the fact that the metabolism and toxicity of APAP in mice and rats are quite different (we are targeting dose-dependent APAP clearance, not metabolite types). For comparable mg/kg doses, the plasma half-life of APAP in rats is approximately 3x-5x that in mice [40]. The dose of APAP objects used in Fig 6B maps to 300 mg/kg in mice. Using the median value of 4x, we determined that dose should scale to approximately 75 mg/kg in rats. Twice that dose of APAP objects should map to approximately 150 mg/kg, and half that dose should map to approximately 37. 5 mg/kg in rats. Supporting S1B and S1C Fig show that characteristic dose-dependent pharmacokinetics is observed. Attaining Target Attributes from different species increases credibility and documents achieving the second capability under Requirement 2. There are qualitative but not quantitative time-course data for relative amounts of necrosis in mice. Necrosis is not detectable during the first hour following a toxic but non-lethal APAP dose of 300 mg/kg, and there is no evidence of significant further induction of necrosis 12 hours after dosing [3]. Covalent adduct formation and histological evidence of hepatocyte damage always precedes measurable necrosis by tens of minutes to a few hours. That time interval maps to Death Delay. To specify the Death Delay feature, we needed a Target Attribute estimate for percent total necrosis as a function of time post-dose. Target range estimates (see Supporting S5 Fig) were provided by coauthor Kaplowitz for percent total necrosis at several times post-dose. The prespecified Similarity Criterion was that cumulative mean Death events fall within the estimated Target Attribute range. That Target Attribute was achieved in two different ways (Supporting S5 Fig): draw each Hepatocyte specific Death Delay value from a uniform (a pseudo-random draw from uniform [1. 2,12] hours) or a normal distribution. Although the second seemed more realistic, we used the first because it is more parsimonious and doing so enabled achieving the Target Attribute range. Necrosis is the target wet-lab phenomenon. Our working hypothesis is that necrosis is a process that is triggered by early damage and that the wet-lab measurement is detecting a late stage in that process. Death Delay maps to the process of necrosis. The time of an Hepatocyte Death event = time of that Hepatocyte’s Death trigger event + a pseudo-random draw from [1. 2,12]. In Fig 5, we elected to focus on Death trigger events rather than time of an Hepatocyte Death events for three reasons. 1) Although there is no current wet-lab method to measure its short transient duration, the trigger event is a key causal event. 2) We know that final Lobular locations of trigger events and Dead Hepatocytes are the same. 3) To make the virtual experiment as much like a wet-lab counterpart as possible, we strive to measure the analog in the same way one measures the referent. Conventional models are not measured; rather, they simply output variables (e. g. concentrations of species) each time step. In contrast, Mouse Analog is measured. Doing so during executions requires keeping Death event measurements separate from Analog Mouse Mechanisms. Measuring Death events requires a high frequency polling effort, where during one poll the observer agent records a Hepatocyte in a particular location as not Dead, and the next poll records that Hepatocyte as Dead. Such polling increases considerably the time required to complete each simulation experiment. Mouse Analogs are treated as a form of data, using both the implicit schema of Java, JavaScript, and R and the explicit schema of the configurations. Mouse Analogs and configuration data are maintained, archived, and released using the Subversion version control tool in two repositories, one private (Assembla) for rapid and prototyping development with project partners and another public for collaboration [40,41]. Input-output (I/O) data is handled separately. Smaller data sets are stored in simple CSV. I/O data is tightly coupled to experiments and configuration details, requiring a common versioning system, aggregating all three data types. Versioned I/O data is archived as downloadable packages. The entire toolchain, including the operating system, used for Mouse Analogs, configurations, and I/O handling is open-source, thereby ensuring repeatability. Similarly, all project generated and released data is available to be licensed as open data. Mouse Analogs are built for, maintained, and executed in a cloud environment (e. g. , Google Compute Engine) to ensure platform and infrastructure repeatability across experiments, project team members, partners, and the wider community. Regression and unit tests are special cases of canonical use cases (e. g. , a single pass Liver perfusion experiment). While use cases are instances of the class of experiments to which Mouse Analogs are being applied (e. g. , AILI), they also provide the measures by which the software and methods are maintained. Each toolchain iteration can execute the canonical use cases so that current results can be compared to prior results for which some degree of credibility has been documented. Significant variations are documented, investigated, and explained. While the majority of variations are the result of model iterations, the process does catch artifacts introduced by changes. Unexpected variations in results during an Iterative Refinement Protocol cycle first trigger a software verification process designed to test the toolchain from the most general and widely used, up to the most specific tool: machine, OS, compiler, libraries, simulator, and model. Absent unexpected variations, that same process is applied at least yearly. So doing is rigorous and avoids wasteful iterations that result from assuming that all variations are caused by the most specific layer. Should anomalies occur at the Mouse Analog layer, they would trigger a source, data, configuration setting, and trace review in comparison with the design of experiments. Any unexpected variations that survive the verification process are then submitted to an Iterative Refinement Protocol based falsification battery designed to invalidate the Mouse Analog as compared to wet-lab observations. Results satisfying canonical use cases, but yielding unexpected variation as a result of changes to Mouse Analog or toolchain, provide material that can be used to formulate a useful hypothesis.
summary: Acetaminophen-induced liver injury in mice is a model for drug-induced liver injury in humans. Challenging an explanatory mechanism in mice is problematic because variables determining causes and effects cannot be controlled adequately. We circumvent that impediment by performing virtual experiments that challenge the prevailing scientific explanation for the characteristic spatiotemporal pattern of early acetaminophen-induced hepatic necrosis. Virtual mice utilize a biomimetic software liver. Results of virtual experiments provide compelling evidence that the prevailing explanation is insufficient. Without further studies in mice, we discovered a new, marginally more complex explanatory mechanism that met stringent tests of sufficiency. We argue that this virtual causal mechanism and the actual mechanism in mice are strongly analogous within and across multiple biological levels. Variants of the virtual mechanism provide possible explanations for the considerable variation in 24-hour necrosis scores among 37 genetically diverse mouse strains following a single toxic acetaminophen dose.
keywords: biotechnology
  cell death
  medicine and health sciences
  liver
  pathology and laboratory medicine
  engineering and technology
  cell processes
  enzymology
  vertebrates
  mice
  biomimetics
  animals
  mammals
  signs and symptoms
  enzyme metabolism
  pharmacology
  enzyme chemistry
  research and analysis methods
  bioengineering
  necrotic cell death
  animal cells
  mathematical functions
  mathematical and statistical techniques
  hepatocytes
  pharmacokinetics
  necrosis
  biochemistry
  rodents
  diagnostic medicine
  cell biology
  anatomy
  biology and life sciences
  cellular types
  amniotes
  sine waves
  organisms
title: Competing Mechanistic Hypotheses of Acetaminophen-Induced Hepatotoxicity Challenged by Virtual Experiments