Supplementary for: Modeling of Xenobiotic Transport and Metabolism in Virtual Hepatic Lobule Models

Computational models of normal liver function and xenobiotic induced liver damage are increasingly being used to interpret in vitro and in vivo data and as an approach to the de novo prediction of the liver’s response to xenobiotics. The microdosimetry (dose at the level of individual cells) of xenobiotics vary spatially within the liver because of both compound-independent and compound-dependent factors. In this paper, we build model liver lobules to investigate the interplay between vascular structure, blood flow and cellular transport that lead to regional variations in microdosimetry. We then compared simulation results obtained using this complex spatial model with a simpler linear pipe model of a sinusoid and a very simple single box model. We found that variations in diffusive transport, transporter-mediated transport and metabolism, coupled with complex liver sinusoid architecture and blood flow distribution, led to three essential patterns of xenobiotic exposure within the virtual liver lobule: (1) lobular-wise uniform, (2) radially varying and (3) both radially and azimuthally varying. We propose to use these essential patterns of exposure as a reference for selection of model representations when a computational study involves modeling detailed hepatic responses to xenobiotics.


Histological studies of mouse liver
Mice were housed in the AAALAC certified animal facility at Indiana University, Bloomington, Indiana. All animals were maintained in accordance to the Guide for the Care and Use of Laboratory Animals from the Institute of Laboratory Animal Research [1] and treatment protocols were approved by the Institutional Animal Care and Use Committee at Indiana University (IACUC)(permit number 14-012). Mice, after 12 hrs fasting, were treated with APAP at doses of 0 (control), 10, 100, or 250 mg/kg BW by IP injection. Mice were sacrificed 0.25 to 48 hrs after treatment in accordance with animal care procedures. A portion of the liver was fixed in 10% neutral buffered formalin and sectioned for histopathology analysis by hematoxylin and eosin (H&E) staining. The rest of the liver was snap frozen in liquid nitrogen and stored at -80C.
Images of H&E stained liver sections are shown below in S1 Fig 1. In a H&E stained liver section, cytoplasm appears pink and cell nucleus purple. In normal liver, hepatocytes have large round nuclei (S1 Fig   1A). In contrast, following treatment with a toxic dose of acetaminophen, hepatocytes near the central vein undergo necrosis as evidenced by a loss of nuclei (S1 Fig 1B).

Confocal imaging of rat and mouse liver sinusoids
Rats were maintained using protocols approved by the IACUC (permit number 10971). Rats were euthanized by inhalation of anesthetic gases (isoflurane) followed by bilateral pneumothorax, as well as exsanguinations as physical means to ensure death and as required by IUCAC. This minimizes animal distress and is quick, effective and efficient. This method is consistent with the recommendations of the AVMA Guidelines on Euthanasia.
A portion of a rat liver lobule sinusoid network is shown in S1 Fig 2. An animated file, giving better depth perception of this sinusoid network, is in file Rat liver sinusoids-1.mpg.
We labeled 300 µm thick sections of rat liver that had been perfusion fixed in 4% paraformaldehyde in PBS with FITC conjugated lens culinaris agglutinin (Vector Laboratories, Burlingame, CA) according to manufacturer instructions. Following labeling we optically cleared the tissue sections using benzyl alcohol-glycerol [2]. We acquired 3D mosaic images using a Leica SP2 MP microscope system (Leica Microsystems, Buffalo Grove, IL). We stitched these images, then segmented the vasculature and acquired tissue vasculature metrics using FIJI ImageJ [3,4].
The volume fraction of sinusoids in this rat lobule, excluding the CV and PT, is 10.4% ± 1.1%. Note that 15.3 ± 3.9% [5] and for rats 10.6 ± 4.5% [6], though it is not always clear if the number refers to the lumenal volume or the lumen+endothelial volume.
Ex vivo sinusoid networks show a small variation in sinusoids diameters, with the largest diameters occurring near the PT and CV and smaller diameters in the mid-lobular region. In addition, the ex vivo Supplement samples have a somewhat more densely connected network in three dimensions and it is not uncommon to observe nodes with four or more connections. We also measured the average sinusoid-hepatocyte interfacial area per unit parenchymal volume S int , which influences transport of xenobiotic in our model. S int , in this rat liver, was found to be 0.143 ± 0.019 µm 2 /µm 3 .

