The liver is the central organ for detoxification of xenobiotics in the body. In pharmacokinetic modeling, hepatic metabolization capacity is typically quantified as hepatic clearance computed as degradation in well-stirred compartments. This is an accurate mechanistic description once a quasi-equilibrium between blood and surrounding tissue is established. However, this model structure cannot be used to simulate spatio-temporal distribution during the first instants after drug injection. In this paper, we introduce a new spatially resolved model to simulate first pass perfusion of compounds within the naive liver. The model is based on vascular structures obtained from computed tomography as well as physiologically based mass transfer descriptions obtained from pharmacokinetic modeling. The physiological architecture of hepatic tissue in our model is governed by both vascular geometry and the composition of the connecting hepatic tissue. In particular, we here consider locally distributed mass flow in liver tissue instead of considering well-stirred compartments. Experimentally, the model structure corresponds to an isolated perfused liver and provides an ideal platform to address first pass effects and questions of hepatic heterogeneity. The model was evaluated for three exemplary compounds covering key aspects of perfusion, distribution and metabolization within the liver. As pathophysiological states we considered the influence of steatosis and carbon tetrachloride-induced liver necrosis on total hepatic distribution and metabolic capacity. Notably, we found that our computational predictions are in qualitative agreement with previously published experimental data. The simulation results provide an unprecedented level of detail in compound concentration profiles during first pass perfusion, both spatio-temporally in liver tissue itself and temporally in the outflowing blood. We expect our model to be the foundation of further spatially resolved models of the liver in the future.
The liver is the main site of metabolization and detoxification of xenobiotics in the body of mammals. Compounds delivered by blood flow through the portal vein and the hepatic artery are continuously processed within hepatic cells, such that foreign and potentially harmful compounds can be cleared from the blood. Metabolization by enzyme-catalyzed biotransformation produces chemical alterations of the original compounds, thereby enabling their elimination. A second, complementary process is the active secretion to the bile duct from which the compound is further transported to the gastrointestinal tract. In pharmacology and medicine, plasma clearance is used to quantify the rate by which a compound is eliminated from the body [1]. Plasma clearance describes the overall detoxification capacity of the organism and summarizes individual clearance rates from all eliminating organs with the largest contribution coming from the kidney and the liver. While renal clearance can be measured by urinary secretion, a quantification of liver detoxification capacity is difficult since the different physiological functions cannot be assessed directly. In particular, the relative contributions of the different underlying physiological functions such as metabolization or biliary secretion cannot be adequately differentiated since the liver is frequently rather considered as a net systemic sink. While hepatic turnover can be indirectly quantified with drugs following a known pharmacokinetic profile, the local, time-resolved distribution of compounds within the whole organ can generally not be analyzed even with distinguished measurement techniques. This holds in particular for the first pass of drug perfusion in a liver directly after drug administration, when hepatic tissue is exposed to a novel xenobiotic for the first time. Due to the pivotal role of the liver in drug pharmacokinetics and detoxification, several models quantifying hepatic clearance have been developed before [2]. These include tube models for representative hepatic sinusoids (single tubes, in parallel or in series) [3], dispersion liver models [4], fractal liver models [5], circulatory models [6], [7] including zonal models [8], [9], and distribution-based models describing statistical variations in transition times [10]. For a more detailed overview, we refer to [11]. Generally, PK models are well suited to investigate distribution and clearance of compounds in the body. In compartmental PK modeling, a limited number of rather generic compartments is usually used to simulate plasma concentration and drug clearance. Following a complementary approach, physiologically based pharmacokinetic (PBPK) models describe biological processes at a large level of detail based on prior physiological knowledge [12]. This involves amongst others organ volumes and blood flow rates, such that physiological mechanisms underlying absorption, distribution, metabolization and excretion of compounds can be explicitly described. While most approaches consider intravenous or oral application of therapeutic compounds, PK models describing further uptake routes such as inhalation have also been developed [13], [14]. Organs in PBPK models are divided in several subcompartments. So-called distribution models are used to describe the mass transfer between these subcompartments which are quantified based on physicochemical compound information such as lipophilicity or the molecular weight. Basic PBPK models can be extended to include enzyme-mediated metabolization or active transport across membranes [15]. Coming from the field of toxicology [16], PBPK models are nowadays routinely used in drug discovery and development [17]. They are for example applied in pediatric scaling [18], model-based risk assessment [19], as well as for multiscale modeling by integrating detailed models of intracellular signaling [20] or metabolic networks [21] into the cellular subcompartment, thus allowing for the analysis of cellular behavior within a whole-body context [22]. Notably, each organ in PBPK modeling is generally represented by few well-stirred subcompartments, thus allowing a quantitative description of drug pharmacokinetics once an equilibrium between the vascular system and the surrounding tissue has been reached. However, a spatially resolved description of drug perfusion in the whole organ covering particularly the first instants after drug administration is impossible due to the inherent well-stirred assumption. To mechanistically describe first pass perfusion, we here present a spatio-temporal model of drug distribution and metabolization in a mouse liver. The model represents liver lobes at the spatial length scale of lobuli such that physiological and molecular processes can be simulated in great detail. Our combined spatially resolved model (cf. overview in Figure 1) involves three main building blocks. These comprise physiological vascular structures, an organ-scale perfusion model describing blood flow and advection of compounds, and finally pharmacokinetic models translated to the spatially resolved tissue. Geometrically accurate models of murine hepatic vascular structures were obtained by using in vivo micro-CT imaging [23]. The mass balance within the tissue was quantified based on equations coming from PBPK modeling. Our combined model was inspired by [24], where a lobule-scale perfusion model in more physical detail and also allowing for deformation of the porous medium is introduced. A model for cardiac perfusion using very similar modeling techniques is presented in [25]. It considers multiple geometric scales, but no draining vascular system and no metabolization. Our spatially resolved model covers several scales of biological organization displayed at varying levels of resolution. The scales range from the organ level to the cellular space where metabolization and molecular transport take place. The vascular systems form the scaffolds which links the hepatic in-flow to the sinusoidal space and thereby to the lobulus level. The model considers one supplying and one draining vascular system (denoted by SVS and DVS, respectively), with a homogenized hepatic space (denoted by HHS) in between as a tissue representation. The HHS includes in particular the sinusoids which are not explicitly resolved as vascular structures. In the latter, blood flow is represented by a fluid transport model [26]. Microcirculation and microanatomy [27] are only considered in effective form. While more detailed perfusion and metabolization models on the lobular scale [28], [29] or the tissue-level [2] have been developed before, the representation of the HHS as a porous medium [30] is sufficient for our needs. Experimentally, the combined model corresponds to an isolated perfused liver [31], [32], since recirculation of blood is not considered. The resolution of the model allows calculating local concentration profiles within the tissue which can for example be used to address heterogeneous phenomena such as spatially varying distributions of lipid droplets in steatotic livers. Spatial heterogeneity of pharmacokinetic parameters such as metabolic capacity can be taken into account. Likewise local exposure profiles of toxic compounds can be simulated such that off-target effects of drug therapy can be analyzed at a high level of resolution. Applications of spatially resolved perfusion and metabolization modeling include optimized design of therapeutic treatment where spatio-temporal perfusion effects are of relevance, e. g. hypothermic machine perfusion [33] or islet cell transplantations [34]. Moreover, such models can be used for planning drug delivery for which spatially heterogeneous distribution is an important property and crucial for administration design itself. Two major examples are intrahepatic injection of chemotherapeutics or radionuclides (see e. g. [35]), in particular in combination with optimization of catheter placement [36], and targeted drug delivery [37] where drugs are injected in bound form and released at the desired location by mild hyperthermia. Likewise, our model may support data processing and interpretation in imaging or diagnostics. We expect the spatially resolved model presented here to be the foundation stone of further mechanistic models describing the spatial organization of the liver in an unprecedented level of physiological detail.
The animal experiment was reviewed and approved by the local authorities (NRW LANUV, permit number 10576G1) according to German animal protection laws. In the methods to be presented, we follow a geometric view from coarse to fine, i. e. from (1) the body scale (providing organ input and output) via (2) the vascular structures on the organ scale (perfusion only) to (3) the tissue scale (perfusion and metabolization). Models for steatosis and -induced liver necrosis are subsequently presented to demonstrate how our spatially resolved simulations can be used for the analysis of pathological states influencing drug distribution in the liver. Finally, some aspects of computational resolution are addressed. Our spatially resolved model distinguishes between the supplying vascular tree, the HHS and a draining vascular tree which are considered in series (Figure 2). For reasons of simplicity, the supplying vascular system comprises both the portal vein and the hepatic artery. Since an isolated perfused liver was considered here which explicitly excludes recirculation through the body, the respective contributions of the portal vein and the hepatic artery were not distinguished and only the total blood inflow was taken into account. The draining vascular system represents the hepatic vein. Bile ducts were not considered in our model, since their geometric structure could not be resolved experimentally. In our combined model, the HHS is composed of several subspaces in analogy to the liver compartment in PBPK models. The latter is divided in four subcompartments, i. e. red blood cells (rbc), plasma (pls), interstitium (int), and liver cells (cell). Those four subcompartments are also considered as subspaces of the HHS, in addition a fifth remaining subspace (rest) is taken into account. This remaining subspace comprises all those parts inside the liver that are not considered for perfusion, metabolization and active transport, in particular the larger and explicitly resolved vascular structures and the bile ducts. The plasma subspace only refers to the blood plasma. For notational convenience, the sinusoidal subspace (sin) is defined as the combination of red blood cells and plasma subspace, thus representing a whole-blood compartment. The sinusoidal subspace is subject to advection, thereby reflecting blood flow through the tissue. The actual metabolization takes place only in the cellular subspace and is part of the PBPK equations that also model the exchange between the HHS subspaces. The vascular trees are resolved separately and considered for pure advection. The volume fractions of the subspaces relative to the overall liver volume, also needed for the mass balance in the compartment models, are denoted by. For our simulations in a mouse liver, we use the values (1) Volume fractions were obtained by setting to cover the fraction of the vascular volume inside the segmented liver volume, both determined from the experimental image data as described below. The values for the other subspaces from [38] were then scaled accordingly by so that all five volume fractions sum up to. From Equation 1, we immediately obtain. For the perfusion part in our model, we address how molar concentrations of compounds change due to advection through the vascular systems and the sinusoidal HHS subspace. The exchange with the remaining HHS subspaces and cellular metabolization are considered as a separate contribution to our combined model. The body scale determines the total perfusion of the liver in mice [38]. The blood flow into and out of the root edges of the supplying and draining vascular system, respectively, is the connection of the HHS to the surrounding organism. More precisely, inflowing and outflowing molar concentrations of the compounds of interest are the main model input and output quantities. The HHS is modeled as a porous medium representing the effective behavior of the whole liver volume, with the sinusoidal subspace being perfused according to Darcy' s law [51]. The perfusion is split between the red blood cells and plasma subspaces (see Figure 2) proportional to their respective volumes. Using appropriate flow source and sink terms in the HHS corresponding to the exchange with the vascular systems, flow velocities for 3D advection in the HHS are obtained. The advection of concentrations can then be calculated from the given flow velocities using appropriate discretizations described in this section. Technical details about the discretization and implementation in the model are presented in Section 2 in Text S1. Finally, the pointwise exchange in the spatially resolved model between different HHS subspaces and the actual metabolization are modeled by equations coming from PBPK modeling. The final spatially resolved model can also be used for analysis of pathophysiological states of the liver which have not been taken into account during model establishment itself. We here considered the case of steatosis leading to changes in intracellular lipid content as well as carbon tetrachloride () -induced liver necrosis affecting the spatial composition of the organ. Describing pathophysiological changes in spatially heterogeneous patterns is a key strength of our approach. A comparison of the simulation results with experimental data allows to evaluate model validity, thereby providing targeted suggestions for model extensions and modifications. Steatosis is a common liver disease often caused by obesity or alcohol abuse [57]. It is characterized by lipid accumulations in the cellular subspace [58], the influence of which can be structurally represented in the model. We here analyzed to what extent steatosis affects hepatic clearance following changes in intrahepatic drug distribution. For our simulations we consider data reported from rats in [59, Table 8], namely steatosis extents of about and (mean SD) in the left lateral and median liver lobe, respectively, obtained after two weeks of a specific diet. We proceed assuming that similar steatosis patterns can also exist in mice. Let be the ratio of lipid accumulation per liver volume at position, corresponding to the steatosis percentages in [59, Table 8]. Using the lobe decomposition of our mouse liver dataset (cf. Figure 3), we consider two states of steatosis. First, we use a homogeneous lipid accumulation throughout the liver. This value is obtained as the average for in the left lateral lobe and in the remaining lobes as the left lateral lobe in our case has of the total volume. Second, we assign a pseudo-randomly varying value uniformly distributed in (left lateral lobe) and (remaining lobes) to obtain a spatially heterogeneous steatosis pattern. To avoid unphysiologically large local variations, we generated random numbers [60] on a grid four times coarser than the computational resolution and interpolated multilinearly at the nodes actually used. The two steatotic cases are visualized in Figure 4. In this setting, corresponds to the healthy state [12]. We quantify the impact of steatosis (an increased intracellular volume fraction of lipids) by changing in cellular partition coefficient in the PBPK model. Any other effects of steatosis are explicitly omitted here for the sake of simplicity of this proof of concept for a spatial heterogeneity. The cellular partition coefficient is calculated according to the formula [12] (10) with a constant specific for the respective compound and (Table 1). The values are substituted in Equations 7 and 8a/8b. As varies spatially, also and thus intracompartemental exchange and metabolization vary accordingly. This is in contrast to commonly used PBPK models that, due to their compartmental organ representation, cannot distinguish between the homogeneous and heterogeneous case as they only use one constant value of. As a second example for pathophysiological changes in the liver, we consider the case of -induced hepatic injury. Administration of in animal models is a frequently used experimental protocol to investigate the processes underlying toxic liver damage [61]. Inducing hepatic injury by leads to necrotic cell death in the pericentral zone, similar to acetaminophen overdoses [62]. We analyzed the impact of pericentral necrosis on hepatic metabolization capacity. In our spatially resolved model, -induced necrosis was represented by replacing the cellular space by interstitial space in the necrotic volume. The latter was set to be the of the liver volume closest to the DVS terminal edges (see Figure 4), where the percentage is based on the area analysis in [63, Table 1], observed one day after administration. Our basis for choosing computational resolutions is the actual size of lobuli in mice. From a cross-section area of, a radius of (assuming a regular hexagonal shape), both from [63, Table 1], and assuming the same elongation (length divided by diameter) of as for human lobuli [64, Chapter 2. 5], a mouse lobulus has a volume of approximately, the total liver volume of thus corresponds to lobuli. By definition, a lobule is the volume drained by one terminal edge of the hepatic vein, so a fully resolved vascular tree has approximately that many leaf nodes.
As a first application example without metabolization, we considered the distribution of the tracer CFDA SE within the liver. Since adequate pharmacokinetic data for CFDA SE were not available for mice, a PBPK model could not be validated in detail. Instead, only the basic physicochemical parameters (and) were estimated and subsequently used to calculate the parameters quantifying passive mass transfer in the PBPK model (Table 1). The pharmacokinetic behavior of CFDA SE was described by passive exchange as given in Equation 7. We considered an intravenous dose of per kg body mass [68] corresponding to an inflowing concentration of for a duration of for a mouse. Note that the concentration of the compound in the inflowing blood encompasses the corresponding equilibrium concentrations in the red blood cells and the plasma, respectively. The model for CFDA SE was in particular used as a proof of concept for the general performance of the spatially resolved model. We could show with this model that overall mass conservation is satisfied, see Table 1 in Text S1. Since metabolization of CFDA SE was not considered here, concentrations of CFDA SE in the in- and the outflow alone could be used for this essential step in model validation. The outflow curves for the spatially resolved model (Figure 6) show two effects, a temporal delay and a more smeared-out form of the peak from the spatially resolved simulation compared to the PBPK compartment simulation. The reasons for these observations become clearer when considering the temporal development of the concentrations in the four hepatic subspaces. The spatially resolved model no longer considers mean concentrations in well-stirred compartments but rather calculates heterogeneous distributions of these compounds. Likewise, the transition times needed to flow from the supplying to the draining vascular geometry are heterogeneous due to the different routes taken. We next visualized the total CFDA SE concentration in the HHS (Figure 7) obtained as the weighted average of the concentrations in the different subspaces, (11) Note that this is the quantity one observes in general for CT or MRI contrast agents by 3D imaging. In Figure 7 and in a Video S3, the different phases of the first pass of drug perfusion and distribution are shown. Also, the subsequent wash-out of the compound can be observed once the incoming pulse has passed through the liver. Notably, our spatially resolved model describes drug passage as a continuous process which may be used to complement experimental image data at discrete time points. Finally, we simulated steatotic cases where lipid accumulation in the cellular space of the liver influences the distribution behavior of compounds. In particular, we considered whether our spatially resolved simulations may be useful to support diagnostics and imaging. Concentration changes of CFDA SE due to spatially homogeneous and spatially heterogeneous states of steatosis are shown in Figure 8. As a pharmacokinetic application including intracellular metabolization, we next considered the distribution and metabolization of the sedative midazolam. For model establishment and parameter identification, we used previously published PK data [66] for mice obtaining an intravenous dose of per kg body weight. Metabolization of midazolam by CYP3A was quantified by using gene expression data as a proxy for tissue-specific protein abundance within a whole-body context [15]. This also involves a specific quantification of the hepatic metabolization capacity which is an essential prerequisite for the consecutive parametrization of mass transfer in the HHS. The PBPK model of midazolam was pre-parametrized with the physicochemical compound parameters molecular weight, fraction unbound and lipophilicity. Subsequently, the compound parameters as well as the metabolization parameters were fine-tuned with respect to the experimental PK data [66] (Table 1). The simulated plasma time curves obtained with the thus established PBPK model are in good agreement with the experimental PK data in mice (Figure 5). For the midazolam PBPK model in Figure 5, a concordance correlation coefficient was found, see also Figure 3 in Text S1. We next used the model parameters identified in the midazolam PBPK model for the spatially resolved model. As before, mass transfer of midazolam within the liver was described by passive exchange between the sinusoidal and interstitial subspace as well as the interstitial and cellular subspace as given in Equation 7. In addition, a nonlinear cellular metabolization according to Equation 8b was considered in this model. Values for the parameters in the equations are listed in Table 1. We considered a dose of per kg body mass, corresponding to an inflowing concentration of for a duration of. Outflow concentration time curves from the draining vascular system for the healthy state are shown in Figure 9. The total molar amounts (concentrations integrated over the whole liver) of compounds contained in the red blood cells, plasma, interstitial and cellular subspaces are plotted in Figure 9. In the simulations, we again observe a delayed and more smeared-out peak in the spatially resolved model. After simulated time, our model predicts a metabolization of approximately of the injected midazolam (healthy state), the rest having flown out from the model or still being present in the HHS and vascular systems. For midazolam metabolization, we also considered steatosis and -induced liver necrosis. In the homogeneous and heterogeneous steatotic state, an increase of the metabolization compared to the healthy state by and, respectively, can be observed, again after simulated time. For liver necrosis following intoxication [61] our simulation predicts a decrease of of the metabolized midazolam amount after. Finally we considered a model for the antibiotic spiramycin for which experimental data for an isolated liver were available in the literature [32]. For model establishment and parameter identification, we again used previously published PK data [67] for mice obtaining an oral dose of per kg body weight of spiramycin. Intravenous PK data are generally necessary for PBPK model development in order to identify systemic clearance capacity and distribution behavior without overlaying processes in the gastro-intestinal tract during oral absorption. Since intravenous PK data, however, were not available for mice, intravenous monkey PK data [69] were used for establishment of the fundamental model structure (Figure 1 in Text S1). We considered a linear metabolization term and pre-parametrized the distribution model with the physicochemical compound parameters (MW, ,). Based on the structure of the intravenous PBPK model, we then established a model for oral administration of spiramycin in mice [32]. Subsequently the model parameters were adjusted with respect to the experimental data [67] (Table 1). As before for midazolam, the spiramycin PBPK model provides a quantitative description of hepatic clearance capacity. The simulation time curves with the mouse PBPK model for intravenous spiramycin administration are in good agreement with the experimental plasma concentrations (Figure 5). For the PBPK model for spiramycin, we obtained a concordance correlation coefficient, see also Figure 3 in Text S1. Based on the validated mouse PBPK model for spiramycin we parametrized the spatially resolved model which is structurally identical to that of midazolam, except for the (now linear) metabolization kinetics. The spatially resolved model was then used to simulate experimental data for administration of spiramycin in an isolated liver [32]. The model structure of our spatially resolved model corresponds entirely to the experimental setup of the ex vivo assay, the availability of such highly specific data provided the opportunity to further validate our model. In the experiments [32], perfusion was performed using a buffer not containing red blood cells. The volume fractions from Equation 1 were hence changed to and. Moreover, a total perfusion of was used, which changes the flow velocities in our model and requires using a smaller time step (). Passive exchange between plasma, interstitial, and cellular subspaces was again modeled as in Equation 7, mass transfer involving red blood cells, however, was set to zero to take into account the specific experimental setup [32]. Due to the unphysiologically high flow rate, the local effective permeability parameters between interstitial and cellular space were adapted to and. An inflowing spiramycin concentration of for a duration of minutes was used as inflow condition reproducing the inflowing concentration profile in the experimental setup [32]. For a comparison to the experimental data reported in Figure 2 (wild-type) in [32], the outflowing rate of spiramycin was computed and plotted in Figure 10, again for the healthy state and the two steatotic states described above. Comparing experimental outflow concentrations and those simulated using the spatio-temporal model for the healthy reference case, a concordance correlation coefficient is obtained. Complementarily, volume renderings were generated at different time points after the end of the inflow for minutes (Figure 10) and show the spatial distribution of the spiramycin concentration immediately. This comparison illustrates very nicely how our spatially resolved model can be used to relate macroscopic observations in the plasma to distribution processes at the tissue scale.
We here present a spatially resolved model which describes the perfusion, distribution, and metabolization of compounds within the liver. The model structure is based on mass transfer equations obtained from PBPK modeling and vascular structures generated from micro-CT imaging. Our model excludes in particular any recirculation through the body such that metabolization and distribution of compounds can be considered without any overlaying effects. After the end of the initial administrations, a bi-phasic behavior can be observed which is initially governed by the distribution within the tissue and a slow release afterwards. Note that wash-out after the end of injection is additionally determined by advection in the blood flow. Comparing outflowing concentrations from our spatially resolved simulations to those from PBPK compartment models showed a temporal delay, both for CFDA SE and midazolam (Figures 6 and 9). This is because the compound now needs to pass sequentially through the supplying vascular system, the homogenized hepatic space and the draining vascular system. Different paths through the liver model require different transit times, hence the peaks are more smeared-out in the spatially resolved simulations. This is further emphasized by the temporal development of the concentrations in the four hepatic subspaces for the CFDA SE simulations (Figure 6). The spatially resolved model no longer considers mean concentrations in well-stirred compartments but rather calculates heterogeneous distributions of the concentrations. Likewise, the transition times needed for the compounds to flow from the supplying to the draining vascular systems are heterogeneous due to the different routes taken. This shows the general performance of the spatially resolved model where mass flows follows the physiological architecture of hepatic tissue governed both by vascular geometry and the composition of the connecting hepatic tissue. While this temporal delay only plays a role during first pass perfusion or similar sudden incidents, results from the spatially resolved model can nevertheless be used to revise PBPK model parameters by comparison with targeted experimental data [32]. Previous approaches already described macroscopic effects such as transit time distribution [10], [11], this can also be reproduced using our model. In addition, our approach provides a mechanistic interpretation and visualization of the underlying processes. Our model allows for example a physiology-based description of the liver, thus providing more insight into drug distribution and underlying clearance processes. Likewise, in contrast to fractal models [11] translating the vascular branching to effective pharmacokinetics parameters, we consider the actual anatomical geometry of the organ and its vascular structures. A highly resolved representation is indispensable for models that can also describe individual, potentially heterogeneous, pathologies of the liver. One major drawback, however, of the spatially resolved model is the highly increased computational effort required to run the simulations, see Table 1 in Text S1. To initially validate our spatially resolved model, we compared simulation results for spiramycin to experimental data obtained ex vivo with an isolated liver. The outflow concentrations simulated using the spatially resolved model and the experimental measurements in [32] are not in full agreement. Note, however, that the simulations of the isolated perfused liver are actually a prediction, since the original equations in the PBPK model were initially adjusted with respect to in vivo PK data [67]. In the light of this workflow it should be noted that the PBPK model represents only an intermediate step before the final spatially resolved model is ultimately established. It is only in this subsequent step that the liver model is integrated in the spatially resolved model, in this case to simulate ex vivo data from an isolated perfused liver [32]. Our approach hence extrapolates in vivo results obtained in a whole-body context to ex vivo data generated in an isolated liver as such supporting a structural transfer of knowledge. Hence, the setup of an isolated perfused liver is a suitable test case. The drawback of this prediction approach is the necessity of integrating experimental data coming from different sources which may partly explain the deviations in the stationary phase during the first 15 minutes during the onset of perfusion. While deviations between experimental data and simulated concentrations can be attributed to large experimental standard deviations or limitations of in silico to ex vivo transferability, a general agreement between the spatially resolved model and experimental data can be observed (Figure 6). In particular, the clearance rate after the interruption of the spiramycin inflow is in good agreement with the experimental data. This illustrates how our spatially resolved model can be used to relate macroscopic observations in the plasma to distribution processes at the lobulus scale. When applying the combined model to the case of steatosis it was found that already a spatially homogeneous change in the tissue composition leads to spatially heterogeneous differences in the distribution (Figure 8). The observed behavior showing the effect of an increased intracellular lipid content is actually a zonation effect on the length scale between terminal edges of the supplying and the draining vascular system. As discussed above, the qualitative result and the overall clearance are correct even though the flow distance between the two vascular trees is not the real hepatic lobule size. It was also found that the increased lipid content of the cells leads to longer intracellular retention times since the bound and thus immobile drug fraction increases. In turn, this leads to higher metabolization rates since lipid binding protects the compound from a fast wash-out due to increased retention times in steatotic livers. Differences between spatially homogeneous and heterogeneous steatotic states were also analyzed (Figure 8). It was found that the difference in lipid accumulation between different lobes and within the lobes had an observable influence on the concentrations as retention times in the cellular subspace are longer in case of higher lipid accumulation. This heterogeneous effect is only visible in spatially resolved modeling (Figure 9). The model thus provides a mechanistic description of pathophysiological states of the liver and can moreover distinguish between different spatial patterns of the pathology. Let us point out that both the temporal delay of the outflowing peak and differences between different steatotic states are inherent properties of the spatio-temporal model that the original compartmental PBPK model cannot describe. In contrast, the spatially resolved model can capture these effects in a qualitatively plausible way. Our model predicts increased metabolization in steatotic livers, but decreased metabolization following -induced liver necrosis. The simulations hence provide testable predictions which can be compared to previously published experimental data [70], [71]. For steatotic livers, drug lipophilicity has been related to intrinsic elimination clearance in rats with nonalcoholic steatohepatitis (NASH) and control rats, respectively [70]. From this study, an increased clearance of approximately can be estimated for midazolam () in steatotic animals. Even though this relationship has been established in rat livers perfused in situ and cannot be translated directly to our model, it nevertheless confirms qualitative validity of our simulations, since our model predicts an increased metabolization between and. For more detailed comparisons, simulations and experimental measurements would need to be performed for the same experimental setup and in particular in the same species. However, the significantly increased clearance found experimentally in steatotic animals already points to the necessity of a more refined diseased model of steatosis since lipid accumulation alone is obviously not sufficient to explain the observed decrease in metabolic capacity. Possible model extension include, amongst others, previously discussed changes in microcirculation [72] and intracompartmental permeability [70]. For -induced necrosis, our computational findings are also qualitatively validated by experimental observations, where a decreased metabolization of midazolam after pretreatment has been found in rats [71]. Comparing the experimental findings with our current model structure indicates that decreased cytochrome levels in -treated animals need to be considered as future model extensions. Here, our spatially resolved model could in particular be used to differentiate the contributions of enzymatic depletion and volumetric extension of necrosis on the decrease of metabolic capacity. Despite the performance of the newly developed spatially resolved model, several limitations need to be addressed, which represent excellent opportunities for future model refinement. On the technical side, a more detailed geometric vascular model and flow simulation [73], not only using constant velocity in each cylinder could be considered. However, all this will drastically increase computational costs with little benefit as the intravascular flow patterns are largely irrelevant for what happens in the HHS. Deformations of the organ as in [24], [74] could also be taken into account. Likewise, changes of the effective blood viscosity [75] besides those due to the Fhrus-Lindqvist effect [44] could also be considered in the model. The hepatic artery as the second supplying vascular system with other inflow concentration could become part of the model if its geometry and the local mixing of blood provided by portal vein and hepatic artery is known for the concrete situation considered [76]. This would also allow for more realistic flow velocities in the SVS. More generally, perfusion heterogeneity could also be considered as well as geometric scales of the perfusion [25]. A more detailed sensitivity analysis than merely one with respect to the vascular geometry (Figure 7 and Table 1 in Text S1) should be performed. For this purpose, known variations as well as measurement uncertainty of both PBPK model parameters and physiological/geometrical data need to be quantified, see e. g. [77]. For a physiologically relevant simulation output, such a sensitivity analysis will require substantial experimental and computational effort and should be part of a future study. Other implementations of the advection-PBPK simulation in the HHS should be investigated as well as the influence of computational resolution on the results. Comparing such fundamentally different implementations, however, is beyond the scope of this article. When considering other metabolization processes, additional compounds, e. g. products formed by the metabolization or compounds only stored in the cellular HHS subspace can easily be included in the model. The exchange across membranes, in Equation 7, can also be extended easily by active or other nonlinear processes. As discussed above, comparing our computational simulations of pathophysiological states of the liver to experimental data [70], [71] suggests several model extensions. For steatosis, these include, but are not limited to, a significant increase in liver weight as observed in [59, Table 6] as well as changes and spatial variations in the effective permeability in Equation 4 and the volume fractions, , and, as a significant decrease of functional capillary density (sinusoidal length per area) was reported [59, Figure 17]. Sinusoidal flow velocities, however, were not observed to change significantly [59, Figure 16]. Other studies indicate that a change in the microcirculation should be taken into account in a more realistic model of steatosis, see [72]. Moreover, changes of the intracompartemental permeability [70, Figure 4A] as well as the activity of drug metabolizing enzymes due to steatosis as discussed in [78] may affect the cellular metabolization of compounds. For -induced liver necrosis, changes in cytochrome levels [71] need to be considered in addition to necrotic changes in organ geometry. Here, our spatially-resolved model together with targeted liver histology could be used to differentiate between the different contributions to the decrease in metabolic capacity. Such integrative studies will allow further systematic analyses including iterative model testing and refinement in the future. More general pathological situations can be considered if one has solid knowledge of their spatial heterogeneity and their influence on the model parameters. In case of drugs being administered, also temporal changes of the parameters are possible and can be included in our model. A sensitivity analysis of the spatially resolved model with respect to such parameter perturbations could help to quantify their influence on the heterogeneity of drug distribution. The model in general is not specific for mice, so it can be applied to other species provided the geometry information and PBPK parameters are available. Possibly other connectivity patterns between larger vascular structures and sinusoids depending on the species [79] (or, closely related in the simulation, diffusion of compounds through vascular walls) need to be taken into account. The vascular tree geometries used in the model are easily exchanged if more detailed experimentally [80] or algorithmically [42] determined data is available. Similarly, more detailed information about the geometric shape of lobuli (as in [81] for human livers) could be taken into account. In particular, in vivo imaging with a slightly higher level of detail than used here will allow running simulations for patient-specific vascular geometries, thus providing great promises for imaging and diagnostics in the future. Corrosion casts [82], or other types of ex vivo specimens, also scanned in micro-CT, provide higher resolution as time and high radiation doses are not an issue, but obviously do not permit in vivo imaging. Even higher resolution could be obtained by extracting vascular geometries from optical microscopy images of histological serial sections [80]. This, however, requires a tremendous experimental and image processing effort and again is not applicable in vivo. As discussed above, possible zonation effects are qualitatively correctly observed at the length scale between the two incomplete vascular trees in our model rather than the actual length scale of hepatic lobuli. For correct observations in lobuli, our organ-scale simulations should be complemented by sinusoid-scale [83] or lobule-scale simulations in a multi-scale framework [21], [84]. Since the model can deal with pathological states of the liver and in particular spatially heterogeneous such states, their influence on the intrahepatic distribution of compounds could thereby be simulated pointing to future applications of spatio-temporal modeling in diagnostics. Here, comparison of our continuous simulations with new MRI or CT based image data could support the detection of pathological deviations. Predicting contrast agent distributions may help optimize time points for imaging after injection, benefiting from the much higher temporal and spatial resolution which our simulations can provide. The comparison of simulated and measured contrast agent distributions could therefore be used to identify changes in physiological parameters such that pathologies can be diagnosed. The possibility to simulate heterogeneous distributions provides also important applications for the prediction of toxic side effects. The spatially resolved model allows a location-specific prediction of exposure profiles within the liver. PBPK models have been linked before to models at the cellular scale to predict toxicity responses within hepatic metabolism in response to paracetamol [21]. Together with the spatially resolved model, this can now be used to simultaneously simulate intralobular exposure profiles and the specific cellular response. This allows an in silico prediction of toxic side-effects following the drug administration during the first pass perfusion. Simulations of spatial heterogeneity can also be used to describe local zonation effects within an whole-organ context. PBPK models have been used before to describe genotype-specific differences in hepatic drug uptake [19] and intracellular metabolization [20]. Since the corresponding equations are also used in the spatially resolved model, it also becomes possible to describe first pass effects in a genotype-specific way. Our spatially resolved model could be used for a wide range of technical and medical applications. It could for example be used for hypothermic machine perfusion [33] of livers to be transplanted for which mere static cold storage is ineffective. In this case, recirculation by a perfusion device needs to be considered, for which the influence on relevant compound concentrations can be described based on the existing PBPK models. Moreover, the model could be used to improve treatment planning for islet cell transplantations [34]. Here, mainly the perfusion simulation is needed to predict the distribution of a concentration of cells (not solutes) injected in the portal vein. Similarly, the model could help to improve intrahepatic injection of compounds, as discussed in the introduction. Another application could be optimization of targeted drug delivery [37] where drugs are injected in bound form and released at the desired location by mild hyperthermia induced by focused ultrasound. For this purpose, the model has to be combined with a heat transfer simulation [85]. For in vivo modeling within an organism context, more complicated full-body recirculation needs to be taken into account. This in turn requires our model to be integrated in whole-body simulations, regardless whether or not other organs are implemented at a comparable level of detail. Since our model has a spatially resolved internal state and (depending on the exchange and metabolization kinetics) may behave non-linearly, a transfer function approach [86] is not immediately applicable. Including recirculation in combination with our spatially resolved model will allow to mechanistically describe the distribution kinetics of fast acting drugs shortly after administration, similarly as it is has been done before with other circulatory models [11]. Extending such earlier approaches, our model will additionally use CT-based vascular trees within the liver. While the spiramycin simulations above show general agreement with experimental results in [32], this is just a first step towards an exhaustive validation of our approach. Starting points for the important step of model validation in future studies could be comparing simulated and experimentally measured outflow concentrations similar to what was done in [32] or time-resolved imaging of the distribution of tracers (at least imaged on some slices; see e. g. [87]) for comparison to simulation results as in Figure 7. For the latter, also mean transit times [87] estimated from the results in Figure 6 or from Equation 6 in Text S1 could be used for comparison to experimental results. The setting of a compound not entering the cellular subspace, such as in [88] for MRI contrast agent in rats, could be a starting point with a simpler model. In both cases, the ex vivo setting potentially allows for artificially low and thus slow total perfusion, possibly enhancing CT or MR imaging at multiple time points. Much higher spatial resolution at a single time point could be obtained from histological whole-slide scans for which registration [89] and analysis [90] techniques are available. More generally, validation combined with a parameter sensitivity analysis could also help to narrow down parameter ranges where the model predicts physiologically realistic behavior. In this regard, our model could be used for experimental planning to estimate the required spatial and temporal resolution for imaging. Likewise, the number of animal sacrifices could be minimized by specific design of experiments. The model could furthermore be used to quantify the contribution of first pass effects to the overall bioavailability and the experimental variability. As outlined above for steatosis and -induced liver necrosis, our model can be used in combination with targeted experimental data to iteratively investigate pathological changes in liver physiology. Validation or falsification of computational predictions can thereby support mechanistic insights in underlying processes such that overall model structure can be adjusted accordingly. Due to the large level of detail included in our model, such modifications can be directly assigned to specific pathophysiological changes. It is thus possible to test hypotheses about the behavior of pathological livers or to analyze pharmacokinetic effects such as zonation [65]. To this end, PK data, which are ideally sampled densely in time both in the portal vein and in the hepatic vein need to be compared to specific simulation results. Experimentally, one could for example use isolated, pathological livers from genetically modified mice strains or use PBPK models to correlate plasma PK data in these animals with exposure profiles in the liver. Verifying or falsifying these in silico results can then, in turn, trigger further model refinement. We here present a novel method for spatially resolved simulations of first pass perfusion in the liver based on mass balance equations from physiologically based pharmacokinetic modeling as well as vascular geometries obtained by in vivo imaging. The spatio-temporal description of blood flow through the vascular systems in combination with distribution models used in pharmacokinetic modeling allows a mechanistic yet local description of compound perfusion within the tissue. Our combined model is capable of representing spatial parameter heterogeneity, so that the local impact of pathophysiological changes within the liver can be analyzed. The model was used in the present study to investigate spatio-temporal effects of first pass perfusion for exemplary drugs. Two pathophysiological states, steatosis and -induced necrosis, were considered and were found to influence the distribution and metabolization of the compounds. Future applications of the model include optimized design of therapeutic treatments where spatially heterogeneous distributions or spatio-temporal perfusion effects are of relevance, e. g. targeted drug delivery, islet cell transplantations, or catheter placement for intrahepatic injections. We expect the spatially resolved model to be the foundation for further physiologically highly detailed modeling which will help to address specific spatial aspects of pharmacokinetics in the future.
The liver continuously removes xenobiotic compounds from the blood in the mammalian body. Most computational models represent the liver as composed of few well-stirred subcompartments so that a spatially resolved simulation of hepatic perfusion and compound distribution right after drug administration is currently not available. To mechanistically describe the local distribution of compounds in liver tissue during first pass perfusion, we here present a computational model which combines micro-CT based vascular structures with mass transfer descriptions used in physiologically based pharmacokinetic modeling. In the resulting spatio-temporal model, hepatic mass transfer is governed by the physiological architecture and the composition of the connecting hepatic tissue, such that hepatic heterogeneity and spatial distribution can be described mechanistically. The performance of our model is shown for exemplary compounds addressing key aspects of distribution and metabolization of drugs within a mouse liver. We furthermore investigate the impact of steatosis and carbon tetrachloride-induced liver necrosis. Notably, we find that our computational predictions are in qualitative agreement with previous experimental results in animal models. In the future, our spatially resolved model will be extended by including additional physiological information and by taking into account recirculation through the body.
medicine
computer science
mathematical computing
liver diseases
mathematics
materials science
gastroenterology and hepatology
drugs and devices
drug distribution
porous materials
computerized simulations
drug metabolism
material by attribute
pharmacokinetics
Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver