A physiologically-based flow network model for hepatic drug elimination I: regular lattice lobule model

We develop a physiologically-based lattice model for the transport and metabolism of drugs in the functional unit of the liver, called the lobule. In contrast to earlier studies, we have emphasized the dominant role of convection in well-vascularized tissue with a given structure. Estimates of convective, diffusive and reaction contributions are given. We have compared drug concentration levels observed exiting the lobule with their predicted detailed distribution inside the lobule, assuming that most often the former is accessible information while the latter is not.


Background
The liver is the major site of biotransformation of endogenous and xenobiotic substances including drugs in the body. Its main role is to prevent accumulation of a wide range of chemical compounds in the blood by converting them into a form suitable for elimination. Such vital processes, however, can potentially damage liver tissue and hence its functionality. Examinations of hepatic clearance shows that substance extraction not only can be limited by damaged hepatocytes (liver cells) but also by the intrinsic (enzymatic) ability to eliminate the drug, by the resistance to drug transport from blood to eliminating tissue cells, and by the hepatic blood flow itself. Indeed, perturbations in the hepatic flow patterns e.g. through disease or aging can significantly alter the systemic clearance of these substances. As a result, a quantitative understanding between the liver performance and its structural integrity would be a great utility in the toxicology of newly developed drugs.
To quantitatively assess these interacting processes, numerous models of hepatic clearance have been developed both in vitro and in silico, see an extensive review in [1]. Although in vitro liver models can be considered as the best alternative to analyze actual organ performance, their major shortcoming is to assume the liver as a homogeneous environment. It is known that metabolic functions, such as xenobiotic metabolism, amino acid conversion, cholesterol synthesis, oxidation, etc. are all zone specific. As a result, in vitro models are not necessarily well suited for studying and capturing the physiology of liver with hepatocyte morphology, property and functionality varying from point to point throughout of the organ. Because of this, there has been surge of interest among investigators in developing computational (in silico) models of the liver in order to study inhomogeneous structures in the presence of xenobiotic metabolisms and their effects on the hepatic clearance. A major challenge remains to integrate both physiological aspects and various metabolic processes in a single model.
Generally, in all computational models, nutrient and metabolite containing blood enters through a segment of interest (the lobule) via portal veins (and hepatic arteries) and flows through small channels (sinusoids) that circumvent liver cells (hepatocytes), Figure 1a. These (possibly toxic) chemicals transport from sinusoids to tissue, are converted to other forms by hepatocyte metabolic reactions, transport back to the sinusoids, and are finally excreted from the lobule via the central (hepatic) vein.
Multiple computational models have been proposed: single and multi-compartment models, distribution-based models, agent-based models, interconnected parallel tubes, etc. See [1] for further details. Two simple and commonly employed models are (a) the well-stirred compartment model, and (b) the parallel tube model. The first considers the liver as a homogeneous compartment [2] and the concentration of drug in the liver to be in equilibrium with venous (emerging) blood drug concentration. The second, developed by Bass and coworkers [3,4], regards the liver as a series of parallel tubes and the concentration of drug declines along the length of these tubes, due to enzymes distributed evenly along the length of these tubes. Pang and Rowland [5,6] present an in-depth comparison of the strengths and limitations of each approach. Subsequent developments [7][8][9][10][11] indicate that these two models represent limiting cases, and attempts to generalize each, e.g. by adding a series of well stirred compartments or by considering a distribution of tubes or adding directly dispersive mixing, have tended to make the models more consistent with each other. Further enhancements of these models include consideration of heterogeneity by adding stochastic processes, Michaels-Menten kinetics to describe intrinsic elimination, uniform and non-uniform metabolic processes, and enzyme zonation along the sinusoid by stacking a number of compartments, each with different concentration level and metabolic activities [1]. We emphasize, however, all of the above models are essentially 1D or pseudo 1D approaches and hence ignore the true three dimensional structure of the liver.
In series of papers, we intend to introduce a physiologically-based lattice model of the functional unit of the liver that integrates all structural and physiological aspects discussed above. Our model takes into account parameters such as the distribution volume, permeability, and porosity of the liver vasculature and cells. Furthermore, it includes both flow-limited and diffusion-limited exchange of drug molecules into the extravascular space from the sinusoids; sequestration of the drug molecules within liver cells with enzymatic transformation; and exchange of the metabolized drug molecules back from liver cells to the vasculature. Estimates and consequences of the competing flow processes are given. Furthermore, the enzymatic transformation of the drug can be either simple or saturable. The model allows us to include the effects of the intrahepatic mixing process on the enzymatic transformation of drug molecules at the cellular level.
In this paper, we focus on a simple regular lattice to define a base-case model for our study. Here we explore the dynamics of competing convective, diffusive, and reactive processes acting on an injected drug. Multiple sensitivity simulations are performed and their consequences on drug concentration levels found exiting the lobule as well as their detailed spatial distribution within the lobule are discussed. Future extensions to additional structural and/or physiological inhomogeneities such as non-regular and statistical lattices, enzymatic zonation, and extensions to 2D / 3D hexagonal lattices will be developed, see Figure 1c and 1d.

