- Research
- Open Access

# A physiologically-based flow network model for hepatic drug elimination III: 2D/3D DLA lobule models

- Vahid Rezania
^{1}Email author, - Dennis Coombe
^{2}and - Jack A. Tuszynski
^{3}

**13**:9

https://doi.org/10.1186/s12976-016-0034-5

© Rezania et al. 2016

**Received: **17 September 2015

**Accepted: **22 February 2016

**Published: **3 March 2016

## Abstract

### Background

One of the major issues in current pharmaceutical development is potential hepatotoxicity and drug-induced liver damage. This is due to the unique metabolic processes performed in the liver to prevent accumulation of a wide range of chemicals in the blood. Recently, we developed a physiologically-based lattice model to address the transport and metabolism of drugs in the liver lobule (liver functional unit).

### Method

In this paper, we extend our idealized model to consider structural and spatial variability in two and three dimensions. We introduce a hexagonal-based model with one input (portal vein) and six outputs (hepatic veins) to represent a typical liver lobule. To capture even more realistic structures, we implement a novel sequential diffusion-limited aggregation (DLA) method to construct a morphological sinusoid network in the lobule. A 3D model constructed with stacks of multiple 2D sinusoid realizations is explored to study the effects of 3D structural variations. The role of liver zonation on drug metabolism in the lobule is also addressed, based on flow-based predicted steady-state O_{2} profiles used as a zonation indicator.

### Results

With this model, we analyze predicted drug concentration levels observed exiting the lobule with their detailed distribution inside the lobule, and compare with our earlier idealized models. In 2D, due to randomness of the sinusoidal structure, individual hepatic veins respond differently (i.e. at different times) to injected drug. In 3D, however, the variation of response to the injected drug is observed to be less extreme. Also, the production curves show more diffusive behavior in 3D than in 2D.

### Conclusion

Although, the individual producing ports respond differently, the average lobule production summed over all hepatic veins is more diffuse. Thus the net effect of all these variations makes the overall response smoother. We also show that, in 3D, the effect of zonation on drug production characteristics appears quite small. Our new biophysical structural analysis of a physiologically-based 3D lobule can therefore form the basis for a quantitative assessment of liver function and performance both in health and disease

## Keywords

## Background

The liver is a complex organ that removes chemicals, including drugs, from the blood through metabolic processes. When such liver detoxification processes are incomplete or overwhelmed, liver damage can result leading to liver failure. Drug access to hepatocytes is governed by transport processes in the well-vascularized liver tissue, and so structural variability can impact such transport. Therefore, a quantitative understanding of drug distribution and metabolism in the liver is essential for the ability to predict both liver performance and damage. We envisage a future generation of algorithms for anatomically- and physiologically-based liver models that can be then used for the prediction of optimal doses and scheduling of various drugs prior to clinical administration. This could be especially useful for patients with liver diseases such as hepatocellular carcinoma or cirrhosis.

Several investigators have recently explored computational fluid models of liver lobule function. Ierapetritou et al. [1] give a comprehensive general overview of the modelling approaches and issues that arose up to 2009. More particularly, Rani et al. [2] developed a detailed computational dynamics model of a small portion of a liver lobule, focusing on the non-Newtonian characteristics of blood. They considered the flow along one sinusoid with exit fenestrations fed by a portal vein and hepatic artery segments and exiting via a hepatic vein segment. The contributions of portal vein (PV) versus hepatic artery (HA) flows to the overall pressure drop and velocity profiles were detailed, including regions of Eddie flows and high strain rates near the exit fenestrations. Steady state flows were achieved after 5 × 10^{−5} s.

Yan et al. [3] developed a physiologically-based, multi-agent model of liver lobule performance. They used a Monte Carlo selection method to determine properties of the sinusoidal graph structure, convective-dispersive flow and metabolic interactions. This allowed inclusion of lobule zonation effects plus the metabolic influences of chemical molecular weight, octanol partition coefficients and protein binding to predicted drug-liver lobule interactions. They matched lobule outflow characteristics of 4 cationic drugs (atenolol, antipyrine, labetalol, and diltiazem), using sucrose as a base flow chemical. Wambaugh and Shah [4] compared the role of various lobule (sinusoidal) morphologies on drug propagation again using an agent-based simulation approach. While demonstrating that their model can reproduce traditional coarse averaging results (well-mixed and parallel tubes) under some flow regimes, they emphasized the utility of their random statistical approach for rapidly metabolized chemicals.

