Skip to main content

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.


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.

Figure 1
figure 1

Flow network structure. (a) Schematic diagram of a cross section of hepatic parenchyma consisting hexagonal lobules, portal and hepatic veins. The lobule contains liver cells (sinusoids and hepatocytes). The segment represents a typical area studied in this paper. (b) Homogeneous lattice (segmented area) with high porosity bands (in red) representing sinusoids and lower porosity regions (in blue) representing tissue containing hepatocytes. (c) 2D hexagonal lobule lattice with sinusoids generated via a diffusion limited aggregation algorithm. (d) 3D hexagonal lobule lattice.

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 [711] 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.

The liver and drug kinetics

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.

Table 1 Anatomical parameter values for the liver[12, 13]

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 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).

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]


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 Dp and the porosity ϕtis as follows:


An ideal result where Dp/L = 1 or R/L = 0.5 with L as grid size and R as particle radius gives


Then assuming as we do, L = Dp = 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.

Table 2 Lobule regular lattice flow parameters

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 [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 Jd 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


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 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 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.

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.

Table 3 Paclitaxel kinetic elimination Michaelis-Menten parameters (converted* from Vaclavikova et al., their Table 4 )

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 [2527]. Additionally, STARS has earlier been used to model reactive flow processes in cortical bone [2831] 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 cm3/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).

Figure 2
figure 2

Injected flow (PAC). (a) Steady state flow across the lobule, (b) Short time flow transients.

Figure 3
figure 3

Steady state velocity profile across the lobule. Color bar is in cm/min.

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 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.

Figure 4
figure 4

Non-reactive PAC drug propagation across the lobule, with and without diffusion effects.

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).

Figure 5
figure 5

Non-reactive PAC profiles across the lobule wit and without diffusion. (a) PAC at 0.01 min, with diffusion effect, (b) PAC at 0.14 min, with diffusion effect, (c) PAC at 0.01 min, no diffusion effect, (d) PAC at 0.14 min, no diffusion effect. Color bar is in molfrac.

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 6
figure 6

Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and base case metabolism.

Figure 7
figure 7

Reactive PAC and PAC-OH profiles across the lobule with diffusion effects and base case metabolism. (a) PAC at 0.01 min, (b) PAC at 0.14 min, (c) PAC at 0.50 min, (d) PAC-OH at 0.01 min, (e) PAC-OH at 0.14 min, (f) PAC-OH at 0.50 min. Color bar is in molfrac.

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.

Figure 8
figure 8

Reactive PAC and PAC-OH drug propagation across the lobule, without diffusion effects and base case metabolism.

Figure 9
figure 9

Reactive PAC and PAC-OH profiles across the lobule without diffusion effects and base case metabolism. (a) PAC at 0.01 min, (b) PAC at 0.14 min, (c) PAC at 0.50 min, (d) PAC-OH at 0.01 min, (e) PAC-OH at 0.14 min, (f) PAC-OH at 0.50 min. Color bar is in molfrac.

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.

Figure 10
figure 10

Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and 100-fold reduced metabolism.

Figure 11
figure 11

Reactive PAC and PAC-OH profiles across the lobule with diffusion effects and 100-fold reduced metabolism. (a) PAC at 0.01 min, (b) PAC at 0.14 min, (c) PAC at 0.50 min, (d) PAC-OH at 0.01 min, (e) PAC-OH at 0.14 min, (f) PAC-OH at 0.50 min. Color bar is in molfrac.

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.

Figure 12
figure 12

Reactive PAC and PAC-OH drug propagation across the lobule, with diffusion effects and 1000-fold reduced metabolism.

Figure 13
figure 13

Reactive PAC and PAC-OH profiles across the lobule with diffusion effects and 1000-fold reduced metabolism. (a) PAC at 0.01 min, (b) PAC at 0.14 min, (c) PAC at 0.50 min, (d) PAC-OH at 0.01 min, (e) PAC-OH at 0.14 min, (f) PAC-OH at 0.50 min. Color bar is in molefrac.

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.

Figure 14
figure 14

PAC-OH metabolite production levels from various injected PAC concentrations.


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.

The next timescale is governed by a “pressure diffusion” coefficient


Here Cpeff 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


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


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.

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 Dsin = 2.5 × 10-4 cm2/min and Dtis = 2.5 × 10-5 cm2/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.