Liver architecture
At the macroscopic level, the liver consists of three vascular trees, two supply trees that originate from the portal vein and hepatic artery, and one collecting tree that drains into the portal artery [12]. The vessels bifurcate down to the terminal arterioles and venules, which are organized into portal tracts along with a terminal bile duct. Liver cells, called hepatocytes, radiate outward from the terminal vessels. These plates of hepatocytes are interspersed by sinusoids, which play the role of the capillary in the liver, and the spaces of Disse, which are the extravascular space of the liver [13]. Finally, the blood is collected and removed by the hepatic venules.

Functional unit
The functional unit of an organ is the smallest structural unit that can independently serve all of the organ's functions [14]. Because of its complexity, there is continued debate about what the functional unit of the liver should be. The classic lobule is a hexagonal cylinder, centered around a hepatic venule and with portal tracts situated at the corners. The portal lobule has a similar shape but is centered about a portal tract with the hepatic venules at the periphery [15]. The acinus is another proposed unit and is based on the pattern formed by the cords of hepatocytes between two central venules. Matsumoto and Kawakami [16] suggested that the classic lobule can be divided into primary lobules, which are cone-shaped and each fed by one portal tract and drained by one hepatic venule. Teutsch and colleagues [17,18] performed a morphological study of rat and human liver lobules, and their results support the idea of a secondary unit made up of primary units in what they term as a modular architecture. They conclude, however, that the primary unit is more polyhedral in shape than conical. Other experiments done by Ruijter et al. [19] suggest that the primary unit is needle-shaped and that there are equal amounts of portal and central vein associated with one unit. For this study, the primary unit is taken to be one-fourth of the classical lobule. The relevant anatomical values are listed in Table 1.

Elimination kinetics
Here drug uptake and elimination (i.e. conversion to metabolized product) is viewed as a single-step saturable process following Michaelis-Menten kinetics [20], such that Here C(t) is the local drug concentration and dC(t)/dt is the drug metabolization rate. (In what follows, C(t) is expressed as ρx i with ρ the fluid density and x i mole fractions of i-th species). Note we are explicitly modeling the drug transport to an individual hepatocyte surface via our lobule lattice model and assume an effective one-step reaction transformation beyond that point. We recognize that drug incorporation and elimination is still a multi-step process even once the drug reaches the cell surface however. It is hoped that these approximations ignore processes occurring on a shorter time scale than the experimental resolution. Nonetheless, we explore possible complications via simple sensitivities to the choice of reaction time constant.

Model and method
The primary unit of the liver was approximated by a symmetry element of a 51 × 51 square lattice such that four units make up one lobule. The architecture of the lattice consists of hepatocyte grid cells interlaced by a network of narrower sinusoidal grid cells ( Figure 1a). This "intermediate realism" model represents a practical compromise with simulation run times of seconds compared with precise, high resolution, high realism models which we will explore in future work. These latter models can be expected to have equivalent run times of hours and before proceeding to such simulations, it was deemed prudent to explore a wide range of process variables including both lattice effects as well as metabolic effects with our simpler representation in order to get an overall assessment of liver lobule performance and drug metabolism in healthy and diseased livers. Additionally, realistic, high resolution models by their very nature can contain conflicting mechanisms whereas simpler models allow their separation for analysis. In our model, the vascular supply tree (hereafter termed "injector") is located at top left-corner (grid cell 1,1) while the vascular collecting tree (hereafter termed "producer") is located at bottom right-corner (grid cell 51,51) of this model. A complete lobule would be formed by placing three more equivalent lattices around the central injector but because of symmetry for this ideal lattice model, all 4 sections would produce equivalent profiles and hence be computationally redundant. The diameter of the sinusoid grid cells was taken to be 0.0006 cm, and the diameter of the hepatocyte grid cells was taken to be 0.0024 cm. (Again because of symmetry, the bounding sinusoids are taken as ½ size, 0.0003 cm, as the remaining portion of these sinusoids are included in the adjacent sections). The length of the lattice is thus 0.0750 cm per side. Doubling this value gives a lobule diameter of 0.150 cm, which is consistent with values listed in Table 1.
Convective molar flux is modelled according to Darcy's Law [21]: where J c ik is the i-th component of fluid flux in k-direction, ρ and μ are fluid molar density and viscosity, v k , K k and ∇ k p are the Darcy velocity, permeability, and pressure gradient in direction k, respectively. The blood viscosity μ is taken to be 3.5 mPa-s (3.5 centipoise). Blood molar density is assumed that of water, ρ at 55.4 mmoles/cm 3 . It is emphasized that in this paper, following Darcy's Law and the conventions of flow in porous media, permeability K is defined as a measure of the transmissibility of a grid cell to the flow of a fluid, and is expressed in units of area (e.g. cm 2 ).
Each sinusoid grid cell represents a tubular vessel of diameter 2a. Taking the ratio of the volume of the vessel to the volume of the grid cell yields a porosity of For a cylindrical tube, the permeability can be calculated from the tube radius [14] K sin ¼ Here and in what follows, it is noted that porosity is dimensionless while permeability has units of a characteristic length squared.
Each parenchymal grid cell represents a cellular (hepatocyte) component and an extracellular (space of Disse) component. A ratio of 0.75 to 0.25 was chosen for their respective contributions to the volume, and the porosity of the parenchymal sites was therefore 0.25.
The corresponding permeability of the tissue grid cells is estimated from a Carmen-Kozeny model [21] of flow around a spherical object (the hepatocyte) of radius 12 μm. This states the permeability is proportional to the object's diameter D p and the porosity ϕ tis as follows: An ideal result where D p /L = 1 or R/L = 0.5 with L as grid size and R as particle radius gives Then assuming as we do, L = D p = 24 μm, the Carmen-Kozeny formula yields The Carman-Kozeny expression basically states physically that the order of magnitude of the permeability scales with the particle size squared.
The above analysis is an ideal result as the ECM around the cell particle will further reduce porosity. Thus we could study a range of tissue effective porosities, leading to a range of effective tissue permeabilities. Table 2 summarizes a range of possible values. In this paper, we will utilize the "base case" value for tissue porosity and permeability, while the effects of more extreme choices will be examined in a second paper.
The convective driving force originates from an input site corresponding to a terminal portal venule at one corner of the lattice and an output site corresponding to a terminal hepatic venule at the opposite corner. For simplicity, the hepatic artery blood supply, which is lower in volume and pulsatile in nature, is omitted for the current simulations. The pressure value at the inlet and the outlet are taken to be P in = 103 kPa and P out = 101.8 kPa, respectively. After subtracting the atmospheric pressure, these values are consistent with experimental values quoted by Rappaport [22], who found that the terminal portal venule pressure was in the range 0.59 kPa to 2.45 kPa and that the terminal hepatic venule pressure was 0.49 kPa. As we shall demonstrate, this applied pressure differential results in a convective flow level that is determined primarily by the effective permeability of the lobule. Thus various liver damage scenarios can be expected to affect this flow. This aspect of the modelling is of practical importance and will be explored in more detail in a separate publication.
For multicomponent flow, the model tracks the compositions (molar or mass fractions) of all components in the fluid. In this paper, we focus on the injected drug paclitaxol (PAC) and the Phase I transformed metabolite 6-hydroxypaclitaxel (PAC-OH). In addition to convective transport a diffusive flux contribution of is considered with J d ik being the molar diffusive flux and D ik the diffusion constant of species i in direction k. The estimated diffusion constant in all directions used here is based on a molecular weight rescaling of glucose diffusion. Here, glucose diffusion coefficient in water is taken as a basic reference value for comparison, (E. L. Cussler, "Diffusion: Mass Transfer in Fluid Systems"), With a cubic root of the molecular weight ratio of glucose to PAC used as conversion factor (Factor = (180/854) 0.33 = (1/4.74) 0.33 = 1/1.68), an estimated effective diffusion constant for PAC is Tissue effective diffusion value should be less; here we employ an order of magnitude reduction in the value of D D pac;k ¼ 7:1 Â 10 −11 1:68 Effective diffusion constants for PAC-OH are assumed identical to PAC values. These values are converted to the simulation units of cm 2 /min and also summarized in Table 2.
The drug paclitaxel was used as a reactive tracer, and its Phase I metabolism was modeled using the general formula of one paclitaxel (PAC) molecule being transformed into the metabolite (PAC-OH) by the cytochrome P450 (CYP) isozyme CYP2C8 [23].
The enzyme only exists in grid cells containing hepatocytes so all reaction is localized in these sites. In this paper, saturable Michaelis-Menten kinetics are assumed, defined by a maximum rate v max (in units of molar fraction/min) and half saturation value K m (in units of molar fraction). This reaction proceeds in a linear manner at a rate characterized by k = v max /K m (in units of min -1 ) when injection concentrations are much below the half saturation value.
Reaction parameter values are based on the work of Vaclavikova et al. [24] who measured directly PAC conversion to PAC-OH kinetics without any tissue distribution issues and uptake by the cell itself, as they use microsomes directly as the source of CYP. As such, any bottlenecks associated with drug uptake should imply that reaction rates would be slower than those based on parameters values given by Vaclavikova et al. [24]. We are modeling tissue distribution effects separately based on our lobule model. Table 3 summarizes the reaction parameters.
The simulations were performed using the STARS advanced process simulator designed by the Computer Modelling Group (CMG) Ltd. in Calgary, Alberta, to model the flow and reactions of multiphase, multicomponent fluids through porous media [25][26][27]. Additionally, STARS has earlier been used to model reactive flow processes in cortical bone [28][29][30][31] as well as through the intervertebral disk [32].

Non-reactive flow characteristics
As discussed above, flow is induced on the regular lattice of Figure 1b by applying a pressure difference across the inlet and outlet points. With the chosen lobule flow parameters for porosity, permeability, and blood viscosity, this translates to a steady flow rate 2.1 cm 3 /min as illustrated in Figure 2a. A short timescale of about 2.0 × 10 -5 min needed to establish this pressure gradient is also illustrated in Figure 2b by expanding the time axis. Figure 3 demonstrates the steady state velocity profile throughout the lattice, illustrating both the diverging/converging nature of the flow near the inlet and outlet points at the top-left and bottom-right portion of the grid (ie injector and producer, respectively), as well as the orders of magnitude difference of the flows in the sinusoids and tissues, respectively. (This plot uses a logarithmic colour scale axis).
When blood with a relative composition of 1 micro-gram paclitaxel (1.8 × 10 -8 mole fraction) is infused into the lattice assuming nonreactive hepatocytes, the time Linear rate v max /K m 6.0 × 10 -3 min -1 6.0 × 10 -3 min -1 *Their Table 4 quotes v max = 61 picomol/mg_protein/min. Using their microsomal protein concentration of 1 mg/ml, these numbers convert to v max = 0.06 μM/min. required to traverse the lattice is approximately 1 min as demonstrated in Figure 4. This production profile is convective flow dominated as the addition of diffusion minimally alters the production profile. The evolution of the paclitaxel concentration on the lattice followed a spatially homogeneous progression ( Figure 5), which shows the increasing levels of injected drug after 0.01 min and 0.14 min. By 0.14 min, paclitaxel is being seen at the outlet of the lobule. After 0.5 min, paclitaxel completely covers the lattice (not shown).
If the diffusive flow contribution is removed, however, the paclitaxel profiles on the lattice are significantly different. As is also illustrated in Figure 5, again at 0.01 min and 0.14 min, a distinct two-scale behavior is noted, whereby the sinusoids are first infused with the drug, and only at later times do the drug levels in the tissue approach injected concentration levels. This behavior reflects the convective levels of flow in the sinusoids and tissues noted earlier (Figure 3). A further comparison of Figure 5 cases reveals that the sinusoid drug concentration levels in the two cases are similar, however, explaining   the similar drug production characteristics note in Figure 4, as paclitaxel is produced directly from the sinusoids.