Hoehme et al. [5] employed a combined experimental and computational study of lobule structural restoration after CCl_{4} damage over several days as a protocol for liver regeneration. They used agent-based representations of both hepatocytes and sinusoids in their computational model and demonstrated that the regeneration process is characterized by a hepatocyte-sinusoid alignment process to correctly restore liver microarchitecture. Schliess et al. [6] extended this approach to consider ammonia detoxification during liver damage and regeneration. Here they first proposed a simple metabolic model considering ammonia, urea, and glutamine components utilizing three mass balance ordinary differential equations (ODEs) for chemical reactions and two compartments (periportal and perivenal) such that urea generation is maximized in the periportal region while glutamine regeneration occurs perivenally. The resulting concentrations were interpolated on their spatial-temporal grid as the sizes of the damaged zones change over the regeneration process. They also considered a multi-lobule pattern to reduce boundary effects. Drasdo et al. [7] summarized both modeling approaches in a review article.

In this paper, we extend our biophysical structural analysis to a full lobule situation, a hexagonal based model with one central blood in-flow (portal vein) and six corner located out-flows (hepatic veins), see Fig. 1. We also introduce the sinusoid network using a sequential diffusion-limited aggregation (DLA) algorithm (as presented and discussed below, Fig. 3) to gain insight on a more realistic resulting morphological sinusoid structure, see Fig. 1c. In a nutshell, sequential means a series of DLA generated pattern steps (5 steps) to generate the desired sinusoid structure variations in 2D. Thereafter we repeat this algorithm for various layers to generate a representative 3D lobule pattern. As a result, individual sinusoid layers, although generated by the same DLA algorithm, are not identical but exhibit some variability as observed in real lobules. A sensitivity analysis is conducted by observing drug concentration levels exiting the lobule with their predicted detail distribution inside the lobule.

### The liver lobule functional unit

It is commonly accepted that a single lobule serves as the functional unit of the liver, the smallest structural unit of the liver with ability to perform all hepatic functionalities [10]. 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 [11]. We shall invoke this second point of view as the basis for our model development.

Even among the approximately 1.5 million lobules (assuming a liver size of 1500 cm^{3} and a lobule size of 1 mm^{3}) that make up the human liver, structural variability of lobule units is the rule. Teutsch and colleagues [12] illustrate this specific microarchitecture variability and diseased states can be expected to add additional variability. Here we will attempt to quantify the consequences of such variability via our computational model.

## Methods

### 2D hexagonal lobule construction

The flow equations describing reactive-convective-diffusive flow in the liver lobule remain unchanged from our first two papers. Figure 1a and b illustrate the network structure of our base case model discussed in papers I [8] and II [9] (which uses a regular sinusoid pattern). Figure 1c represents our model that we study here, which consists of a 2D square lattice clipped by hexagonal boundary with DLA constructed sinusoid network. Figure 1d demonstrates a 3D extension of the lobule model that will be discussed later.

*N*

_{par}) with particle size (

*dp*×

*dp*in a pixel

^{2}). Figure 2a and b show result of DLA runs for

*N*

_{par}= 14,000 and

*dp*= 1 pixel and for

*N*

_{par}= 7000 and

*dp*= 2 pixel, respectively.

- i)A zero
*M*_{dim}×*M*_{dim}matrix (field) has been created and the DLA algorithm has been called to create a DLA pattern with specified*N*_{par}and*dp*. A particle, valued 1, has been located at a random position in the field by the DLA algorithm and then performed a random walk (with a desired step size, here denoted by*dp*) toward the center of the field. The newer particle will do the same until it hits another (older) particle, sticks to it and stops. This will continue for all*N*_{par}= 9500 particles. See Fig. 3a. Here*M*_{dim}= 256 and*dp*= 1. - ii)
As shown in Fig. 3a, the central region is very crowded. We overwrite the central region (here we chose half-the original field) using the DLA algorithm and with smaller number of zero-valued particles (blue dots). Here we use

*N*_{cen1_par}= 1200,*M*_{dim}= 128 and*dp*= 1. See Fig. 3b. - iii)
To fill the central region with sinusoids again (red dots), the DLA algorithm will be called for the third time with fewer particles, i.e.

*N*_{cen2_par}= 800. Here*M*_{dim}= 128 and*dp*= 2. See Fig. 3c. - iv)
In this step, we determine the six corners of the hexagon to create the out-flowing wells at these locations. Again the DLA algorithm will be used to create a DLA pattern around each corner. See Fig. 3d. Here we use the domain size of 1/3 of the original domain, i.e.

*M*_{dim}/3, and fewer particles,*N*_{cor_par}= 500. Here*dp*= 1. - v)
At the final stage, any point outside the bounding hexagon is removed. Figure 3e shows the final result if the DLA hexagonal lobule. As shown, hepatocytes are blue islands that are encompassed by red sinusoids.

This sinusoid generation algorithm is similar in spirit but not in details to that employed by Wambaugh and Shah [14]. They have used the term “sinusoid morphology” to characterize these investigations, which we will also employ. In recent work, Hoehme et al. [15] captured realistic lobule sinusoid patterns from image analysis of confocal microscopy of liver tissue. To check whether our simulated sinusoid pattern is comparable with their result, we perform a population density analysis for several 2D simulated and real lobules (Wisk S, Rezania V: A comparison between real and DLA simulated liver lobules using a population density analysis, Submitted). We calculated the ratio of the area covered by sinusoids to the total area of the lobule for both simulated and real 2D lobule images to optimize values of *N*
_{par} and *dp* for corresponding *M*
_{dim}. Based on these calculations, we find *N*
_{par} ≈ 37 *M*
_{dim} (~9500) and *dp* = 1 will produce comparable results (~45 % sinusoid area/total 2D area). A comparison of our Fig. 3d with their results demonstrates that our DLA algorithm with the chosen parameters generates equally physiologically reasonable sinusoid morphology patterns (in particular, see their Fig. 1d or further details in their SI Appendix, their Fig. 3d).

### Extension to a 3D lobule structure

In order to generate an even more realistic survey, we next study a 3-dimensional lobule. Again our morphological sinusoid generation algorithm is similar in spirit to that employed by Wambaugh and Shah [14] in their 3D models.

- i)
A vertical dimension

*M*_{z-dim}will be specified. - ii)Starting from the very bottom layer that has a 2D DLA structure, we skip
*M*_{skip}layers to introduce the next layer with a 2D DLA structure (a DLA layer), see Fig. 4a. This procedure will be repeated to reach to the very top layer which again is a DLA layer. Here, each DLA layer represents the sinusoidal network. Two of these layers encompass hepatocyte cells (the skipped layers). The value of*M*_{skip}is chosen based on the average thickness size of a hepatocyte cell. - iii)
Then, two DLA layers vertically will be connected at several random locations to approximate the shape hepatocyte cells, see Fig. 4b.

- iv)
The center and the six corners of the hexagon in all layers will be connected to represent arterial and portal veins.

Since the skipped layers represent the hepatocyte tissue, the value of *M*
_{skip} will be determined based on the typical volume of a lobule and hepatocytes. The typical volume of a human liver lobule is somewhat less than 1 mm^{3} [12]. Furthermore, a hepatocyte has a diameter of 12 – 24 μm and thickness of 25 μm in average and mean sinusoid thickness is about 5 – 7 μm. Here, we chose a sinusoidal cell as a cube with dimension of 6 μm and a hepatocyte cell as rectangular cube with dimensions 6 μm and 12 μm, with thickness of 6 μm in the DLA layers and 12 μm in the skipped layers, respectively. The area of a hexagon with diameter of *d* is given by \( S=3\sqrt{3}{d}^2/8 \).

In our case, *d* = 256 × 6 μm = 1536 μm that leads to *S* = 1.53 × 10^{6} μm^{2}. The 3D lobule is then composed by stacking up 33 sinusoidal layers and 64 tissue layers (*M*
_{z-dim} = 97) that gives an approximate model volume of 1.5 mm^{3}. Here, every two tissue layers are sandwiched by two sinusoidal layers (*M*
_{skip} = 2). This choice of tissue/sinusoid grid thickness represents a compromise between simulation runtime speed and numerical discretization error.

Figure 4 shows the 3D DLA hexagonal lobule for *M*
_{z-dim} = 97 and *M*
_{skip} = 2. Figure 4a demonstrates a 2D DLA layer representing a sinusoidal layer. Figure 4b shows a hepatocyte layer. The red points connect all hepatocyte layers between two DLA layers. These points approximately determine the boundary of the hepatocyte cells. To construct the third dimension, a DLA layer is created for the bottom, then two hepatocyte layers (*M*
_{skip} = 2), and then a new DLA layer will be produced. The procedure repeats over the whole *z*-dimension. Figures 4c and d show two different 3D views.

### Flow calculation methods

Base case flow and metabolism parameters

Parameter | Characteristic (SI) Unit | STARS Unit |
---|---|---|

Sinusoid Porosity | 0.7854 | 0.7854 |

Sinusoid Permeability | 1.125 μm | 1.140 Darcy |

Sinusoid Effective Diffusion | 4.2 × 10 | 2.5 × 10 |

Tissue Porosity | 0.2382 | 0.2382 |

Tissue Permeability | 7.35 × 10 | 7.45 × 10 |

Tissue Effective Diffusion | 4.2 × 10 | 2.5 × 10 |

Maximum Rate | 0.06 μM/min | 1.08 × 10 |

Half Saturation Constant | 10.0 μM | 1.8 × 10 |

Linear Rate | 6.0 × 10 | 6.0 × 10 |

Blood Viscosity | 3.5 × 10 | 3.5 cpoise |

## Results and Discussion

### 2D drug distribution – comparison to our earlier studies

^{−6}cm

^{3}/min with (0.4077, 0.3011, 0.1950, 0.9682, 0.2251, 2.394) × 10

^{−6}cm

^{3}/min for ports A to F, respectively, as illustrated in Fig. 5. (For steady flows, the inflow rate equals the total outflow rate.) It is clear that ports D and F have highest outflow among the others. This is due to the random nature of the DLA simulation for this realization.

^{−8}mole fraction) is infused into the lobule through the central inlet. Assuming nonreactive hepatocytes, the time required to traverse the lattice is approximately 1 min or less as demonstrated in Fig. 7. The fastest drug propagation is through ports D and F for this realization. Similar to our Paper I [8], this production profile is convective flow dominated as the addition of diffusion minimally alters the production profile (figure not shown).

Flow rate variability of multiple 2D DLA sinusoid morphologies

