A physiologically-based flow network model for hepatic drug elimination I: regular lattice lobule model
© Rezania et al.; licensee BioMed Central Ltd. 2013
Received: 28 May 2013
Accepted: 22 August 2013
Published: 5 September 2013
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.
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 . 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.
Multiple computational models have been proposed: single and multi-compartment models, distribution-based models, agent-based models, interconnected parallel tubes, etc. See  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  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–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 . 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.
The liver and drug kinetics
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 . 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 . Finally, the blood is collected and removed by the hepatic venules.
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.
where Jc 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/cm3. 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. cm2).
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 Carman-Kozeny expression basically states physically that the order of magnitude of the permeability scales with the particle size squared.
Lobule regular lattice flow parameters
Characteristic (SI) unit
Sinusoid Grid Cell Size
Sinusoid Porosity: ϕsin
Sinusoid Permeability: Ksin
Sinusoid Effective Diffusion: Dsin
4.2 × 10-10 m2/sec
2.5 × 10-4 cm2/min
Tissue Grid Cell Size: 2a
Tissue Porosity (Ideal): ϕtis
Tissue Permeability (Ideal): Ktis
Tissue Porosity (Base): ϕtis
Tissue Permeability (Base): Ktis
7.35 × 10-2 μm2
7.45 × 10-2 Darcy
Tissue Porosity (ECM): ϕtis
Tissue Permeability (ECM): Ktis
6.883 × 10-3 μm2
6.97 × 10-3 Darcy
Tissue Effective Diffusion: Dtis
4.2 × 10-11 m2/sec
2.5 × 10-5 cm2/min
Blood Viscosity: μ
3.5 × 10-3 Pa-sec
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 Pin = 103 kPa and Pout = 101.8 kPa, respectively. After subtracting the atmospheric pressure, these values are consistent with experimental values quoted by Rappaport , 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.
Effective diffusion constants for PAC-OH are assumed identical to PAC values. These values are converted to the simulation units of cm2/min and also summarized in Table 2.
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 vmax (in units of molar fraction/min) and half saturation value Km (in units of molar fraction). This reaction proceeds in a linear manner at a rate characterized by k = vmax/Km (in units of min-1) when injection concentrations are much below the half saturation value.
Paclitaxel kinetic elimination Michaelis-Menten parameters (converted* from Vaclavikova et al., their Table 4 )
Characteristic (SI) unit
Maximum rate vmax
1.08 × 10-9 molefrac/min
Half saturation constant Km
1.8 × 10-7molfrac
Linear rate vmax/Km
6.0 × 10-3 min-1
6.0 × 10-3 min-1
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–27]. Additionally, STARS has earlier been used to model reactive flow processes in cortical bone [28–31] as well as through the intervertebral disk .
Non-reactive flow characteristics
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.
Reactive flow sensitivities
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 cm3 while the complete lobule element volume is 7.344 × 10-7 cm3. Because the steady state flow rate in our model is 2.44 × 10-6 cm3/min (see Figure 2), this means the breakthrough time is 4.77/24.4 = 0.20 min for just sinusoid accessibility and 7.34/24.4 = 0.30 min for the whole lobule-sampled space. 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.
Figure 2 also illustrates this characteristic time. This timescale is essentially a function of fluid properties and lobule structure (through ϕ and K). If we were to consider 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.
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.
Assumed lobule process time constants
Convective transit time (sinusoid network only)
Pressure relaxation time Constant (in sinusoids)
1.5 × 10-5 min
Diffusion relaxation time constant (CYP/CYP-OH in sinusoids)
4.0 × 10-3 min
Diffusion relaxation time constant (CYP/CYP-OH in tissue)
4.0 × 10-2 min
Base case metabolic uptake/elimination time constant
1.0 × 10-4 min
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.
J.A.T. acknowledges funding support for this project from NSERC. The Allard Foundation and the Alberta Advanced Education and Technology.
- Ierapetritou MG, Georgopoulos PG, Roth CM, Androulakis IP: Tissue-level modeling of xenobiotic metabolism in liver: an emerging tool for enabling clinical translational research. Clin Transl Sci. 2009, 2: 228-237. 10.1111/j.1752-8062.2009.00092.x.PubMed CentralView ArticleGoogle Scholar
- Jacquez JA: Compartmental analysis in biology and medicine. 1996, Ann Arbor, MI: BioMedware, 3Google Scholar
- Bass L, Keiding S, Winkler K, Tygstrup N: Enzymatic elimination of substrates flowing through the intact liver. J Theor Biol. 1976, 61: 393-409. 10.1016/0022-5193(76)90026-6.View ArticlePubMedGoogle Scholar
- Bass L, Robinson P, Bracken A: Hepatic elimination of flowing substrates: the distributed model. J Theor Biol. 1978, 72: 161-184. 10.1016/0022-5193(78)90023-1.View ArticlePubMedGoogle Scholar
- Pang KS, Rowland M: Hepatic clearance of drugs I: theoretical considerations of a well-stirred model and a parallel tube model. J Pharmacokin Biopharm. 1977, 5: 625-653. 10.1007/BF01059688.View ArticleGoogle Scholar
- Pang KS, Rowland M: Clearance of drugs I: experimental evidence. J Pharmacokin Biopharm. 1977, 5: 655-680. 10.1007/BF01059689.View ArticleGoogle Scholar
- Roberts MS, Rowland M: A dispersion model of hepatic elimination 1. formulation. J Pharmacokin Biopharm. 1986, 14: 227-260. 10.1007/BF01106706.View ArticleGoogle Scholar
- Roberts MS, Rowland M: A dispersion model of hepatic elimination 2. steady-state considerations. J Pharmacokin Biopharm. 1986, 14: 261-288. 10.1007/BF01106707.View ArticleGoogle Scholar
- Roberts MS, Rowland M: A dispersion model of hepatic elimination 3. application to metabolic formation and elimination kinetics. J Pharmacokin Biopharm. 1986, 14: 289-308. 10.1007/BF01106708.View ArticleGoogle Scholar
- Bass L, Roberts MS, Robinson PJ: On the relation between extended forms of the sinusoidal perfusion and of the convection–dispersion models of hepatic elimination. J Theor Biol. 1987, 126: 457-482. 10.1016/S0022-5193(87)80152-2.View ArticlePubMedGoogle Scholar
- Roberts MS, Donaldson JD, Rowland M: Availability predictions by hepatic elimination models for Michaelis-Menten kinetics. J Pharmacokin Biopharm. 1989, 17: 687-719. 10.1007/BF01062125.View ArticleGoogle Scholar
- The Liver: Biology and Pathology. Edited by: Arias IM. 2001, Philadelphia: Lippincott Williams and Wilkins, 4Google Scholar
- Goresky CA: Uptake in the liver: the nature of the process. Int Rev Physiol. 1980, 21: 65-101.PubMedGoogle Scholar
- Saxena R, Theise ND, Crawford JM: Microanatomy of the human liver – exploring the hidden interfaces. Hepatology. 1999, 30: 1339-1346. 10.1002/hep.510300607.View ArticlePubMedGoogle Scholar
- Bhunchet E, Wake K: The portal lobule in rat liver fibrosis: a re-evaluation of the liver unit. Hepatology. 1998, 27 (2): 481-487. 10.1002/hep.510270223.View ArticlePubMedGoogle Scholar
- Matsumoto T, Kawakami M: The unit concept of hepatic parenchyma – a re-examination based on angioarchitectural. Acta Pathol Jpn. 1982, 32 (suppl2): 285-314.PubMedGoogle Scholar
- Teutsch HF: The modular microarchitecture of human liver. Hepatology. 2005, 42: 317-325. 10.1002/hep.20764.View ArticlePubMedGoogle Scholar
- Teutsch HF, Schuerfeld D, Groezinger E: Three-dimensional reconstruction of parenchymal units in the liver of the rat. Hepatology. 1999, 29: 494-505. 10.1002/hep.510290243.View ArticlePubMedGoogle Scholar
- Ruijter JM, Markman J, Hagoort MM, Moorman AF, Lamers WH: Relative distance: the key to the shape of hepatic building blocks. Image Anal Sterol. 2000, 19 (1): 19-24. 10.5566/ias.v19.p19-24.View ArticleGoogle Scholar
- Marsh RE, Tuszynski JA, Sawyer M, Voss K: Emergence of power laws in the pharmacokinetics of paclitaxol due to competing saturable processes. J Pharm Pharmaceut Sci. 2008, 11 (3): 77-96.Google Scholar
- Dullien FAL: Porous Media: Fluid Transport and Pore Structure. 1992, San Diego, CA, USA: Academic Press, 2Google Scholar
- Rappaport AM: Hepatic blood flow: morphologic aspects and physiologic regulation. Int Rev Physiol. 1980, 21: 1-63.PubMedGoogle Scholar
- Monsarrat B, Chatelut E, Royer I, Alvinerie P, Dubois J, Dezeus A, Roche H, Cros S, Wright W, Canal P: Modification of paclitaxol metabolism in cancer patient by induction of cytochrome P450 3A4. Drug Met Disp. 1998, 26: 229-233.Google Scholar
- Vaclavikova R, Soucek P, Svobodova L, Anzenbacher P, Simek P, Guengerich F, Gut I: Different in vitro metabolism of paclitaxel and docetaxoel in humans, rats, pigs, and minipigs. Drug Metab Dispos. 2004, 32 (6): 666-674. 10.1124/dmd.32.6.666.View ArticlePubMedGoogle Scholar
- CMG Ltd: STARS User’s Guide: Advanced Process and Thermal Reservoir Simulator. 2011, Calgary, AB: Computer Modelling Group LtdGoogle Scholar
- Oballa V, Coombe D, Buchanan W: Factors affecting the thermal response of naturally fractured reservoirs. J Can Pet Tech. 1993, 32 (8): 31-37.View ArticleGoogle Scholar
- Darche G, Grabenstetter JE, Sammon PH: The use of parallel processing with dynamic gridding. SPE Reservoir Simul Symp. 2005, paper 93023-doi:10.2118/93023-MS. Houston, TX, USAGoogle Scholar
- Goulet G, Hamilton N, Cooper D, Coombe D, Tran D, Martinuzzi R, Zernicke R: Influence of vascular porosity on fluid flow and nutrient transport in loaded cortical bone. J Biomech. 2008, 41 (10): 2169-2175. 10.1016/j.jbiomech.2008.04.022.View ArticlePubMedGoogle Scholar
- Goulet GC, Cooper DML, Coombe D, Zernicke RF: Influence of cortical canal architecture on lacunocanalicular pore pressure and fluid flow. Comput Methods Biomech Biomed Eng. 2008, 11 (4): 379-387. 10.1080/10255840701814105.View ArticleGoogle Scholar
- Goulet GC, Cooper DML, Coombe D, Zernicke RF: Poroelastic evaluation of fluid movement through the Lacunocanicular system. Annals of Biomed Eng. 2009, 37 (7): 1390-1402. 10.1007/s10439-009-9706-1.View ArticleGoogle Scholar
- Goulet GC, Cooper DML, Coombe D, Zernicke RF: Validation and application of iterative coupling to poroelastic problems in bone fluid flow. Bull Appl Mech. 2009, 5 (1): 6-16.Google Scholar
- Louman-Gardiner KM, Coombe D, Hunter CJ: Computational models simulating notochordal cell extinction during early aging of an intervertebral disk. Comput Methods Biomech Biomed Eng. 2011, 14 (12): 1071-1077. 10.1080/10255842.2010.508037.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.