Hepatic drug disposition is different in normal and diseased livers. Different disease types alter disposition differently. What are the responsible micromechanistic changes and how do they influence drug movement within the liver? We provide plausible, concrete answers for two compounds, diltiazem and sucrose, in normal livers and two different types of cirrhotic rat livers: chronic pretreatment of rats with carbon tetrachloride (CCl4) and alcohol caused different types of cirrhosis. We started with simulated disposition data from normal, multilevel, physiologically based, object-oriented, discrete event in silico livers (normal ISLs) that validated against diltiazem and sucrose disposition data from normal livers. We searched the parameter space of the mechanism and found three parameter vectors that enabled matching the three wet-lab data sets. They specified micromechanistic transformations that enabled converting the normal ISL into two different types of diseased ISLs. Disease caused lobular changes at three of six levels. The latter provided in silico disposition data that achieved a prespecified degree of validation against wet-lab data. The in silico transformations from normal to diseased ISLs stand as concrete theories for disease progression from the disposition perspective. We also developed and implemented methods to trace objects representing diltiazem and sucrose during disposition experiments. This allowed valuable insight into plausible disposition details in normal and diseased livers. We posit that changes in ISL micromechanistic details may have disease-causing counterparts.
Liver cirrhosis alters hepatic drug disposition, complicating drug therapy management (Le Couteur et al., 2005; Dourakis, 2008). The nature of alterations is dependent on both the cause and the extent of disease. Improved mechanistic insight is needed on two fronts. We need supported concrete theories for 1) how the interaction of a drug with the hepatic microarchitecture contributes to overall measures of disposition and 2) how disease progressively alters those microarchitectural features. Achieving both is complicated by the heterogeneity of hepatic microarchitectural features (for examples and a discussion, see Liu and Pang, 2006) and differences in cirrhosis. We focus on two standard rat models of cirrhosis: chronic treatment with 1) carbon tetrachloride (CCl4) and 2) ethanol (Hung et al., 2002a,b). In advanced stages, CCl4 treatment produces acute hepatocellular injury with centrilobular necrosis and stenosis. In contrast, chronic alcohol treatment produces hepatocellular injury with inflammation and perivenular macrovesicular steatosis (fatty degeneration). Both treatments cause fibrosis.
We report significant progress in achieving both objectives by developing, refining, validating, and experimenting on discrete, object- and agent-based in silico livers (ISLs). ISLs are works that are in progress. Observable micromechanisms in ISLs map directly to wet-lab counterparts, which facilitates the falsification of hypotheses approximately those mechanisms. See Hunt et al. (2009) for complete descriptions of the modeling method. We matched disposition profiles from ISLs and from in situ perfusion experiments in normal rat livers for diltiazem and sucrose. Matches were achieved through cycles of iterative ISL refinement that narrowed the space of plausible, spatially distributed micromechanisms. Many future ISL iterations are likely. Nevertheless, we discovered a limited set of ISL changes that produced profiles that acceptably matched in situ experiments on the two types of diseased livers (Hung et al., 2002a,b). Diseased ISLs were created independently.
Parameterizing detailed multizone versions of familiar physiologically based pharmacokinetic models has proven challenging because of hepatic heterogeneity coupled with interindividual variability (Liu and Pang, 2006). We pursued a fundamentally different approach. A dictum of the physicist Richard Feynman was “what I cannot create, I do not understand.” Uncertainty is large and detailed data are limited for micromechanisms that are responsible for differences in hepatic disposition profiles. Hence, it follows that to gain insight, we need to build extant (actually existing, observable) working mechanisms that exhibit some of those same phenomena. We must construct, validate, and explore analog mechanisms and tissues to better understand the biological counterparts. We built extant biomimetic mechanisms and tissues using object-oriented software tools. In doing so, we were not building a model using published biological facts, ideas, and assumptions. Rather, we constructed reasonably realistic biomimetic mechanistic hypotheses. We then explored and shrank the space of the resulting disposition profiles to achieve a prespecified similarity to wet-lab profiles.
Our focus was on constructing and falsifying plausible biomimetic mechanisms at multiple levels. We argue that the causative, mechanistic details executed in ISLs may have wet-lab hepatic counterparts during diltiazem disposition, as diagrammed in pair A of Fig. 1. Differences in mechanistic details between normal and the two “diseased” ISLs (pair B in Fig. 1) are hypotheses approximately corresponding to differences between the normal and diseased livers. Simulation enables testing those hypotheses. Achieving a degree of profile matching is evidence supporting those hypotheses. The differences in dynamic, multilevel details during execution provide plausible explanations of differences between the two disease models. Simulations of transformations of normal into diseased ISLs stand as abstract theories for disease progression: similar transformations may have occurred in rats during disease progression.
An important advantage of this class of models is that micromechanistic details are observable, unlike real livers. Methods were implemented to trace, measure, and record changes in dynamic spatiotemporal details. Doing so provided a view, not attainable until now, of how and where micromechanistic events combine to influence disposition within ISLs. We posit that changes in micromechanistic details from normal to diseased ISLs may have disease-causing hepatic counterparts. These methods are extensible to whole organisms and, eventually, patients. Hence, they open a door to new experimental means of testing the plausibility of mechanistic explanations. That is expected to facilitate translation of research results to benefit patients.
Materials and Methods
Rationale for Modeling and Simulation Approach.
Classical physiologically based pharmacokinetic (PK) modeling (Hung et al., 2001, 2002a,b) and the approach used herein, although fundamentally different, are complementary approaches to gaining insight into plausible mechanisms responsible for PK data. Classical modeling provides a conceptual generalization; a global description of flow, influx, efflux, binding, sequestration, and metabolism; and relates the resulting parameter values conceptually to observed changes in liver pathophysiology and biochemistry. Neither the model nor its parameters can be made realistic and similar to hepatic features (Rescigno, 2004). The synthetic method of modeling and simulation is fundamentally different (Hunt et al., 2009). The method enables ISLs to use actual mechanisms—events occurring within specified spaces—in a manner that is more consistent with the actual intralobular arrangements, normal and pathological, than is possible using traditional mathematical models. An ISL instantiates a hypothesis (Fisher and Henzinger, 2007; Hunt et al., 2008). Execution and comparison of results to referent counterparts tests the hypothesis. The relationship among ISL spaces, components, mechanisms, and phenomena—mappings 1 to 3 in pair A of Fig. 1—can be made increasingly realistic and similar to hepatic counterparts. Doing so is facilitated by grounding the internals of the ISL relative to each other rather than to absolute metric units such as seconds or meters. Relative grounding makes the causal effects of changes to any one component evident as they ripple through the entire model and are measured by an observer. This capability dramatically increases the extent to which a model can be refined and translated to other contexts. For a discussion of grounding, see Hunt et al. (2009). For convenience, a brief discussion of how grounding influences differences between traditional, inductive, equation-based PK models (left side of Fig. 1) and synthetic, internally grounded analogs such as ISLs is provided in the discussion in supplemental data.
Most ISL form and function details used herein are the same as used previously (Park et al., 2009). A dense, abridged description follows. Some of the tracing features, detailed in supplemental data, are new. Included is a discussion and example of technical issues related to the granularity of the tracing and the granularity of the computation. For convenience, parameter descriptions are also provided in the Appendix. Detailed descriptions of ISL design considerations, including ISL-to-liver mappings are available (Hunt et al., 2006; Yan et al., 2008a,b,c; Park et al., 2009). ISL mechanistic details are presented as plausible approximations of what actually occurred. To avoid conflating ISL details with the biology, and to clearly distinguish in silico components and processes from hepatic counterparts, hereafter we use small caps font when referring to the ISL counterparts.
In Situ Liver Perfusion Studies.
Full details of the original single-pass, liver perfusion experiments along with an explanation for the choice of compounds studied are provided by Hung et al. (2001, 2002a,b). Included are descriptions of two established methods to induce fibrotic, hepatic cirrhosis in rats (150-g male Wistar rats). Normal and two types of diseased livers were studied. Similar pretreatment protocols were used to create the eight diseased livers for each group: one group was produced by chronic CCl4 treatment and the other group by chronic alcohol (ethanol) treatment. Both protocols induced acute hepatocellular injury, but their histologies were different. Each model type reflects different aspects of human disease. Control, normal PK profiles were obtained using livers from matched rats treated identically, absent either CCl4 or alcohol treatment. Several histopathology measures characterized the nature and extent of disease. Nine outflow profiles of coadministered diltiazem and sucrose were analyzed individually using established PK methods. The referent data are presented in fig. 1 of Hung et al. (2002b).
An ISL is a simulation framework: an in silico counterpart to an entire wet-lab experimental system (analytical instrumentation and all). It comprises an experiment agent, a data management module, a statistical observer module (used to analyze data), a parameter manager, data from the perfusion experiments (and interpolations, when needed), RefModel, and lobule. RefModel is the parameterized classical PK model, which was fit by Hung et al. (2002a) to the wet-lab data. Concurrent execution of it and an ISL enabled judging similarities or lack thereof during iterative ISL refinement (discussed below). RefModel is a two-phase stochastic liver model (stochastic PK model hereafter): it predicts the time course of diltiazem in liver effluent during a perfusion experiment that follows a standardized dosing protocol. Its implementation is provided in supplemental data.
A lobule is one Monte Carlo variant of the complete system illustrated in Fig. 2. Three lobule variants (normal and two diseased) were developed in parallel. For simplicity, we began by assuming that anatomical, physiological, and PK characteristics of all hepatic lobules within a specific liver, normal and diseased, are somewhat similar (Hunt et al., 2006; Yan et al., 2008b). Pooled and averaged results from ISL executions of 48 Monte Carlo variants of a single, parameterized lobule represents a single wet-lab outflow profile; those results comprise one ISL experiment.
Figure 2 shows that lobules are abstractions: they are not intended to be accurate, precise descriptions of hepatic physiology. ISL components can be modified and plugged together in different ways as needed to represent different lobular properties or the consequences of disease. ISL experiments help reduce uncertainties about hepatic mechanisms by enabling the formulation and testing of fine-grained mechanistic hypotheses about plausible, fine-grained mechanistic details that may be occurring during drug disposition (Hunt et al., 2009). Software objects represent spatiotemporal aspects of hepatic organization and function. The consequences after execution are measured and studied simultaneously, analogous to how wet-lab experiments are conducted. An ISL—normal or diseased—achieves a degree of validation when the similarity between its outflow profile and a referent perfusion profile is judged adequate based on some quantitative comparison. Once that has been achieved, we state that the mappings marked 2 in pair A of Fig. 1 are plausible. As discussed recently (Hunt et al., 2009), a lobule should be no more complicated than is needed to achieve the stated objective. For this study, the objective has been to discover a minimal set of normal ISL parameter changes that will lead to diseased ISLs that achieve a degree of validation against referent outflow profiles.
When ISL components and their arrangement are judged acceptable, the detailed micromechanisms causing traceable events may correspond to the hepatic micromechanisms; i.e., mappings 1 in pair A of Fig. 1 are plausible. At that stage, the traced diltiazem and sucrose dynamics within and between the six levels of an ISL provide previously unavailable insight into plausible drug disposition details.
A lobule Is a Network of Sinusoidal Segments.
The relative arrangement of hepatic function and blood flow is represented at the lobule level using a directed graph called a sinusoid network. Each Monte Carlo variant maps to a distinct arrangement of flow paths from portal vein tracts (PV) to the central vein (CV) within a portion of an acinus. A sinusoid network is subdivided into three zones. Zonation enables mimicking quantitative and functional differences between periportal and perivenous lobular regions (Miller et al., 1979; Gumucio and Miller, 1982; Gaudio, et al., 1997). It is a means of achieving an adequate variety of PV-to-CV flow paths. sinusoid network structure properties are specified as a topology of nodes and edges. The number of nodes per zone is always zone 1 > zone 2 > zone 3. A graph edge specifies a flow connection from a sinusoidal segment (SS) exit to a downstream SS entrance. Having edges assigned pseudorandomly at the start of each ISL experiment simulates lobular variability within and between livers.
A SS is a software agent (an autonomous object that schedules its own events and interacts with other agents and objects in its environment) that represents all aspects of sinusoid function that can influence drug disposition. One SS is assigned to each graph node. Each is somewhat different, and the stochastic differences are parameter controlled. SSs per zone were zone 1 = 45, zone 2 = 20, and zone 3 = 5. As explained previously (Yan et al., 2008b), those numbers were needed to have sufficient PV-to-CV path variety to reduce fluctuations within outflow profiles. They were connected using a minimum of 109 edges: 39 intrazone 1 edges, 8 intrazone 2 connections (but no intrazone 3 connections), 37 zone 1-to-zone 2, and 25 zone 2-to-zone 3 connections. All zone 3 nodes were connected to CV. There were two constraints: no self-self edges and no two-node cycles (a restriction that was not imposed previously; Hunt et al., 2006; Yan et al., 2008b,c), and if any node was, by chance, not assigned an outgoing edge, it was connected directly to CV. Two SS types were used: direct (larger, shorter; controlled by the DirSin parameters in Table 1) and tortuous (thinner, longer; controlled by the TortSin parameters in Table 1).
A SS consists of a core and three identically sized layered toroidal spaces as illustrated in Fig. 2. Spaces A to C within a SS are identical, but SS sizes (Table 1) are Monte Carlo-specified. The spaces are subdivided into a parameter-specified number of square grid spaces. The core maps to blood flow. It provides a direct PV-to-CV path through which compounds (mobile objects) can traverse. Space A maps to the interface between vascular flow and the endothelial layer. Space B is called the endothelial layer. It maps to easily accessible spaces and cells presumed to be primarily endothelial cells and fenestrae. Space C is called the hepatocyte layer. It maps to less accessible spaces and cells, primarily the space of Disse, hepatocytes, and bile canaliculi. Cells in space B are called endothelial cells; those in space C are called hepatocytes. Parameters allow the resolution of the spaces to be changed as needed. A bile space, not illustrated in Fig. 2, can be added when needed.
Cells contain whatever objects are needed to represent required intracellular processes, such as drug binding, metabolism, transport, and sequestration. Because the in situ perfusions had short durations (<1 h), we assume that cell biology and biochemistry were relatively constant. Consequently, details not needed are abstracted away but can be added easily when needed (Hunt et al., 2006; Yan et al., 2008b; Park et al., 2009). It is known that basic compounds such as diltiazem are sequestered in organelles such as lysosomes and mitochondria, as well as being bound to intracellular material (Hung et al., 2002a; Siebert et al., 2004). However, motivated by parsimony, sequestration and binding are not resolved: everything within a cell that can bind or sequester diltiazem is conflated and represented by some number of identical binding objects (hereafter, simply binders). Within hepatocytes, we do not resolve binding to metabolic enzymes, such as the cytochrome P450 isozymes, and binding to or sequestration by other cell components. binders called enzymes handled binding inside hepatocytes. They use a parameter (metabolizeProb) to determine which of their binding events ends with release of metabolite. When needed, several different objects that produce the same net event sequence can replace a binder enzyme.
As in Park et al. (2009), ISLs were iteratively refined (Hunt et al., 2009) to achieve a set of targeted attributes along with a rigorous degree of prespecified similarity with referent outflow profiles. An illustration of iterative parameter adjustment is provided in supplemental data. No new attributes, other than the additional disease-altered outflow profiles were added. ISL parameters specify structural and functional properties, experiment configuration, and dose. The ISL and its components from earlier work were reused (Yan et al., 2008b; Park et al., 2009).
Monte Carlo Runs.
An ISL experiment averages 48 Monte Carlo trials (lobule). Execution duration in simulation cycles is specified by cycleLimit. The number of steps executed each cycle is stepsPerCycle. Drug disposition is observed at cycle resolution; spatiotemporal activities are traced at step resolution. One simulation cycle maps to 0.5 s, and a step maps loosely to 0.25 s. CycleLimit and stepsPerCycle were set to 200 and 2, respectively. ExperAgent is absolutely grounded, whereas lobules are not. A technical issue dealing with how time is discretized at the simulation cycle and step levels is discussed in supplemental data.
Parameters named in this paragraph govern structural and spatiotemporal properties. GraphSpecFile is used to create a lobule's sinusoidal network at the start of each simulation. DirSinRatio and TortSinRatio specify the ratio of the two SS types. SS circumference and length are generated using the values of DirSinCircMin/Max, TortSinCircMin/Max, DirSinLenAlpha/Beta/Shift, and TortSinLenAlpha/Beta/Shift (Hunt et al., 2006; Yan et al., 2008b,c). A2B/B2A/B2C/C2BJumpProb governs probabilistic compound movement. Simulated blood flow (in the core) and local, biased random walk (spaces A–C) are controlled by CoreFlowRate and SinusoidTurbo; their values can depend on the compound's physicochemical properties. ECDensity and HepDensity specify the fractions of spaces B and C occupied by endothelial cells and hepatocytes, respectively. The numbers of binders and enzymes within each cell is Monte Carlo drawn within the range specified by BindersPerCellMin/Max. SoluteBindingProb designates the probability that a drug will be bound. Each binding event lasts SoluteBindingCycles. The latter two parameters also depend on the compound's physicochemical properties (Hung et al., 2001; Mager, 2006).
MetabolizeProb is the probability that a metabolite, rather than diltiazem, will be released from an enzyme-diltiazem complex. ISLWetLabScaling is a scaling factor used to map compounds exiting CV directly to drug concentration (in perfusate). MembraneCrossing (valued 0,1) specifies whether a particular compound is allowed to enter cells.
Differences between normal and diseased ISLs.
The histopathology data showed that within both CCl4-treated and alcohol-treated livers, significant cirrhotic change occurred at cellular and subcellular levels (Hung et al., 2002a,b). Microscopy evidence showed that microvascular and microcirculation changes occurred (Gaudio et al., 1997). For simplicity, many of the details known about liver fibrosis, including the principal role played by hepatic stellate cells, for example, were not added to the list of targeted attributes. It is tempting to assume that those visible changes must have contributed in important ways to the observed alterations in outflow profiles. Yet, we do not know the dispositional significance of such changes for specific compounds. We elected to make no inferences but seek one of the simpler sets of changes that could provide a plausible explanation. It may be possible, for example, that two different sinusoid networks, with all cellular and subcellular details being unchanged, can provide an explanation for the differences. However, Hung et al. (2002a,b) showed significant correlations between measures of histopathology and changed PK parameter values. We therefore elected to begin by reusing the validated normal lobule network structure the diseased lobule network structure. We then looked for explanatory change at the SS level and below. Note that when required, a validated simpler micromechanism can be made more complicated.
Starting with the validated normal ISL (Park et al., 2009), the protocol of each refinement cycle followed three steps. That ISL provided a similarity measure (SM) > 0.8, discussed below. 1) We identified a subset of parameters to be changed (everything except the zonation, primary flow paths and their arrangement) and then sought focused changes in SS component parameterizations that would alter the outflow profile in ways consistent with the diseased outflow profile. More detail is provided under ISL Implementation and Execution. 2) We tested the hypothesis that a new valuation of the subset, in combination with unchanged values of all other parameters, would yield an ISL that could achieve SM > 0.8. When falsified, we returned to step 1. 3) We fine-tuned the diseased ISL parameter vectors with the objective of achieving SM > 0.9 (within a factor of 0.33 of the referent value, as discussed below). When that failed, we returned to step 1 or 2.
ISL parameter influences are not independent (Hunt et al., 2006). More than one parameter vector can give essentially indistinguishable outflow profiles (Yan et al., 2008b). Although every effort was made to construct a minimal model, the requirement that lobular structures be derived or inferred from published hepatic knowledge rather than induced from particular data sets provides some complexity and overlapping phenomena.
Drug Input, Dosage Time Management, and Similarity Measure.
ISL experiments followed the same dosing protocol used in situ (Yan et al., 2008b,c; Park et al., 2009). As illustrated in Fig. 2, a bolus dose of sucrose and/or diltiazem was injected into a simulated catheter that feeds into PV. Compounds were collected as they entered CV, simulating collection by a fraction collector. Details of dosing within the ISLs are provided in Park et al. (2009) and for convenience in supplemental data.
An ISL outflow profile was accepted as valid—as being indistinguishable experimentally from a profile obtained from a repeat wet-lab experiment—when SM > 0.8. Once that was achieved, it was increased to SM > 0.9. ISL outflow profiles were compared with referent profiles using the quantitative SM used previously (Hunt et al., 2006; Yan et al., 2008b,c; Park et al., 2009). For convenience, details are provided in supplemental data. An ISL profile was accepted as being experimentally indistinguishable when the fraction specified by the SM was within the ± 33% band around the referent values (Fig. 3).
ISL Implementation and Execution.
ISLs were implemented within the high performance, computing infrastructure diagrammed in Supplemental Fig. S1. Its methods were designed to provide improved life cycle management, execution efficiency, tracing quality, and analysis of traced results (Ropella et al., 2003). Each parallel mode was associated with one of the six ISL levels illustrated in Fig. 2 or an experimental requirement. Heterogeneity in parallel execution helped achieve improved performance along with efficient resource management. ISL parallel mode (Supplemental Fig. S1) was supported at group and experiment levels. Group level parallel mode enabled executing multiple experiments in parallel by segregating each and allowing each to run concurrently without interaction. Parallel batch processing and analysis of local execution results were performed using that mode. Experiment level parallel mode enabled executing single experiments in parallel as separate lobule Monte Carlo variants.
Parallel parameter sweeping (Supplemental Fig. S1) was used to construct and explore regions of ISL parameter space for vectors that enabled the ISL outflow profile to meet or exceed the specified SM. The location of such a region was deduced by coupling prior ISL parameterization experience with physiologically based heuristics. The latter helped identify regions that contain abiotic parameterizations: the resulting ISL cannot map to a liver. Occasionally, a region was selected randomly. A global search of the entire parameter space was impractical.
Within the experiment framework, a simulation coordinator managed the overall simulation life cycle by controlling the top-level agent ExperAgent. It also supervised system components including the parallel parameter sweeper, parallel batch processor, parallel model partitioner, and model deployer. Their integration provided a fully automated, high-performance simulation environment within which all experiments were conducted. We built the environment using the tools specified in supplemental data.
Multiscale Event Tracing within ISLs during Simulations.
Tracing disposition events within and between ISL levels was divided into two phases: generating tracing data and then evaluating the data using quantitative measures. First, all spatiotemporal events involving compounds (of the same type) across all ISL levels were stored. Tracing measures were then derived to include temporal changes within a SS, the lobular components each compound encountered plus the path each compound traversed, the length of each traversed path, and each compound's resident time within an ISL or one of its features.
During the first phase, two types of raw tracing data were generated. A trace data file was generated for each SS including PV and CV. It recorded the temporal order of events experienced by all compounds that resided within a particular SS. To trace metabolic events, each SS node also generated a tracing file that listed the compound's identification and type along with the time metabolism occurred.
Collected data were evaluated during phase 2. The tracing process and steps within are diagrammed in Supplemental Fig. S2. We started by tracing changes at levels 1 and 2 (Fig. 2): changes associated with endothelial cells, hepatocytes, enzymes, and binders at level 1 and within the four SS spaces. Tracing results at higher levels were deduced by aggregating results from lower levels. A traverse path provided an image of each SS node visited by a compound during its course through a lobule. Not all injected compounds reached the CV. Some were metabolized. Others were retained because of multiple binding events. All paths started at PV; they ended at either CV or with the SS where the compound was metabolized or remained when the run terminated. Traverse length summed the length of each SS entered. The fine-grain paths taken by a compound within levels 5 and 6 were ignored. Regardless of path taken by a compound within a SS, its traverse length was simply the length (in grid spaces) of that SS. Residence time was the time spent within the ISL before exiting, being metabolized, or having the run terminate. Each component that encountered a compound recorded the number and type of compound with which it interacted.
Initial Iterative ISL Refinement.
We sought outflow profiles produced by normal and two types of diseased ISLs, derived from that normal ISL, that would achieve SM > 0.9. We named the latter two diseasedCCl4 and diseasedAlc. Absent evidence that major features of lobular anatomy (zonation, primary flow paths, and their arrangement) were meaningfully different between the three liver types, we specified that in the ISLs they would be essentially the same. Previously (Park et al., 2009), we described a diseasedCCl4 ISL that validated against a diltiazem outflow profile from a diseased, CCl4-treated liver. It was achieved by adjusting ten parameters of a normal ISL that had validated against a diltiazem outflow profile from a normal liver. Adhering to the parsimony guideline, we sought an alternative parameterization (limited to the same ten parameters) of that normal ISL that would also produce a diseasedAlc ISL that could achieve SM > 0.9. We sampled ISL parameter space many times and found alternative normal ISLs that could be transformed to either a diseasedCCl4 or diseasedAlc ISL (for which a SM > 0.9 could be achieved), but not both. That failure falsified that earlier, normal ISL as a starting model for achieving both diseasedCCl4 and diseasedAlc ISLs. An alternative, validated normal ISL was needed. Note, however, that earlier, normal ISL remained valid when outflow profiles from only normal and CCl4-treated livers were targeted. Again, adhering to the parsimony guideline, we increased the number of parameters available to adjust from 10 to 11 and subsequently to 12, which proved adequate. The two additional parameters that were available for adjustment were endothelial cell density (ECDensity) and SinusoidTurbo, which biases compound random walk in the direction of the CV.
Validation of Disposition in normal and diseased ISLs.
As described under Materials and Methods, we specified that achieving SM > 0.9 within a factor of 0.33 of referent data were an acceptable target: an outflow profile that meets that SM was accepted as being experimentally indistinguishable from its referent profile. Through use of the new parameter sweeping capability (Supplemental Fig. S1), we iteratively adjusted the new normal ISL parameter values (Table 1) to move outflow profile properties toward those of diltiazem and sucrose from CCl4-treated livers. We continued that process until the outflow profile achieved SM > 0.9. We then repeated that protocol with outflow profiles from an alcohol-treated liver to obtain validated diseasedAlc ISLs. Parameter adjustments that enabled achieving SM > 0.9 are diagrammed in Fig. 4. Diltiazem outflow profiles from each of the three ISLs are graphed in Fig. 3. SM values for unsmoothed profiles were 0.92 (normal), 0.92 (diseasedCCl4), and 0.91 (diseasedAlc).
The ISL is nonlinear. Consequently, linear sensitivity studies are less informative and less meaningful than are location changes in lobule parameter space. Individually, the parameter changes in Fig. 4 did not cause statistically distinguishable changes in outflow profiles. In general, a 5% change in any one of the 12 adjusted parameters will produce an imperceptible change in an outflow profile and no change in SM value. However, a 5% change in all 12 parameters can cause a noticeable change in outflow profile.
Note that the diseasedAlc ISL was simpler than the diseasedCCl4 ISL in two ways. 1) Whereas 12 lobular parameter adjustments were needed to achieve diseasedCCl4 ISLs, only six were needed to achieve diseasedAlc ISLs. 2) Except for the C2BJumpProb adjustments, the magnitude of the adjustments needed for diseasedAlc ISL was smaller than those of diseasedCCl4 ISL. The fact that diseasedAlc ISLs exhibit fewer, smaller magnitude changes relative to normal ISLs is consistent with less profound and fewer observed pathological changes caused by alcohol pretreatment (Hung et al., 2002a,b).
Hereafter, all results are reported in the order of normal, diseasedCCl4, and diseasedAlc when values for all three are provided and in the order of diseasedCCl4 and diseasedAlc when only diseased ISL values are provided.
Changes in Lobule Properties at Three Levels Can Account for Cirrhosis-Caused Differences in Disposition.
cell density in spaces B and C were controlled by ECDensity and HepDensity. The value of ECDensity (0.65) was unchanged for diseasedAlc ISLs but was smaller (0.6) for diseasedCCl4 ISLs. The values of HepDensity reflect the same change: 0.7, 0.65, and 0.7 (Fig. 4). The latter difference maps to the observed CCl4-induced changes in hepatocyte quantitative ultrastructure (Hung et al., 2002b). Acceptable SM values were obtained without having to change the probability that a diltiazem released from an enzyme (in hepatocytes) would be a metabolite.
Four parameters control compound movement between spaces (Fig. 4, top row). Lower values relative to normal for the first two mean that disease reduced compound access to endothelial cells and hepatocytes. The magnitude and direction of those changes map to and are consistent with fibrotic changes (Hung et al., 2002a). The effective ability of space C to retain compounds (B2CJumpProb/C2BJumpProb) is a major determinant of the metabolic event rate. That ratio in diseasedAlc ISLs (2.6) was larger than that in normal (1.2) or diseasedCCl4 (0.62) ISLs. The larger ratio reduced the total number of diltiazems reaching CV thereby increasing metabolic events within space C. HepDensity's value was smaller in diseasedCCl4 (0.65) than in normal or diseasedAlc (0.70) ISLs and that contributed to decreased metabolic events and diminished diltiazem retention in space C.
The graphs in Fig. 5, A and B, show how stochastic compound movements within the three spaces influenced average compound resident times within lobules. The bar graphs specify the dose fraction having resident times within the indicated 10-second range. The three curves specify the fraction of dose having a lobule resident time equal to or less than the indicated time. Sucrose (Fig. 5A) did not enter cells. The bar graphs (a.1 and a.2) indicate that there were few differences in the spaces accessed and dwelt in by sucrose in normal and diseasedAlc ISLs. However, sucroses dwell times in the diseasedCCl4 ISLs were different in two ways. There was a dramatic reduction in sucroses having longer and also 0- to 10-second dwell times.
Dwell time patterns were quite different for diltiazem (Fig. 5B). For both diseased ISLs, there was a reduction in diltiazems having 0- to 10-second dwell times, and an increase in diltiazems having 10- to 20-second dwell times. The net effect was an increase in diltiazems having longer dwell times: disease makes it harder for diltiazem to move through and exit a diseased liver. Disease type altered dwell time patterns differently (b.2 and b.3). Surprisingly, there were no significant differences in mean (and SD) diltiazem resident times: 47.8 (26.0), 48.8 (25.2), and 45.9 (26.5) seconds. Several factors influenced the differences between resident times. Differences in A2BJumpProb between the three ISL types (0.38, 0.21, and 0.35) caused fewer diltiazems to move into space B in diseased relative to normal lobules. Similarly, differences in B2CJumpProb (0.55, 0.34, and 0.65) caused fewer of the diltiazems that did reach space B in diseasedCCl4 relative to diseasedAlc to move on to space C.
Changes in dwell time patterns coupled with changes in access to spaces A to C influence metabolism appreciably. The fraction of dose that was metabolized (that underwent a metabolic event) by a given time is graphed in Fig. 5C.
Tracing Compound Path Lengths and Spatiotemporal Binding Patterns.
The path length curves in Fig. 6, A and B, measure cumulative lengths of all SSs entered by each compound. All three ISL types had 70 SSs distributed among the three zones. A small subset of both diltiazem and sucrose path lengths in normal and diseasedAlc ISLs were long (second peak in Fig. 6, A and B), but they were essentially absent in the diseased ISLs. The SSs in zone 1, in combination with intrazone edges map to the interconnections between sinusoids that are most numerous in the periportal region, yet are absent in the perivenous region of normal lobules. Microscopy evidence suggests that some of those interconnections are lost in CCl4-treated cirrhotic lobules (Gaudio et al., 1997).
There are no wet-lab methods to measure which lobular subspaces are visited by a compound passing through the liver. We recorded each compound's traverse path (in grid spaces) for 100 seconds, until it either exited the lobule, was metabolized, or the run ended. Sucrose had shorter path lengths because it did not enter cells and so was more likely to reach CV before the run terminated. The shorter mean path for diseased ISLs shows that both disease types made it easier for compounds to move closer to CV as time advanced. It is evident from Fig. 6, A and B, insets that diseasedCCl4 ISLs had a more narrowly distributed variety of path lengths. Note also that the diseasedCCl4 ISLs had fewer of the shortest paths (0–25 grid spaces) than either normal or diseasedAlc ISLs.
Each component type that encountered and interacted with diltiazems and sucrose was recorded. The time diltiazems spent associated with components was in the order enzymes > binders > (unbound in) hepatocytes > (unbound in) endothelial cells. For all components, influential factors included the population density of diltiazems in each space, traverse lengths, and the values of these three parameters in Fig. 4: SoluteBindingProb, SoluteBindingCycles, and BindersPerCell. Despite those differences, the data in Fig. 6, C and D, show that the fraction of compounds that was in a lobule at a particular time and attached to a binder eventually reached a similar steady state ratio of approximately 0.8 in diseasedCCl4 and diseasedAlc lobules. However, the relative fractions bound in endothelial (space B) and hepatocyte layers (space C) were different.
Compound Dynamics within and between Zones.
Observing diltiazems percolating through individual SS within different zones provides an informative perspective. A video recording of compounds percolating through SSs in an earlier ISL is provided as a supplement by Yan et al. (2008b). Selected results are graphed in Fig. 7 and Supplemental Fig. S3 for one similarly sized SS in each of the three zones. Figure 7 shows the fraction of diltiazem dose within core and spaces A to C of the selected SS, regardless diltiazem's state or location. The differences in relative trends between zones are striking. Zone 2 peaks occur later than those in zone 1. The amount of diltiazem trickling through zone 3 was still increasing at 40 seconds for all three ISLs. The relative patterns for normal and diseasedAlc SS within the same zone are similar and clearly different from those in the diseasedCCl4 SS. More diltiazem was in spaces A and B of all three diseasedCCl4 SS compared with normal and diseasedAlc SS and that may map to fibrosis retarding diltiazem's access to spaces more distant from blood flow.
A portion of the data in Fig. 7 is shown subdivided further in Supplemental Fig. S1 into bound and unbound diltiazem within the endothelial and hepatocyte layers (spaces B and C). The dramatic increase in the binding of diltiazem within endothelial layers (keeping it away from the hepatocyte layer) maps to the observed fibrotic changes.
Tracing multiscale events facilitated identifying plausible micromechanistic differences regarding where and how normal and diseased ISLs interacted differently with diltiazem and sucrose. It is currently infeasible to obtain comparable wet-lab data, but advances in intravital microscopy methods are moving us closer. The outflow profiles in Fig. 3 achieved a degree of validation against wet-lab counterparts. That achievement allows us to conjecture that mappings exist between micromechanistic causes determining diltiazem disposition in the three types of ISLs (mappings 1 and 2 in Fig. 1) and corresponding causes in normal and diseased perfused livers. diltiazem and sucrose dynamics within and between ISL levels provide previously unavailable insight into plausible disposition details. However, the ISLs are not yet sufficiently refined to make precise predictions. With additional rounds of refinement and validation, their progeny can be expected to provide increasingly useful scientific predictions, as well as deeper insight into causal micromechanisms.
Similarity was achieved without having to use different normal and diseased sinusoid networks. However, it is premature to assign biological significance to that lack of difference. The three ISLs are specific examples of a number of similar ISLs that would also validate (analogous to how livers from matched rats are different, but for experimental purposes can be treated as being the same). To emphasize that we are early on the path to strong validation and trust, we have confined our rhetoric to describing and contrasting what occurred during ISL executions, and have drawn attention when those phenomena were consistent (or not) with reported wet-lab observations. Additional, more detailed observations on results are included in supplemental data.
The take-home message has five parts. 1) Multilevel ISL changes provide plausible mechanistic explanations for differences in compound disposition between normal and diseased livers. 2) Lobules changed at three levels: SS geometry; cell density and compound movement between SS spaces; and intracellular processes. 3) As illustrated in Fig. 8, those changes may have counterparts during disease progression. Twelve (of 29) ISL parameters were involved in normal-to-diseasedCCl4translation. However, six of those were invariant for normal-to-diseasedAlc translation. 4). Counterintuitively, three of the parameter changes required for translation changed in opposite directions. 5) Tracings of micromechanistic events below the validation level were required for a meaningful discussion of the results. So doing increased the heuristic value of ISL models.
All of the wet-lab data in Fig. 4 are coarse-grained, end-of-experiment, whole liver measures. Each ISL parameter, in contrast, influences fine-grained, spatiotemporal events occurring during simulation. Consequently, the two sets of data are not easily compared. For convenience, descriptions of all measures, along with conjectures on how ISL parameter values may relate to them, are provided under the discussion in supplemental data. None of the wet-lab measures in Fig. 4 were included as targeted attributes. There were two reasons. 1) We do not know whether or how the phenomena measured might actually influence the disposition of diltiazem or sucrose. 2) With the exception of extent of diltiazem metabolism, which can map to intrinsic clearance and P450, there were no measurable ISL features that map logically to the other six histopathology measures. Nevertheless, neither the wet-lab data in Fig. 4 nor other histopathology data reported in the original citations (Hung et al., 2002a,b) falsify any ISL parameter changes in Fig. 4.
ISL metabolic events map to intrinsic clearance. The ratio of metabolic events in Fig. 5 for diseased and normal ISLs maps to differences in intrinsic clearance for the different livers. No ISL parameter maps to P450. BindersPerCell does not; it maps to all cellular material that binds or sequesters diltiazem including the metabolizing enzymes. However, it is not surprising that wet-lab P450 and intrinsic clearance changes correlate. We can thus say that the ratio of metabolic events in diseased to normal ISL also maps to P450 changes. Because diseased ISL parameterizations were iteratively refined starting from the normal ISL, without regard to metabolism, the metabolism data in Fig. 5 can be taken as predictions that validated.
Assume that future data show that sinusoid networks are different between normal and diseased livers in specific ways. Inclusion of that data in the targeted attribute list would invalidate the current diseased ISLs. However, it would be straightforward to iteratively revise them until validation is again achieved. It is even feasible to automate that process.
Deep tracing of the type presented in Figs. 5 to 7 can present methodological challenges and problems because synthetic models are designed to realize an explicit separation between generators (micromechanisms) and phenomena. ISLs are completely observable, unlike livers. An asset is the ability to reach into the ISL and watch everything that happens. Doing so, however, may bias a researcher into trusting an ISL more than is warranted. For example, the ISL is relatively grounded in terms of runs, cycles, and steps of its various macro-, meso-, and microconstituents. The referent data against which we validated requires a very coarse measure, the outflow fraction profile. Hence, any measures, such as the tracings in Figs. 5 to 7 that are finer grained, although clearly relevant to the ISL, may not be valid in relation to the referent. It remains an untested hypothesis because we do not have data against which to validate.
The SS graph is not explicitly spatial. It is an abstract network of SS nodes that have no concrete mapping to the spatial structure of sinusoids within referent liver lobules. However, in aggregate it performs enough like a lobule to produce matching outflow profiles: quantifiable mappings from ISL to referent traverse paths may exist. Because outflow profiles are our validation data, there is no need to make lobules explicitly spatial. Adding such a requirement merely for the purpose of making the structure of the lobule more intuitive or “more like” the referent, absent any fine-grained spatial data against which to validate, runs counter to the parsimony guideline. There is some reasonable pressure to make the lobule spatially explicit, however. Lobular zonation data exist (e.g., Gebhardt, 1992; Jungermann, 1995; Oinonen and Lindros, 1998; Christoffels et al., 1999; Ohno et al., 2008) and are explicitly spatial. Because current lobules can contain graphs that cannot be projected onto a three-dimensional vector space, validation against such zonation data are difficult. Doing so would be more straightforward if the ISL lobules were spatially explicit. We are exploring strategies for making them so for future experiments.
We suggest that the incremental parameter changes that are necessary (and sufficient) to transform a normal into a diseased ISL (the changes on the right side of Fig. 1) may correspond abstractly to molecular, cellular, and sinusoid level transformations responsible for the pathogenesis from normal into diseased livers during CCl4 and alcohol treatment as illustrated in Fig. 8. The general consistency noted above between diseased ISLs during execution and the cited histopathology evidence supports the hypothesis. We thus have a tentative, yet promising, in silico model that enables us to visualize, abstractly from the perspective of diltiazem, how the consequences of cirrhosis may have progressed. That new capability represents an important step toward unraveling the complex influences of disease on drug disposition. Being able to transform one validated model into another is also important: it is evidence that the approach can facilitate rational translation of research results to useful applications and that may open doors to development of strategies for tailoring drug choices to help reverse disease conditions.
Having independently validated normal and diseased ISLs allows one to explore plausible drug disposition consequences of intermediate levels of disease and even disease that is more advanced. Because of individual differences in disease progression, conducting wet-lab experiments to document the former would be problematic, and the latter may be deemed unethical. By assuming disease progression corresponds to gradual change from normal to diseased ISL parameterizations, we can simulate a liver that has progressed halfway, for example, along the path to the currently documented disease state. We can project further parameter changes to explore plausible consequences of more advanced level of disease. Corresponding explorations of intermediate and advanced disease states would be infeasible using traditional inductive, physiologically based, mathematical models.
The methods and approach have been designed to enable the eventual development of horizontally and vertically integrated whole organism—whole patient—analogs: virtual organs, virtual patients. Once that goal has been achieved, it will be feasible before treatment to use in silico experimentation to anticipate narrow, plausible ranges of drug PK properties in patients with liver disease.
We thank members of the BioSystems Group for helpful discussion and commentary. S.P. and S.H.J.K. developed and implemented tracing methods. C.A.H. and G.E.P.R. designed and G.E.P.R. implemented and verified the ISL. C.A.H. designed experiments. S.P. and S.H.J.K. conducted experiments and analyzed results. C.A.H. developed and finalized the manuscript and figures using contributions from S.P., G.E.P.R., and S.H.J.K. C.A.H., G.E.P.R., S.H.J.K., and S.P. edited and revised the manuscript. M.S.R. proved the wet-lab data and editing suggestions.
Appendix: Descriptions of Key ISL Parameters in Table 1
monteCarloRuns. The number of lobule simulations averaged to produce results for one ISL experiment.
cycleLimit. Provides the simulation with a stopping criterion.
StepsPerCycle. Specifies the number of iterations the RefModel executes in a single ExperAgent cycle. The ExperAgent cycle is grounded to the default time scale in the validation data.
nominalProfile. Specifies which model (e.g., RefModel) to use as the nominal when calculating SM values.
GraphInputFile. Specifies the file to read if the SS graph is to be specified by an explicit graph (file format is GML).
GraphSpecFile. Provides the lobule graphical specification, and specifies the base file name (extension.ls) to be used if the graph is to be specified according to the lobule specification file.
Graphspeciterates. Tells the framework to modify lobule specification and run a Monte Carlo set (consisting of N runs) for each different lobule specification. Set to 1, it runs 1 set and provides 1 set of outputs. Set to 5, the first run uses the current contents of lobule.ls; it then runs 4 more sets, slightly modifying the Lobule specification each time, resulting in five sets.
SSTypeRatio = DirSinRatio/TortSinRatio. DirSinRatio specifies the percentage of SS that are type SA (“direct;” shorter and wider); TortSinRatio specifies the percentage of SS that are SB(“tortuous;” longer and thinner).
DirSinCirc. Circumference for each SAtype SS. Although fixed for this study, SASS circumference can be randomly drawn from a uniform distribution having the range DirSinCircMin and DirSinCircMax.
TortSinCirc. Circumference for each SBtype SS. Although fixed for this study, SBSS circumference can be randomly drawn from a uniform distribution having the range TortSinCircMin and TortSinCircMax.
DirSinLenAlpha, DirSinLenBeta, DirSinLenShift. These values set the parameters for a pseudorandom number generator using a modified form of the gamma distribution; the modification consists of a left-right shift of the distribution, allowing the user to clip off the front of the distribution. The pseudorandom number drawn is the length of a specific shorter, wider SS.
TortSinLenAlpha, TortSinLenBeta, TortSinLenShift. These values set the parameters for a pseudo-random number generator using a modified form of the gamma distribution; the modification consists of a left-right shift of the distribution, allowing the user to clip off the front of the distribution. The pseudorandom number drawn is the length of a specific longer, thinner SS.
A2BJumpProb. Specifies the probability that, when given the option, a compound will jump from space A to space B.
B2AJumpProb. Specifies the probability that, when given the option, compound will jump from space B to space A.
B2CJumpProb. Specifies the probability that, when given the option, Compound will jump from space B to space C.
C2BJumpProb. Specifies the probability that, when given the option, compound will jump from space C to space B.
CoreFlowRate. The number of slots (grid spaces) compounds in the SS core move forward during each step.
SinusoidTurbo. The PV to CV bias applied to the otherwise random walk for compounds in grid spaces of an SS. Smaller turbo means greater tendency of any one compound to wander sideways or backward. Larger turbo means a stronger flow from the input to the output of the SS.
ECDensity. Specifies the relative endothelial cell density; given the grid dimensions of space B of a given SS, it specifies the percentage of spaces that index an endothelial cell.
HepDensity. Specifies the relative hepatocyte density; given the grid dimensions of space C of a given SS, it specifies the percentage of spaces that index a hepatocyte.
BindersPerCell. In this study, the number of binders inside each cell. They are simple binders for endothelial cells and enzymes for hepatocytes. The ISL provides the option to randomly draw the number of binders for each cell from a uniform distribution having limits of BindersPerCellMin and BindersPerCellMax.
MetabolizeProb. Probability that an enzyme will transform diltiazem to a metabolite before it releases it.
SoluteBindingProb. Probability that, when a binder and compound make contact, the compound will be bound.
SoluteBindingCycles. Number of simulation cycles a binder holds a compound. This value maps to observing that the fraction bound within the same small region of a lobule is unchanged over an interval.
isMembraneCrossing. Specifies whether the compound can cross into (and out of) a cell or not; diltiazem can; sucrose cannot.
This work was supported by CDH Research Foundation [Grants CDHRP-06/07, -0014, -0023] (to S.P. and S.H.J.K.) and by Australian National Health and Medical Research Council [Grant 252871/9937600] (to M.S.R.).
Article, publication date, and citation information can be found at http://jpet.aspetjournals.org.
- in silico liver
- cytochrome P450
- central hepatic vein
- terminal portal vein tracts
- similarity measure
- sinusoidal segment.
- Received March 23, 2010.
- Accepted April 16, 2010.
- Copyright © 2010 by The American Society for Pharmacology and Experimental Therapeutics