Layer no. | SS flow rate (×10 |
---|---|

Lyr1 | 1.824 |

Lyr4 | 1.366 |

Lyr7 | 1.829 |

Lyr10 | 0.776 |

Lyr13 | 1.781 |

Lyr16 | 1.716 |

Lyr19 | 1.399 |

Lyr22 | 0.889 |

Lyr25 | 1.555 |

Lyr28 | 1.225 |

Lyr31 | 1.684 |

Lyr34 | 1.416 |

Lyr37 | 1.251 |

Lyr40 | 1.044 |

Lyr43 | 1.421 |

Lyr46 | 1.983 |

Average | 1.447 |

STDEV | 0.3511 |

^{−4}cm

^{2}/min in sinusoid and 2.5 × 10

^{−5}cm

^{2}/min in tissue, respectively. In the absence of diffusion, only a small “pressure-difference-driven (convective)” transfer from sinusoids to tissue through the space of Disse is occurring. See our Paper I [8] for details.

Figures 8 and 9 show the increasing levels of injected drug from 0.01 min to 5 min (0.5 min in case with diffusion). As shown, by 0.01 min, PAC almost reaches the outgoing port F of the lobule.

In case with no diffusion, it takes almost 5 min to fill the lobule (Fig. 8), while PAC completely covers the lattice after 0.5 min when diffusion is on as demonstrated in Additional file 1: Figure S1.

### Reactive flows in 2D

Now we consider the effects of PAC drug metabolism by hepatocytes. Here the base case reaction parameters of Table 3 of Paper I [8] are employed, and the same injected PAC 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 equation.

Figure 9 illustrates injected drug and produced drug and metabolite production for this case. The reaction rate is 6 × 10^{−3} min^{−1}. 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 PAC is converted to metabolite by the lobule hepatocytes. The production profile of PAC-OH here is similar to the production profile of PAC in the non-reaction case, as shown in Fig. 7, for most of the ports. That is, the reactive conversion of PAC to PAC-OH occurs reasonably quickly. However, for the two ports with higher flows (ports D and F), less conversion of PAC to PAC-OH is seen to occur.

^{−6}min

^{−1}and 6 × 10

^{−3}min

^{−1}, with and without diffusion. (The latter rate is our base case which has the higher conversion rate.) As shown, at each rate a greater concentration of PAC-OH is generated when diffusion is on, as this provides an additional mechanism to bring PAC molecules to the reactive hepatocyte sites.

### 3D drug distribution – extension of behaviors

As stated above, the 3-dimensional lobule is constructed by stacking up 33 sinusoidal layers and 64 hepatocyte layers (*M*
_{z-dim} = 97) with the approximate volume of 1.5 mm^{3}. The sinusoidal layers have similar structures (randomly generated though) as the 2D DLA lobule structure discussed in previous section. The hepatocyte layers are generated by randomly selected points to mimic hepatocyte cell distribution in a lobule. All layers are connected through the central region (the portal vein) as well as six corners of the hexagon (hepatic veins). In this section we are interested in studying the effect of the third dimension as the blood will be injected from the central vein in all layers simultaneously.

^{−3}cm

^{3}/min as shown in Fig. 12a, as all volumes have been upscaled from our 2D slice models. As seen in Fig. 12b and c, in this realization of the 3D lobule, ports A and F have highest outflow (they are essentially superimposed), then ports B and D (they are essentially superimposed), with C and E having the lowest outflow rates. This is again due to the random nature of the DLA simulation for this realization. However the port-to-port variation in flows is significantly less extreme in 3D, as the third dimension provides additional connecting flow paths to smooth the overall flow pattern.

The same two views of non-reactive PAC injection with diffusion added are shown in Additional file 3: Figure S2 and Additional file 4: Figure S3. In those figures the sinusoid-tissue transfer rate is much more rapid and PAC drug propagates in a more uniform fashion as it enters the lobule. Clearly, diffusion has a very significant effect on the drug propagation details although at later times both cases result in a uniform drug coverage of the lobule.

### Reactive flows in 3D

Now we consider the effects of PAC drug metabolism by hepatocytes across the 3D lobule. Similar to the 2D case, the same injected PAC concentration (1.8 × 10^{−8} mole fraction) is considered.

^{−3}min

^{−1}. As before, both PAC and PAC-OH have assumed equal diffusive flow contributions.

Figures 14a and b demonstrate the levels of PAC and PAC-OH, respectively, without including diffusion effects. As shown, both PAC and PAC-OH have similar concentrations after 1 min in all outlets (i.e. approximately half of injected PAC is metabolized to PAC-OH). The diffusion, however, alters the latter behavior significantly as depicted in Fig. 14c and d. All the PAC converted to PAC-OH almost immediately. The variability of the produced PAC-OH concentration per port for this reactive case with diffusion mirrors the variability of the produced PAC per port for the nonreactive case (Fig. 11, with or without diffusion).

As shown in Additional file 6: Figure S4, with diffusion but still utilizing a base case metabolic rate, the in-situ conversion of PAC to PAC-OH is much more rapid and uniform, and results in an almost total drug conversion to PAC-OH throughout the 3D lobule. There remains only a small inlet region of unconverted PAC.

As in the 2D cases presented earlier, diffusion has greater impact on drug distribution for the reactive cases. Similar behaviors can be observed for other reaction rates (e.g. reaction rate 6 × 10^{−1} min^{−1}) with and without diffusion (not shown).

### 3D drug distribution with zonation

Zonation is a well-known feature of many metabolic processes in the liver lobule [15], some processes are up-regulated in the periportal region, while others are up-regulated in the perivenal region. Here we consider zonation of drug metabolizing enzymes (in particular cytochrome P450 – CYP2C8), such that higher CYP levels are found near the perivenal region.

Zonation has been attributed primarily to a non-uniform distribution of O_{2} across the lobule [16]. Here we will utilize this experimental observation to predict relative CYP levels based on a calculated O_{2} distribution. Appendix A presents details of the O_{2} convective-diffusive-reactive flow problem employed to generate the zoned-CYP distribution.

The resultant initial non-homogeneous distribution of CYP enzymes is imported into the PAC reaction model in a manner analogous to that employed in our earlier idealized model of zonation [9]. Here we consider enhanced perivenal enzyme expression (the “reversed” case of reference [9]) as a reference, where CYP expression is active primarily near the lobule outlet zone. In our present work, two scenarios are envisioned: (a) a CYP distribution generated from “normal” levels of O_{2} injection, and (b) that generated from a “low concentration” O_{2} injection level. The result is shown in Figures A3 and A4, which can be contrasted to Additional file 3: Figure S2b of our earlier idealized model [9]. These distributions are used to modify the PAC metabolic rates non-homogeneously in a manner similar to our earlier paper.

Spatial plots of PAC and PACOH distributions at various time points for the no diffusion cases with zonation show only subtle differences from the base case no-zonation profiles of Fig. 16 and thus are not shown explicitly. This is perhaps not surprising for the low O_{2} concentration generation case with a very uniform zoned CYP distribution, but is more puzzling for the normal O_{2} concentration generated case. We believe this is a result of the high PAC to PACOH conversion rate chosen as our example for this paper coupled to the downstream distribution of CYP. This example appears to be quite robust to details of the exact CYP zonation pattern.

^{−6}min

^{−1}). Again only ports A and C are shown. For these cases, the low concentration zonation production is again essentially identical to the no zonation case, but the normal zonation production differs more substantially from the other two cases. This behaviour is reflected in the spatial distributions of PAC and PACOH at various time points as well. Figure 18 compares the final (at 1 min) spatial distributions for PAC and PACOH for the normal and low-oxygen concentration generated CYP distributions. (The no zonation case is essentially identical to the latter case plots, not shown). The normal oxygen generated CYP distribution results in less PACOH conversion.