S1 Fig 2.
Rat liver lobule from confocal microscopy. Portion of a 3D rat liver lobule showing the sinusoid network. A central vein (CV) is located near the upper right portion of the image. The scale bar is 100µm and the typical sinusoid is ≈ 8µm diameter. The thickness of the imaged region was 50 µm.

Construction of the PIPE virtual liver lobule
The PIPE virtual liver lobule consists of six parallel arrays of hepatocytes along six parallel "pipe"-like sinusoids (S1 Fig 3). The virtual liver lobule has dimensions of 240 µm × 120 µm × 20 µm. PT and CV are not explicitly constructed in PIPE lobule. Instead, we used an extended "Wall" regions of the lobule to represent area that should be occupied by PT and CV (S1 Fig 3B). We modeled no diffusive or active transport processes between SINs and "Wall", and no metabolism within the "Wall". We applied periodic boundary conditions in two dimensions perpendicular to flow direction.

Construction of the NET virtual liver lobule
In the Main Text, we described the steps of sinusoids network construction including determination of network node locations and creation of network edges. Here, we use a schematic to illustrate the processes (S1 Fig 4). We started by constructing the hexagonal boundary of a CV-centered virtual liver lobule including six PTs and peripheral sinusoids. We then added the central vein at the center of hexagon and drew virtual concentric rings surrounding central vein with radial spacing of one hepatocyte size (S1 Fig 4A). Next, on each virtual ring, we randomly selected a point as a node location and determined other node locations by "walking" on the ring with a step of one HEP size. We also subdivided each of six peripheral sinusoids by unit of one HEP size to determine perimeter node locations (S1 Fig 4B). Then, we started from the innermost ring and connected each node on the ring to the nearest node on the next outer ring to form an edge. Nodes on the the outermost ring connected to the nearest nodes on peripheral sinusoids or to PT nodes (S1 Fig 4C).
Lastly, we determined locations of cell centers of HEPs by selecting the middle point between any two nodes on the same virtual rings, which will later become the source voxels of HEPs in CC3D (S1 Fig 4D).
After creation of the sinusoid network, we converted the model into a 3D spatial model in CC3D. PTs and the CV are converted into pixelated cylinders as are the sinusoid segments. The sinusoids have a diameter of 8µm. The thickness of the virtual liver lobule slice is the thickness of one hepatocyte, 20 µm. Individual heptocytes developed from the one-voxel cell centers created in the layout algorithm, to fill the remaining space between sinusoid network. A typical NET virtual liver lobule in the model is shown in S1   passive transport, active transport and metabolism. This section presents a more detailed schematic of these processes (S1 Fig 6). The spatial configuration of hepatocytes and sinusoids in a zoomed view of a virtual liver lobule is shown in S1 Fig 6A. The passive (diffusive) transport of the xenobiotic is modeled bidirectionally at the interface between sinusoid blocks (SINs) and hepatocytes (HEPs) and at the interface between neighboring HEPs. Active transport (transporter mediated) of the xenobiotic is modeled unidirectionally at the interface between SINs and HEPs, from the former to the latter. Metabolism of xenobiotic is modeled within individual HEPs. The topological configuration of conveyor belt blocks (CBs) and junctions associated with the sinusoids network are shown in S1 Fig 6B. Advection of blood-borne xenobiotic is modeled unidirectionally from an upstream CB to the downstream one on the same segment. Converging segments receive the total mass of the xenobiotic in the last CBs of the upstream segments (e.g. segment ik and segment jk) into the first CB of the downstream segment (e.g. segment kr).
SINs, the structural elements of sinusoid segments, have equal length; CBs, the functional element of xenobiotic advection, have different sizes on different segments proportional to blood flow velocity in that S1 Fig 5. A NET virtual liver lobule. A NET virtual liver lobule constructed in CC3D following procedures described in S1 Fig 4. The image on the left is a 2D cross-sectional view through the plane of sinusoid network. The image on the right is a partial cutaway of the 3D lobule slice. The hexagonal lobule has an edge size of 100 px corresponding to 200 µm and a thickness of 10 px corresponding to 20 µm.
sinusoid segment. Those CBs on a segment that carry blood flow with greater flow velocity have larger size and consequently the operation of "conveying" the xenobiotic one block downstream reflects faster "advection" (proportional to linear blood flow velocity).
To coordinate the advection part of the model, which is calculated at the CB level, with passive/active transport and metabolism, which are calculated at the SIN/HEP level, the model converts the mass of the xenobiotic between SINs and their associated CBs. Each SIN is associated with one or more CBs, which is guaranteed by selecting a sufficiently small time step ∆t. In a simulation cycle of advection, transport and metabolism, mass conversion between CBs and SIN occurs twice. Following advection step, the model adds up the masses of the xenobiotic in the associated CBs to give the mass of the xenobiotic in a particular SIN.
Then the model calculates passive/active transport among SINs and HEPs and metabolism within HEPs.
Following metabolism step, the model adds (or subtracts) the same percentage change of the xenobiotic mass in a particular HEP to (or from) its associated CBs. In the next simulation cycle, all CBs have updated xenobiotic mass that reflects the post-transport level in the blood.