Base case reactive flows
The effects of paclitaxel drug metabolism by hepatocytes are next considered. Here the base case reaction parameters of Table 3 are employed, and the same injected paclitaxel concentration (1.8 × 10 -8 mole fraction) is considered. With the employed reaction half saturation constant value of 1.8 × 10 -7 mole fraction, this injection level implies the Michaelis-Menten model reduces to an almost linear reaction scheme. Figure 6 illustrates injected drug and produced drug and metabolite production for this case. Again it is emphasized that both PAC and PAC-OH have assumed equal diffusive flow contributions, as these are components of very similar size. Essentially at this reaction rate, all injected paclitaxel is converted to metabolite by the lobule hepatocytes. The production profile of PAC-OH here is identical to the production profile of PAC in the non-reaction case, as shown in Figure 4. Figure 7 shows the PAC and PAC-OH profiles across the lobule lattice at 0.01 min, 0.14 min, and 0.50 min, respectively. The PAC concentrations in the sinusoids and the PAC-OH concentrations in the tissue are equivalent to the PAC concentrations in both sinusoids and tissue for the non-reacted case ( Figure 5). Figure 7 also shows most clearly there is an inlet distance over which the reaction conversion time is not fast enough to convert the injected paclitaxel. Figure 8 illustrates injected drug and produced drug and metabolite production for the same case except that diffusive transport has been removed. In contrast to Figure 6 with diffusion, there is now only a limited amount of conversion of PAC to PAC-OH even at long times. The PAC and PAC-OH profiles at various times (0.01 min, 0.14 min, and 0.50 min) as shown in Figure 9, confirm this behavior where it is shown that the PAC concentration in the sinusoids propagates throughout the lobule, while the PAC-OH concentrations in the tissue increase less rapidly and up to a lower level.

Reactive flow sensitivities
In this section, the consequences of the chosen reaction parameters are illustrated. Figure 10 shows production behavior with a 100-fold reduction in maximum reaction  rate and with diffusion effects included. The small level of produced PAC indicates that reaction rates must be reduced to about this level before any significant change in drug production behavior can be expected. Drug distribution in the lobule for this case is shown in Figure 11 for the times 0.01 min, 0.14 min, and 0.50 min. This figure should be contrasted with Figure 9. Here the early time results and upstream results for PAC distributions at longer times are quite different, reflecting the reduced reaction rate. However, the later time and downstream results for PAC-OH distribution resemble  quite closely the faster reaction limit. Here the propagation time across the lobule gives enough time to compensate for changes in reaction rate. In summary, reaction rates larger than the base case or even 10-fold reduction from base case can be expected to produce very similar drug production behavior and differences only in the inlet region of the lobule are to be envisioned.
In contrast, once a critical time-scale is crossed, much more significant changes in drug distribution behavior can be expected, both internally throughout the lobule and in terms of produced profiles. Figure 12 shows drug metabolite production behavior with a 1000-fold reduced metabolic rate, and including diffusive mixing. Here almost equal levels of PAC and PAC-OH are seen exiting the lobule. At 0.01 min almost no PAC-OH is converted in the lobule at this rate (see Figure 13), while at later times (0.14 min and 0.50 min), converted PAC-OH starts to be seen at the outlet regions at levels similar to PAC. Essentially, the inlet region behavior occurring at faster reaction rates now covers the whole lobule region.
Finally, sensitivities to injected PAC concentrations were explored, utilizing injection concentrations of 1.8 × 10 -7 and 1.8 × 10 -6 mole fractions (i.e. clearly above that of the base case 1.8 × 10 -8 ). In these runs, the base case reaction parameters were maintained. In particular, the half saturation value of 1.8 × 10 -7 was employed, indicating that the linear, intermediate, and saturation levels of the Michaelis-Menten expression were being probed with the three injected concentration levels. As illustrated in Figure 14 for the runs without diffusion, the production profiles of PAC-OH remained unchanged for each case, as long as the production maxima were rescaled to the corresponding injection concentrations. Apparently, with fast reaction rates, the Michaelis-Menten form had little impact on production behavior.