## Conclusions

### Generalized conclusions relative to our earlier work

- (a)
In 2D due to randomness of the sinusoid patterns, individual producing ports (hepatic veins) respond differently (i.e. at different times) to injected drug. This is similar to the randomness we generated via multiple realizations to our individual well pair models with our previous idealized modelling work. If we sum the individual hepatic vein responses we get an overall 2D lobule drug response (see Fig. 5), the net result of which is an increased spreading (effective diffusion) of the average produced profile.

- (b)
Multiple 2D morphological sinusoid realizations can be explored to ascertain the effects of structural variations, as indicated in Table 2. Additionally a 3D model viewed as interacting stacks of 2D sinusoid realizations represents an enhanced method to consider such variability.

- (c)
In 3D, because the individual producing ports are completed throughout the vertical extent of the lobule, the variation of response of individual producing ports is seen to be less extreme, although the production curves themselves are more diffuse in 3D than in 2D. Furthermore, the average lobule production summed over all hepatic veins is also more diffuse. Thus the net effect of all these variations is to make the overall response more smooth (i.e. more robust to individual variations).

- (d)
When the effect of drug processing (CYP) zonation is included, the effect on drug production characteristics appears quite small and only when the further (smoothing) effects of diffusion are ignored. However, within a lobule we expect there are conditions under which there can be noticeable differences in drug distribution due to zonation. These observations are consistent with our earlier idealized models of zonation. Finally we should emphasize that this treatment of zonation of drug metabolism can be readily extended to other metabolic zonation phenomena occurring in the liver, such as carbohydrate metabolism [21], nitrogen metabolism [22], etc. O

_{2}distribution is implicated as a fundamental cause of such zonation in all cases.

### Follow-up work and future directions