Comparison of Flow Calculations in CC3D with NetFlow of Priese and Secomb
To calculate blood flow rates in the virtual sinusoid network, we wrote a Python module in CompuCell3D (CC3D) [7]. To check its accuracy, we compared calculated flow rates with results using an existing tool NetFlow developed by Pries, Secomb and coworkers for blood flow calculations in capillary (sinusoid) diameter vessels [8,9]. We found that the correlation between our results in CC3D and those from NetFlow was very good (S1 Fig 7). The dispersion about the regression line is due to the inclusion of hematocrit variation in the NetFlow calculation, whereas the CC3D calculation does not take hematocrit into account.
The correlation between the NetFlow and CC3D results is essentially a perfect if the NetFlow calculation is carried out with a hematocrit of zero. In addition to using NetFlow to verify our flow calculations we also used it to explore the hematocrit partitioning across the representative lobule. As shown in S1  [8,9]

Structure
To explore the effects of anastomotic patterns of sinusoid network on structural and hemodynamic quantities of the virtual liver lobule as well as hepatic exposure to the xenobiotic, we examine the structural and hemodynamic features of six total lobule slices created using the same construction algorithm but with different anastomotic patterns (S1 Table 1 and S1 Fig 9). All of these constructed liver lobules show similar results in terms of mean hepatocyte volumeV hep , sinusoids volume fraction η sin , average sinusoid-parenchymal interfacial area per unit parenchymal volume S int , mean blood flow velocityv, and total lobular volume inflow rateQ. Virtual lobules with different anastomotic patterns show appreciable, though small, difference in hepatic exposure to the xenobiotic, as reflected by different lobular-wise average xenobiotic concentrationc lob and two descriptors of spatial variation in xenobiotic micro doses γ cp and γ af . S1  and CL int = 0.1/s for comparison of xenobiotic exposure. Quantitative results from these simulations are summarized in S1 Table 1.