Analysis
The results we have presented can be rationalized by a comparison of process timescales. Calculation of breakthrough times can be based on two concepts: either only the sinusoids are accessible or the whole lobule (tissue + sinusoid) is accessible to injected species. The sinusoid pore volume in our element is 4.77 × 10 -7 cm 3 while the complete lobule element volume is 7.344 × 10 -7 cm 3 . Because the steady state flow rate in our model is 2.44 × 10 -6 cm 3 /min (see Figure 2), this means the breakthrough time is  These should be viewed as two limiting vertical lines on the time history plots as two ideal limits without any diffusion or mixing effects (physical or numerical). It is noted for example that our time history plot of PAC production with no reaction and with or without physical diffusion (see Figure 4) has a produced concentration of 0.9 × 10 -8 (i.e. half of the injected 1.8 × 10 -8 concentration) at 0.19 min, about what is expected. The main point here is that most of the production behavior differences for our various cases should lie between these two ideal "half-value" limits.
The next timescale is governed by a "pressure diffusion" coefficient Here C p eff is an effective compressibility accounting for both fluid and tissue structure effects. Fluid (water compressibility) is of the order of 5 × 10 -7 kPa -1 . For liver (soft tissue) structural compressibility, we have chosen 1.8 × 10 -5 kPa -1 . Using these choices and the base case parameters of this study from Table 2, we obtain D pres ≅1:5 Â 10 4 cm 2 min ð13Þ This parameter essentially describes the time taken for pressure to come to a steady state distribution as follows. Utilizing a characteristic distance d = 0.15 cm (the lobule element size for the pressure calculation), this time is then  pulsatile flow effects caused by hepatic artery inflow, this timescale would be much more important to the general process description and a more precise definition of compressibility might be warranted. For the present, these numbers just indicate that steady state pressure is achieved more quickly than other process effects.
The third timescale is determined by particle diffusion. Here we have chosen diffusion constants based on paclitaxel size and simple estimates of tortuosity. The parameter choices used here are D sin = 2.5 × 10 -4 cm 2 /min and D tis = 2.5 × 10 -5 cm 2 /min. As seen from our simulations, for well-perfused tissue convective effects operate globally over the whole domain while diffusion smooths concentration profiles locally. Thus with a choice of characteristic distance d = 10 μm = 1 × 10 -3 cm, the times required for particles to diffuse are Our diffusion values should be viewed as highly optimistic. In particular, pactlitaxol is normally not molecularly dissolved, but rather it is some type of micellar complex with Cremophor EL surfactant, so the effective diffusion coefficient for this complex is probably one or more orders of magnitude smaller than what has been estimated. Thus the limits of diffusion and non-diffusion cases are meaningful extremes of what might be expected, for small molecules and large nanoparticles, respectively.
The final timescale is reaction rate. The base limiting reaction rate 6 × 10 -3 min -1 , also utilizing a model tissue volume to bulk volume scale factor of SF~0.76 × 10 -6 , converts to a reaction time (reaction half-life) of This is essentially seen as more rapid than or comparable to the other timescales considered (see the early time PAC-OH tissue concentration level appearance in Figures 7d  and 9d). Reducing this basic reaction rate by a factor of 100 or 1000 causes the reaction process to be more similar to the other timescales and different production profiles of PAC and PAC-OH result, as has been shown. Table 4 summarizes the relevant assumed timescales.

Conclusion
In pharmacokinetics, lattice models are introduced to address the non-heterogeneity of the organs on the drug distribution that has a significant impact on drug propagation throughout the body as shown by analysis of clinical data.
Here we utilize the interpretation of the liver as an ensemble of islands of metabolic activity and focus on the liver lobule itself. In contrast to most of earlier studies that assumed that drug molecules only randomly propagate (ie diffuse) through the system, we have emphasized the dominant role of convection in well-vascularized tissue with a given structure. We have utilized an idealized representation to analyze the factors affecting drug propagation and metabolism. The lobule is divided into hepatocyte cells that are interlaced with narrower sinusoidal grid cells. These cells are connected by constant permeability throughout the entire system. The drug molecules convectively flow through the sinusoidal along with blood (water here) and diffuse to the hepatocyte where metabolisms are taking place. A sensitivity analysis of convective, diffusive and reaction parameters are performed and estimates of their contributions are presented. We have compared the drug concentration levels observed exiting the lobule with their predicted detailed distribution inside the lobule, assuming that most often the former is accessible information while the latter is not. As such, we establish how traditional pharmacokinetic analysis might be reflective of the spatial distribution of the drug in the lobule, and situations when this might not hold. Interestingly, our network models including dispersive effects often correspond to the "well-stirred" compartment models, such that relatively uniform steady-state concentration levels occur throughout the lobule (if one ignores the smaller inlet mixing zone). Conversely, simulations on our network models without explicit dispersive mixing often correspond to modified "parallel tube" models, such that observed concentration profiles change along the length of the tubes (ie sinusoids). Here our modified tube network structure allows cross sinusoids as well. These comments reflect Figures 3, 5, and 9. This is the first paper of series of papers on physiologically-based lattice models for liver. In this paper, we consider an idealized lobule lattice in order to understand the basic functionality of the unit and underlying mechanisms through simulations and also to set a basis for future studies. In following papers, we will expand the analysis to include drug propagation & metabolic sensitivities associated with variations in lobule structure, which could reflect extents of liver damage, and how our modelling approach might be used to generate flows on realistic, reconstructed images of lobule structure.