Since various liver diseases can be thought to produce structural variations in the lobule, our analysis also gives insight into the role of disease on liver function and performance. In our previous work on variable lattice lobule models [9], we utilized random sinusoid permeability models and percolation concepts to capture some of the structural effects of hepatitis and cirrhosis. We next plan on extending such analysis to our more physiologically realistic 3D lobule model.

It should be mentioned that there exist other computer-generated liver models that attempt to provide a realistic representation of the liver vasculature for blood flow and metabolic reaction simulations. Schwan et al. [23] used an image of a mouse liver produced with a CT scanner. Based on this image data, they reconstructed the structure of the fine branches of the liver vessel system. The liver was then split into 50,000 small blocks. The results of the simulation show that blood flow and metabolic reactions can be tracked in detail on the computer monitor. In a follow-up publication [24] a multi-scale model was generated that links together four scales of modeling. Thus the focus of these authors is the full liver scale (organ model) and not the lobule scale (tissue level) which is the subject of our current paper. In addition, we have also recently developed a spatial model of the full liver (White D, Rezania V, Coombe D, Tuszynski J: Building a 3D virtual liver: a preliminary example – methods for describing vascular generation, blood flow calculations, and hepatic clearance, submitted) utilizing a computational algorithm for generating vasculature and where we have considered upscaling from our lobule model and some comparisons with the work of Schwen et al. [23]. Eventually, our overall modeling approach is intended to provide input data regarding liver metabolism for multi-compartment physiologically-based pharmacokinetic models such as those reviewed by Jones et al. [25]. This could lead to more accurate and individualized predictions for the pharmacokinetic analyses of drugs and drug combinations. It is also of interest to explore the possibility of fractal vascular structures and their consequences on the scaling laws in liver function as reviewed extensively by Pang et al. [26].

Our focus will continue to be on the structural and metabolic changes induced by liver disease and cancer, and its remediation via drug treatments. More particularly, we will focus on drug combination therapy for cancer (e.g. paclitaxel with doxorubicin) [27].

## Declarations

### Acknowledgements

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

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- 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(3):228.View ArticlePubMed CentralGoogle Scholar
- Rani HP, Sheu T, Chang TM, Liang PC. Numerical investigation of non-Newtonian microcirculatory blood flow in hepatic lobule. J Biomechanics. 2006;39:551.View ArticleGoogle Scholar
- Yan L, Ropella GE, Park S, Roberts MS, Hunt CA. Modelling and simulation of hepatic drug disposition using a physiolocally based. Multi Agent In Silico Liver Pharamceutical Res. 2007;25:1023.Google Scholar
- Wambaugh J, Shah I. Simulating microdosimetry in a virtual hepatic lobule. J Biomechanics. 2010;6(4):e1000756.Google Scholar
- Hoehme S, Brulport M, Bauer A, Bedawy E, Schormann W, Gebhardt R, et al. Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration. PNAS. 2010;107:10371.View ArticlePubMedPubMed CentralGoogle Scholar
- Schliess F, Hoehme S, Henkel SG. Integrated metabolic spatial-temporal model for the prediction of ammonia detoxification during liver damage and regeneration. Hepatology. 2014;60:2040.View ArticlePubMedGoogle Scholar
- Drasdo D, Hoehme S, Hegstler JG. How predictive quantitative modelling of tissue organization can inform liver disease pathogenesis. J Hepatology. 2014;61:951.View ArticleGoogle Scholar
- Rezania V, Marsh RE, Coombe D, Tuszynski JA. A physiologically-based flow network model for hepatic drug elimination I: regular lattice lobule model. Theor Biol Med Model. 2013;10:52. doi:10.1186/1742-4682-10-52.(PaperI).View ArticlePubMedPubMed CentralGoogle Scholar
- Rezania V, Marsh RE, Coombe D, Tuszynski JA. A physiologically-based flow network model for hepatic drug elimination II: variable lattice lobule models. Theor Biol Med Model. 2013;10:53. doi:10.1186/1742-4682-10-53.(PaperII).View ArticlePubMedPubMed CentralGoogle Scholar
- Saxena R, Theise ND, Crawford JM. Micro-anatomy of the human liver – exploring the hidden interfaces. Hepatology. 1999;30:1339.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–7.View ArticlePubMedGoogle Scholar
- Teutsch HF. The modular microarchitecture of human liver. Hepatology. 2005;42:317.View ArticlePubMedGoogle Scholar
- MATLAB: www.mathworks.com.
- CMG. Ltd. STARS User’s Guide: advanced process and thermal reservoir simulator. Calgary, AB: Computer Modelling Group Ltd; 2014.Google Scholar
- Gebhardt R. Metabolic zonation of the liver – regulation and implications for liver function. Pharmacol Ther. 1992;53(3):275.View ArticlePubMedGoogle Scholar
- Jungermann K, Thurman R. Oxygen: modulator of metabolic zonation and disease in the liver. Hepatology. 2000;31(2):255.View ArticlePubMedGoogle Scholar
- Rowinsky EK, Wright M, Monsarrat B, Lesser GJ, Donehower RC. Taxol: Pharmacology, metabolism, and clinical implications. Cancer Surv. 1993;17:283–304.PubMedGoogle Scholar
- Huizing MT, Misser VH, Pieters RC, ten Bokkel Huinink WW, Veenhopf CH, Vermorkem JP, et al. Taxanes: a new class of antitumour agents. Cancer Invest. 1995;13:381–404.View ArticlePubMedGoogle Scholar
- Vaclavikova R, Soucek P, Svobodova L, Anzenbacher P, Simek P, Guengerich F, et al. Different in vitro metabolism of paclitaxel and docetaxoel in humans, rats, pigs, and minipigs. Drug Metab Dispos. 2004;32(6):666–74.View ArticlePubMedGoogle Scholar
- Monsarrat B, Chatelut E, Royer I, Alvinerie P, Dubois J, Dezeus A, et al. Modification of paclitaxol metabolism in cancer patient by induction of cytochrome P450 3A4. Drug Met Disp. 1998;26:229–33.Google Scholar
- Jungermann K, Thurman R. Hepatocyte heterogeneity in the metabolism of carbohydrates. Enzyme. 1992;46:33.PubMedGoogle Scholar
- Haussinger D, Lamers W, Moorman A. Hepatocyte heterogeneity in the metabolism of amino acids and ammonia. Enzyme. 1992;46:72.PubMedGoogle Scholar
- Schwen LO, Krauss M, Niederalt C, Gremse F, Kiessling F, et al. Spatio-temporal simulation of first pass drug perfusion in the liver. PLoS Comput Biol. 2014;10(3):e1003499. doi:10.1371/journal.pcbi.1003499.View ArticlePubMedPubMed CentralGoogle Scholar
- Schwen LO, Schenk A, Kreutz C, Timmer J, Bartolomé-Rodríguez MM, Kuepfer L, et al. Representative sinusoids for hepatic four-scale pharmacokinetics simulations. PLoS ONE. 2015;10(7):e0133653. 1–39.View ArticlePubMedPubMed CentralGoogle Scholar
- Jones HM et al. A novel strategy for physiologically based predictions of human pharmacokinetics. Clin Pharmacokinet. 2006;45(5):511–42.View ArticlePubMedGoogle Scholar
- Pang KS, Weiss M, Macheras P. Advanced pharmacokinetic models based on organ clearance, circulatory, and fractal concepts. AAPS J. 2007;9(2):E268–83.View ArticlePubMedPubMed CentralGoogle Scholar
- Vos KJ, Martin AG, Trimbol MG, Forestell L, Barakat K, Tuszynski JA. A multi-compartment pharmacokinetic model of the interaction between paclitaxel and doxorubicin. EPJ Nonlinear Biomed Phys. 2014;2:13.View ArticleGoogle Scholar
- Davidson AJ, Ellis MJ, Chaudhuri JB. A theoretical method to improve and optimize the design of bioartificial livers. Biotech Bioeng. 2010;106(6):980.View ArticleGoogle Scholar
- Davidson AJ, Ellis MJ, Chaudhuri JB. A theoretical approach to zonation in a bioartificial liver. Biotech Bioeng. 2011;109(1):234.View ArticleGoogle Scholar