Flow Calculations and Microdosimetry Simulations with Varied Lobule Size
To explore the effects of lobule size on structural and hemodynamic quantities of the virtual liver lobule, we construct a lobule with a larger PT-to-CV distance of 400 µm. Overall, the descriptors of structural features are very consistent between the default and the larger-size virtual liver lobules, includingV hep , η sin and S int (S1 Table 2 Spatial pattern of xenobiotic showed some difference as well in an example simulation with D = 1 × 10 −6 cm 2 /s, α = 10, β = 10/s and CL int 0.1/s respectively. In the simulation with a 400 µm lobule, the xenobiotic was even more predominantly absorbed and metabolized in the periportal region, which led to a very small γ cp and very large γ af .

Flow Calculations and Microdosimetry Simulations with Varied PT-CV Pressure Drop
We also investigate the effects of larger a PT-to-CV pressure difference on hydrodynamic features of the virtual liver lobule as well as hepatic exposure to the xenobiotic (S1 Table 3). As discussed in the main text, we justified the selection of ∆P = 0.6 mmHg as PT-to-CV pressure difference in the virtual mouse liver lobule by comparing calculated average flow velocity and total inflow/outflow volume rate with experimental values. However, it's interesting to examine if the virtual liver lobule with the same architecture (therefore S1  the same flow resistance) can represent a rat liver lobule by using a greater PT-to-CV pressure difference.
Not surprisingly, both average flow velocity and volumetric flow rates increase proportionately to increase in ∆P (S1 Table 3). Interestingly, a PT-to-CV pressure difference of ∆P = 0.6 mmHg (reported value for rat liver lobule) results in mean flow velocity ofv = 337.62 µm/s, within the range of reported experimental

19/45
Supplement values [10][11][12][13]. In addition, flow distribution is essentially identical if we scale the velocities by pressure differences (images not shown), because the flow resistances at individual sinusoids segment are the same in these simulations. Therefore, both average flow velocity and volumetric flow rates increase proportionately to change in ∆P .
In terms of hepatic exposure to the xenobiotic, we ran an example simulation with D = 1 × 10 −6 cm 2 /s, α = 10, β = 10/s and CL int = 0.1/s. Simulations with greater ∆P show smaller lobular-wise average xenobiotic concentrations and more uniform lobular distribution of xenobiotic, as both radial and azimuthal discrepancies are closer to unity, indicating smaller hepatic extraction of the xenobiotic due to shorter blood transit time. S1 In the main text, we discussed the effect of varying active transport parameters α and β on radial and azimuthal distribution of the xenobiotic within a NET lobule in PULSE simulations (S1 Fig 17). In particular, we observed greater peri-central xenobiotic concentration when α and β were small. Here, we ask if this pattern found in NET lobule was also present PIPE lobule. In contrast to NET lobule, PULSE simulations within the PIPE lobule with small α and small β showed lobular-wise uniform distribution of xenobiotic (S1 Fig 11). This striking comparison emphasized the interpretation in the main text that heterogeneous pattern of blood flow and diverse blood flow transit times contributed to observed preferential peri-central distribution of the xenobiotic within the NET lobule in PULSE simulations.

Hepatic NET Exposure to PULSE Xenobiotic
To explore the simplest case of xenobiotic uptake across the simulated lobule we examined the effect of varying the passive transport rate constant, in the absence of active transport and metabolism, for a narrow pulse (∆t = 0.1s) of xenobiotic. In passive transport-only PULSE simulations, xenobiotic was absorbed into PLOS 20/45 Supplement S1 Fig 11. Heat maps of central-to-peripheral ratios within the PIPE lobule at t = 20 s in active transport-only PULSE simulations. (A) Central-to-peripheral ratios γ cp for simulations with 12 different (α, β) pairs. γ cp defines the ratio of the average concentration within radially the most peri-central zone to the average concentration within the most peripheral zone. Warmer color reflects greater peri-central than peripheral concentration. "PP" and "PC" refer to peri-portal and peri-central, respectively. "P 1 "-"P 4 " correspond to "N 1 "-"N 4 " (see S1 Fig 17) in terms of (α, β) values. (B) Spatial snapshots of xenobiotic concentration within the PIPE lobule at 20 seconds in simulations corresponding to parameter domains marked "P 1 "-"P 4 ".
the hepatocytes and eventually cleared from the cells back to the blood (S1 Fig 12). It should be noted that in a particular simulation the same diffusive trans-membrane rate constant was applied to both cell-blood and cell-cell transfers, but this did not guarantee the same diffusive flux, which depends also on variables such as interface area and blood flow velocity.
For a small diffusive trans-membrane rate constant D = 10 −7 cm 2 /s (top row S1 Fig 12) extraction of xenobiotic from the blood is insignificant and most of the pulse exited via the central vein with little extraction into the parenchymal space. We observed penetrating "spike" patterns of high blood concentration of xenobiotic around the portal triads at early time points (t = 1s) due to the higher axial blood flow, particularly with small passive rate constants. Slow passive transport resulted in overall low levels of xenobiotic uptake but also reduced the rate that absorbed material was washed back out of the hepatocytes after the end of the input pulse (S1 Fig 12 upper right panel).

21/45
Supplement of xenobiotic was notably enhanced compared to the slow passive transport situation. The majority of xenobiotic diffused into hepatocytes before eventually being cleared out and carried away by blood, evidenced by high peripheral concentration at t = 0.5s and 1s and a well-maintained radial gradient at t = 1s. It is also clear that the diffusive flux between adjacent hepatocytes was insufficient to effectively reduce the gradient. The "spike" patterns were also observed early in the simulation but faded out as the concentration equilibrated between blood and hepatocytes, which ultimately contained xenobiotic of approximately equal concentration during the wash-out phase.
For large diffusive trans-membrane rate constant, D = 10 −5 cm 2 /s (bottom row S1 Fig12), uptake of xenobiotic dominated the advection in blood, indicated by rapid "closure" of the zero-concentration "hole" in the central region via cell-to-cell diffusive transport. Remarkably, unlike simulations with smaller diffusive rate constants, this simulation showed earlier initiation of xenobiotic clearance, likely before t = 1s, which implies that the interval between first entrance and first exit could be shorter than the blood transit time at very high cell-to-cell diffusion rates. Furthermore, the "spike" patterns disappear quickly, suggesting that variability in transit time imposed by the network structure had little influence at high diffusive flux. It appears that blood flow was the limiting factor in this situation, since little change in xenobiotic distribution was observed between t = 2s and t = 5s.
One characteristic that differentiates the PIPE model from the NET model is the possibility of non-uniform azimuthal distributions of xenobiotic in the latter case. Azimuthal distribution have an angular dependence, often corresponding to the 6-fold angular symmetry of an ideal liver lobule. As mentioned earlier, the NET layout has variable path lengths for PT (axial) to CV paths opposed to facial to CV paths.
The differences in path lengths affect both the blood flow velocity and the number of hepatocytes that a particular flow path passes. To explore the radial versus azimuthal distributions we divided the NET lobule into 8 equally spaced radial zones and 12 azimuthal zones, each 30 degrees wide. Radial and azimuthal distributions of xenobiotic at t = 0.5s and t = 20s within the lobule are shown in S1 Fig 13 and S1 Fig 14, respectively. The first time point was selected to explore the uptake phase, while the latter to explore the washout phase. At t = 0.5s, a smaller diffusive trans-membrane rate constant resulted in a steeper radial gradient and more axially biased azimuthal distributions. In addition, mean xenobiotic concentrations within radially or azimuthally divided zones were lower in simulations with a smaller diffusive rate constant, indicating lower xenobiotic extraction (S1 Figs 13A, 13C). At t = 20s, simulations with an intermediate diffusive trans-membrane rate constant displayed the largest reversed radial gradient of xenobiotic (S1 Fig   14A), because at the peripheral region, hepatocytes and inflow xenobiotic-free blood had the greatest concentration difference and these peripheral cells washed out more quickly. As blood became increasingly visually intact leading pulses in the sinusoids channels at t = 0.5, 1, 2s (top two rows in S1 Fig 15). The "spike" patterns, introduced by distinct transit times across the lobule, were evident at t = 1s. Interestingly, although hepatocytes along the six porto-central axes seemed to experience pulses of xenobiotic sooner, they generally failed to absorb more, as indicated by lower concentrations (cooler color) along these axes at t = 20s. Also strikingly and somewhat counter-intuitively, peri-central hepatocytes, on average, accumulated more xenobiotic than peripheral hepatocytes. Both phenomena can be interpreted as a consequence of lobule-wide small uptake rate but different cumulative uptake durations. In PULSE simulations where the pulse size is as narrow as ∆T = 0.1s, peri-central hepatocytes have the advantage of encountering a greater number of derived pulses split from the original incoming pulse due to variable transit times imposed by the network structure. In contrast, peri-portal cells experience only a single, short duration pulse. As a result, peri-central cells tended to accumulate more xenobiotic.
Axial and facial viability in xenobiotic absorption resulted from the distinct local flow velocities. Facial hepatocytes, despite fewer encounters with pulses compared to axial cells, had significant duration of xenobiotic in the local blood with some pulses persisting for at least 5 seconds due to the slower blood velocity (t = 5s top two rows in S1 Fig 15). Apparently, increasing pulse width gradually diminishes the relative difference in uptake durations between peri-portal and peri-central hepatocytes, as will be discussed in more detail later. For α = 1 and β = 1/s, the distribution was much different and the pattern was inverted. Peripheral hepatocytes absorbed more xenobiotic than peri-central cells. In the periphery, axial cells accumulated more xenobiotic than facial cells did. For α = 10 and β = 10/s, the distribution was more extreme, and xenobiotic was preferentially absorbed into peripheral, especially peri-portal, hepatocytes.
In order to quantitatively depict the regional accumulation of the xenobiotic in active transport only predominantly in peri-central hepatocytes for smaller αs (0.01 and 0.1) but in peripheral area for larger α (1 and 10) (S1 Fig 16A-B). In terms of the azimuthal distribution, the xenobiotic primarily accumulated in cells along facial paths for smaller α (0.01 and 0.1) but along porto-central pathways for larger α (1 and 10) (S1 Fig 16C-D). Therefore, in PULSE simulations with β = 1/s, a critical transition between two patterns We also employed two descriptors to characterize radial and azimuthal distributions of xenobiotic in the active transport-only simulations at t = 20s, as illustrated using heat maps in S1 Fig 17. The descriptor for radial distribution, termed the central-to-peripheral ratio γ cp , defines the ratio of the average xenobiotic concentration at the peri-central-most radial zone to that at the peripheral-most radial zone (S1 Fig 17A).
Scenarios with γ cp > 1.25, 0.75 < γ cp ≤ 1.25 and γ cp ≤ 0.75 reflect preferentially peri-central, radially uniform and preferentially peripheral exposure to the xenobiotic, respectively. The descriptor for azimuthal distribution, termed the axial-to-facial ratio γ af , defines the ratio of the average xenobiotic concentration at the six axial azimuthal zones to that at the six facial azimuthal zones (S1 Fig 17B). Scenarios with γ af > 1.25, 0.75 < γ af ≤ 1.25 and γ af ≤ 0.75 reflect preferentially axial, azimuthally uniform and preferentially facial exposure to the xenobiotic, respectively.
For example, for the simulation with α = 1 and β = 1, the rightmost data point of the yellow curve divided by the leftmost data point of the same curve in S1 Fig 16A gave the central-to-peripheral ratio γ cp for this (α, β) pair. The mean value of the six data points of the yellow curve overlapping with dashed lines divided by the mean value of the other six data points of the same curve in S1 Fig 16C gave the axial-to-facial ratio γ af . In S1 Fig 17A-B, the vertical axis from top to bottom reflects increasing maximum uptake rate and the horizontal axis from left to right reflects increasing half-saturating concentration of active transport. Both central-to-peripheral and axial-to-facial heat maps showed a "two-domain" pattern with the boundary between the two domains being the secondary diagonal of the (α, β) parameter matrices (S1 Fig 17A-B). This result echoes the xenobiotic distribution in select simulations as depicted in S1 Fig 16. Very small α alone, or combination of small α and small β, resulted in greater peri-central concentration radially and slightly greater facial concentration azimuthally. In contrast, simulations with large α and large β showed greater peripheral accumulation radially and greater axial accumulation azimuthally. An alternative description of the latter case is simply greater peri-portal accumulation. S1 Fig 17C shows corresponding spatial snapshots of xenobiotic distribution for four select domains "N 1 "-"N 4 ".
In contrast to NET simulations, PIPE simulations didn't show higher peri-central xenobiotic concentrations in any combination of α and β (Supplement 2.5 ). In addition, As briefly mentioned earlier, we also examined whether (α, β)-domain specific pre-centrally preferential xenobiotic distribution was preserved when the pulse size was increased, as detailed in Supplement 2.6. We found that an increase in PULSE size alone was sufficient to alter the spatial pattern observed in active transport only PULSE simulations with weak active transport, from predominantly centri-lobular to lobular-wide uniform xenobiotic distribution. While this section focuses on describing spatial patterns of hepatic exposure to the xenobiotic, Supplement 2.7 further discusses the overall hepatic exposure as reflected by lobular partition of the

29/45
Supplement xenobiotic in the hepatocytes versus in the blood within the virtual liver lobule. S1 Fig 17. Heat maps of central-to-peripheral and axial-to-facial ratios at t = 20 s in active transport-only PULSE simulations. (A) Central-to-peripheral ratios γ cp for simulations with 12 different (α, β) pairs. γ cp defines the ratio of the average concentration within radially the most peri-central zone to the average concentration within the most peripheral zone. Warmer color reflects greater peri-central than peripheral concentration. "PP" and "PC" refer to peri-portal and peri-central, respectively. (B) Axial-to-facial ratios γ af for simulations with 12 different (α, β) pairs. γ af defines the ratio of the average xenobiotic concentration within azimuthally the six axial zones by the average concentration within the six facial zones. Warmer color reflects greater axial than facial concentration. "FA" and "AX" refer to facial and axial, respectively. Blocks marked "N 1 "-"N 4 " in (A) and (B) correspond to (α, β) pairs used for simulations in S1 Fig 15. (C) Spatial snapshots of xenobiotic concentration at 20 seconds in simulations corresponding to parameter domains marked "N 1 "-"N 4 ". Refer to S1 Fig 15 for color bars.

Descriptors of Hepatic Exposure within a NET lobule in Active Transport Only PULSE simulations with Varying Pulse Sizes
In the main text, we discussed spatial patterns of xenobiotic exposure in active transport only PULSE simulations with a PULSE size of ∆T = 0.1s (S1 Fig 17). A striking feature of the spatial pattern observed for weak active transport and short pulse duration was that xenobiotic preferentially accumulated in centri-lobular region. We interpreted this phenomenon as a consequence of variation in cumulative uptake durations of cells located in different lobular regions, and postulate that an increase of PULSE size would diminish the difference in cumulative uptake durations and thus alter radial variation in hepatic exposure to the xenobiotic.
We tested this prediction by gradually increasing the pulse size from 0.2s to 10s in PULSE simulations and created heat maps of the radial descriptor γ CP (S1 Fig 18). As predicted, with pulse size increasing from 0.2s to 10s, (α, β) domains of active transport that led to higher peri-central xenobiotic concentration gradually regressed to a smaller size at the top left corner of heat maps, and displayed smaller magnitudes of discrepancy. Interestingly, the dark red color in the heat maps, that represents a large discrepancy in hepatic exposure to the xenobiotic where peri-central xenobiotic accumulation is at least twice as large as peripheral accumulation, disappeared as the pulse size increased from 1s to 2s, which coincided with typical blood transit time of 1s ∼ 2s (for example, compare t = 1s and t = 2s snapshots at top row in S1 Fig 15). These results implied that in active transport-only PULSE simulations with lobular-wise weak active transport, a peri-central hepatocyte likely accumulated at least twice as much xenobiotic as a peripheral cell did on average in the virtual liver lobule, if the pulse size of the inflow xenobiotic was narrower than the typical transit time of blood flow in the simulated sinusoid network.
PLOS 31/45 S1 Fig 18. Heat maps of central-to-peripheral ratios within a NET lobule at t = 20 s in active transport-only PULSE simulations with different pulse sizes. Sizes of the pulse, ∆T p , indicated below each heat map, range from 0.2 to 10 seconds. Refer to S1 Fig 17 for details on the range of (α, β) pairs. "PP" and "PC" refer to peri-portal and peri-central, respectively.

Lobular Partition of the Xenobiotic in PULSE Simulations
The main text discusses the spatial patterns of hepatic exposure to the xenobiotic via simulations of microdosimetry of xenobiotics with diverse strengths of transport and metabolism. These simulations also produced relative quantities of the xenobiotic in the hepatocytes versus in the blood within the virtual liver lobule, which was an indicator of overall lobular exposure to the xenobiotic.
For passive transport only PULSE simulations, xenobiotic quantity partitioned in the parenchyma had a two-phase temporal pattern: uptake phase and clearance phase (S1 Fig 19A). Simulation with larger diffusive trans-membrane rate constant D displayed faster uptake and typically slower clearance. An obvious exception was the simulation with smallest diffusion coefficient, which showed slowest uptake and slowest clearance. Another feature of these curves was peak magnitudes observed for different diffusion coefficients.
The largest three diffusion coefficients seemed to have the same peak magnitudes despite different peak times, suggesting identical equilibrium states of blood-parenchymal partition of xenobiotic. For the two smallest diffusion coefficients, peak magnitude was lower due to weaker diffusive transport.
For active transport only PULSE simulations, only inward uptake was permissive and all simulations arrived at and maintained at a plateau of xenobiotic concentration within hepatocytes (S1 Fig 19B-D). For a particular β, larger α caused greater extraction of the xenobiotic (S1 Fig 19B). For a particular α, larger β resulted in greater extraction and earlier equilibration of xenobiotic (S1 Fig 19C-D). For extremely large α and β, all incoming xenobiotic ended up within the hepatocytes (S1 Fig 19D). PERSISTENT simulations showed strong resemblance to active transport only PULSE simulations except that weak active transport didn't cause predominant distribution of xenobiotic in the former case (S1 Fig 20).
These PERSISTENT simulations were effectively the extreme cases of simulations with varying pulse sizes, as discussed above (S1 Fig 18). To echo the conclusions above, these simulations suggested that preferential centrilobular accumulation of xenobiotic occurred only if pulse size was so narrow that centrilobular cells would acquire cumulatively longer duration of exposure to xenobiotic due to spatial arrangement of hepatocytes and vessels within the NET lobule. Increasing pulse size in PULSE simulation or introducing xenobiotic continuously in PERSISTENT simulation could partially or entirely eliminate this radial discrepancy.

Variability in Calculated Hepatic Exposure from Replicate Simulations
Simulation of xenobiotic transport and metabolism is stochastic in the model because in each time step the processing order of the cells is random. To test the model's sensitivity to this stochastic component, we ran 50 replicates using the same lobule configuration and examined the variability in lobular-wise mean concentration of the xenobiotic (c lob ) and central-to-peripheral ratio (γ cp ) at 20 seconds using select parameter sets from PERSISTENT simulations. These three parameter sets were; "simA" which had only active transport with α = 1 and β = 1/s; "simAP" which had active/passive transport with α = 1, β = 1/s and D = 10 −6 cm 2 /s; "simAPM" which had active/passive transport and metabolism with α = 1, β = 1/s, D = 10 −6 cm 2 /s and CL int = 1/s.
Average values and standard deviations from the 50 replicate simulations for the three parameter sets are listed in S1

Lobular Partitioning of the Xenobiotic in PERSISTENT Simulations
In the main text, we discussed how the interplay between active/passive transport and metabolism affected spatial patterns of hepatic exposure to the xenobiotic. We continue to discuss here their effects on lobular partition of the xenobiotic. One marked difference in temporal pattern of xenobiotic partition in these PERSISTENT simulations compared to the PULSE simulations is the linearly incremental dotted lines that represent continuous inflow of xenobiotic with constant concentration.
In simulations without metabolism (S1 Fig 21), for particular β and D, larger α resulted in greater xenobiotic partition in the parenchyma at a given time point. For particular α and D, larger β resulted in greater xenobiotic partition in the parenchyma at a given time point. For particular α and β, larger D resulted in smaller xenobiotic partition in the parenchyma at a given time point. S1 Fig 22 and S1 Fig 23 depict

Animations
We describe below parameter sets used for animations (S1 Animations) listed in S1 Table 5. We named the animation files by concatenating simulation type, structure type and process type(s). For instance,

Source Files
The model is implemented in CompuCell3D package (http://compucell3d.org/). Tutorials for using CompuCell3D can be downloaded from the website. We include in Supplementary Information the source files for a NET simulation with parameter values α = 1, β = 1/s, D = 10 −6 cm 2 /s and CL int = 1/s as a compressed file named S1 Code (see the subfolder "[CODE] alpha1 beta1 D1e-6 CLint1.zip"). A list of files in this compressed file are briefly described in S1 Table 6. All simulations described in the Main Text can be derived by adjusting parameter values in "pDictLiverXeno.csv" as described in S1 Table 7. Source files for other simulations appearing only in Supplementary Information are available upon request. S1 Table 6. Files included in the zip file S1 Code.