Table 4 Assumed lobule process time constants


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.


  1. 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.

    Article  PubMed Central  CAS  Google Scholar 

  2. Jacquez JA: Compartmental analysis in biology and medicine. 1996, Ann Arbor, MI: BioMedware, 3

    Google Scholar 

  3. 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.

    Article  CAS  PubMed  Google Scholar 

  4. 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.

    Article  CAS  PubMed  Google Scholar 

  5. 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.

    Article  CAS  Google Scholar 

  6. Pang KS, Rowland M: Clearance of drugs I: experimental evidence. J Pharmacokin Biopharm. 1977, 5: 655-680. 10.1007/BF01059689.

    Article  CAS  Google Scholar 

  7. Roberts MS, Rowland M: A dispersion model of hepatic elimination 1. formulation. J Pharmacokin Biopharm. 1986, 14: 227-260. 10.1007/BF01106706.

    Article  CAS  Google Scholar 

  8. Roberts MS, Rowland M: A dispersion model of hepatic elimination 2. steady-state considerations. J Pharmacokin Biopharm. 1986, 14: 261-288. 10.1007/BF01106707.

    Article  CAS  Google Scholar 

  9. 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.

    Article  CAS  Google Scholar 

  10. 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.

    Article  CAS  PubMed  Google Scholar 

  11. 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.

    Article  CAS  Google Scholar 

  12. The Liver: Biology and Pathology. Edited by: Arias IM. 2001, Philadelphia: Lippincott Williams and Wilkins, 4

    Google Scholar 

  13. Goresky CA: Uptake in the liver: the nature of the process. Int Rev Physiol. 1980, 21: 65-101.

    CAS  PubMed  Google Scholar 

  14. Saxena R, Theise ND, Crawford JM: Microanatomy of the human liver – exploring the hidden interfaces. Hepatology. 1999, 30: 1339-1346. 10.1002/hep.510300607.

    Article  CAS  PubMed  Google Scholar 

  15. 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.

    Article  CAS  PubMed  Google Scholar 

  16. Matsumoto T, Kawakami M: The unit concept of hepatic parenchyma – a re-examination based on angioarchitectural. Acta Pathol Jpn. 1982, 32 (suppl2): 285-314.

    PubMed  Google Scholar 

  17. Teutsch HF: The modular microarchitecture of human liver. Hepatology. 2005, 42: 317-325. 10.1002/hep.20764.

    Article  PubMed  Google Scholar 

  18. 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.

    Article  CAS  PubMed  Google Scholar 

  19. 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.

    Article  Google Scholar 

  20. 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.

    CAS  Google Scholar 

  21. Dullien FAL: Porous Media: Fluid Transport and Pore Structure. 1992, San Diego, CA, USA: Academic Press, 2

    Google Scholar 

  22. Rappaport AM: Hepatic blood flow: morphologic aspects and physiologic regulation. Int Rev Physiol. 1980, 21: 1-63.

    CAS  PubMed  Google Scholar 

  23. 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.

    CAS  Google Scholar 

  24. 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.

    Article  CAS  PubMed  Google Scholar 

  25. CMG Ltd: STARS User’s Guide: Advanced Process and Thermal Reservoir Simulator. 2011, Calgary, AB: Computer Modelling Group Ltd

    Google Scholar 

  26. Oballa V, Coombe D, Buchanan W: Factors affecting the thermal response of naturally fractured reservoirs. J Can Pet Tech. 1993, 32 (8): 31-37.

    Article  CAS  Google Scholar 

  27. 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, USA

    Google Scholar 

  28. 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.

    Article  PubMed  Google Scholar 

  29. 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.

    Article  CAS  Google Scholar 

  30. 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.

    Article  Google Scholar 

  31. 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 

  32. 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.

    Article  CAS  Google Scholar 

Download references


J.A.T. acknowledges funding support for this project from NSERC. The Allard Foundation and the Alberta Advanced Education and Technology.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Jack Tuszynski.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

VR: Primary investigator who conducted the majority of the simulations and model development. RM: Initial investigator who conducted the preliminary simulations and model development. DC: Investigator providing technical support for the simulations and model development. JT: Lead investigator who supervised the work content and model development. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Authors’ original file for figure 4

Authors’ original file for figure 5

Authors’ original file for figure 6

Authors’ original file for figure 7

Authors’ original file for figure 8

Authors’ original file for figure 9

Authors’ original file for figure 10

Authors’ original file for figure 11

Authors’ original file for figure 12

Authors’ original file for figure 13

Authors’ original file for figure 14

Authors’ original file for figure 15

Authors’ original file for figure 16

Authors’ original file for figure 17

Authors’ original file for figure 18

Authors’ original file for figure 19

Authors’ original file for figure 20

Authors’ original file for figure 21

Authors’ original file for figure 22

Authors’ original file for figure 23

Authors’ original file for figure 24

Authors’ original file for figure 25

Authors’ original file for figure 26

Authors’ original file for figure 27

Authors’ original file for figure 28

Authors’ original file for figure 29

Authors’ original file for figure 30

Authors’ original file for figure 31

Authors’ original file for figure 32

Authors’ original file for figure 33

Authors’ original file for figure 34

Authors’ original file for figure 35

Authors’ original file for figure 36

Authors’ original file for figure 37

Authors’ original file for figure 38

Authors’ original file for figure 39

Authors’ original file for figure 40

Authors’ original file for figure 41

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Rezania, V., Marsh, R., Coombe, D. et al. A physiologically-based flow network model for hepatic drug elimination I: regular lattice lobule model. Theor Biol Med Model 10, 52 (2013).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: