Spatio-Temporal Simulation of First Pass Drug Perfusion in the Liver

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.

Previously published PK data [1] for monkeys obtaining an intravenous dose of 50 mg spiramycin were used for basic model establishment. Intravenous PK data are necessary to identify systemic clearance capacity without overlaying processes in the gastro-intestinal tract during oral absorption. A linear hepatic metabolization term was considered in the model which was pre-parametrized using physicochemical compound parameters (MW, f u , log P ).  or the physiology of the organism (V max , K m ) is perturbed to show its impact on the simulation behavior. Experimental PK data [2] are shown (stars) such that the simulation results (lines) can be compared to experimental observations.

Implementation of Advection in the Vascular Trees
For the discretization and implementation of the vascular advection model, we use a 1D ELLAM scheme going back to [3]. The PDE in Equation 2 in the article is written in weak form using special test functions ϑ chosen to be constant along the characteristic curves of Equation 2 in the article.
An edge E discretized in space by equidistant grid points x j = x 0 + jh and in time by t k = t 0 + kτ satisfying the condition i.e. the velocity may be at most the width of one grid cell per time step. Note that it is not necessary to satisfy a CFL condition [4] here, ELLAM can handle larger time steps, but it will be convenient to limit as we will see below. Basis functions for spatially discretizing the molar concentration profile are standard 1D piecewise affine-linear FE basis functions ψ(x). The test functions ϑ are defined in terms of ψ as , forming a partition of unity.
For the time interval [t k , t k+1 ], the weak formulation also taking into account the cross-section area A = A ∅ (E) can be stated as follows Substituting the spatial discretization c(x, t k ) = j C k j ψ j (same for k + 1) and test functions ϑ k+1 , we integrate by parts and observe that the space-time integral and other terms vanish due to the choice of ϑ being constant along characteristic curves. For the mathematical details omitted here, we refer to [3, Section 3] (with diffusivity D = 0). The system of equations for computing C k+1 E m = E × {t m } and terms for j − 2 and j ± 1 only present if the indices are in the respective range. Inflow and outflow at the spatial boundary will be discussed below. This results in a system of equations of the form with M being a finite element mass matrix the cross-section area A ∅ ) and a quadridiagonal matrix T with one upper and two lower diagonals. In fact, we use a lumped mass matrix M for reasons discussed at the end of this section. Note that the factor A ∅ in Equation 2 does not make a difference when only considering a single edge as it appears on both sides of Equation 4, but will come in handy when we now consider branchings because it automatically takes care of mass conservation of the compounds considered. Moreover, let us point out that we always index grid points such that indices increase along the direction of flow.

Inflow Boundary Conditions.
Inflow boundary conditions (depending on time) are mapped to spatial data on a virtual extension of the corresponding edge by one grid cell in upflow direction, essentially resulting in an additional term on the right hand side of Equation 4. bifurcation in the bifurcation in the supplying vascular system draining vascular system Branching Points. As for branching points, let us consider the case of the supplying vascular system first. In this case, the end point of parent edge coincides with the start points of the daughter edges, only one degree of freedom (DOF) as part of the parent edge is assigned to that point. Basis functions ψ b corresponding to such DOF are supported in multiple edges. Test functions ϑ k+1 b (x, t k+1 ) at the new time step still coincide with the corresponding ψ b , tracing them back along characteristics results in discontinuous ϑ k+1 b (x, t k ) with modified slope on the upflow (i.e. parent) edge due to the discontinuous velocity, see Figure 4. Let us point out that, despite the sketch showing only a bifurcation, this scheme is easily adapted to multifurcations.
In case of the draining vascular system, the end points of daughter edges coincide with the start points of the parent edge. Basis functions ψ b are the same as before, test functions at the old time step coincide with the corresponding ψ b . However, tracing them back in upflow direction now results in continuous test function ϑ k+1 b (x, t k ) at the old time step, again with different slope on the upflow (this time daughter) edges, see Figure 4.
Grid spacings h for the different edges need not be the same (which is not even possible because edges have arbitrary length and are to be discretized with an integer number of grid points), but the time step is a common one. If we require at least four grid points per edge and Equation 1 for each edge, at most one bifurcation needs to be considered for adapting the ψ and ϑ. This is one reason for the restriction in Equation 1.
In the analog of Equation 2, we now need to integrate over multiple edges E and take into account the respective cross-section area A ∅ = A ∅ (E). When considering a whole vascular tree, the value vectors C k in Equation 4 then have block structure with one block per vascular edge. The corresponding analogs of matrices M and T in Equation 4 also have a block structure, mainly consisting of individual such M and T on the diagonal blocks with individual additional contributions from ψ and ϑ supported on two adjacent edges. For a more detailed discussion of these couplings, we refer to [5,Section 4.3].
Another observation in [5,Section 4.3] is that using Equation 4 may lead to artificial oscillations at draining bifurcations introducing over-and undershoots, i.e. not preserving the ranges of the molar concentrations. This well-known effect both ELLAM simulations can be remedied [6] by lumping M with the additional benefit of only needing to invert a diagonal matrix but the drawback of introducing artificial diffusion due to the numerical scheme.

Implementation of Advection in the Homogenized Hepatic Space
Discretization of Equation 5 in the article leads to an analog of Equation 4 with a 3D finite element mass matrix M and a matrix T on the right hand side with a more general sparsity structure. Again and for the same reasons as in the 1D case, lumped M was used.
Test functions ϑ for the ELLAM formulation are obtained from standard multilinear FE tent functions ψ as before by tracing back along the characteristics. This is more technical than in the 1D case because any upflow direction is possible now and quadrature in 3D needs to distinguish more possible cases of overlap of supports of ϑ and ψ, also resulting in a more complicated sparsity structure of the resulting matrix T . Similar to Equation 1 in the 1D case, we again restrict the maximal time step so that the maximal velocity is smaller than one grid cell per time step.
Given the discretization of the vascular trees using 1D piecewise linear basis functions and the HHS using piecewise trilinear basis functions, the source terms can be treated by matrix-vector multiplication. For this purpose, the line integrals over terminal edges of products of the two types of basis functions need to be evaluated. As the product is a fourth-order polynomial, we use Gauss quadrature with 3 points [7] providing fifth-order accuracy. Next to bifurcations where the value is stored for the parent edge, this would require interpolation between values stored for different edges. As a slight simplification and reduction of technicalities, we consider piecewise constant basis functions along the edges and assemble the matrix accordingly.
Discrete velocity profiles as the one computed here are not exactly conservative, i.e. have vanishing divergence away from sources and sinks, but this effect can be compensated for [6,Section 5]. Additionally, we are dealing with lower-dimensional sources and sinks, making it necessary to compensate for discretization artifacts. 0. As a preprocessing step, we start with an initial concentration of constant 1 in the HHS and the vascular systems with the same value 1 for the inflow to the supplying vascular system. Computing one coupled time step for the vascular system and HHS advection with appropriate source and sink terms followed by Gaussian diffusion by 1 / √ 2 voxels [8], we determine an additive correction for the HHS advection step.
1. In each simulation time step, we first perform the advection time step in the vascular systems and compute the amount supplied to and drained from the HHS.
2. Next we compute sources (according to the values in the supplying vascular system) and sinks (according to values in the HHS near the terminal edges of the draining vascular system) and perform the HHS advection time step, followed by Gaussian diffusion by 1 / √ 2 voxels. We then apply the correction term from step 0, locally scaled to fit the current molar concentration values. After a second diffusion step by the same amount as before, we compare the update of the contained amount in the HHS to the amount balance of step 1 and rescale accordingly to satisfy conservation of molar amounts.
3. Finally, we compute the PBPK time step in the HHS which may change molar amounts by design as it represents not only exchange but also metabolization. Here, we merely impose a lower limit of 0 for the molar concentrations in case the numerical scheme has computed slightly negative values.
In the literature, more advanced ELLAM schemes have been presented. Instead of separating advection and PBPK time steps, the PBPK could be linearized as in [9] or Runge-Kutta-based higher order time stepping can be used [10]. However, for more involved PBPK terms, it is not straight-forward to estimate the influence of linearization and a decoupling of the time stepping for the two processes permits a more efficient implementation. For a more detailed discussion of spatially varying velocities, we refer to [11], the latter also presenting a multiscale ELLAM approach, and [12] for a numerical upscaling procedure.
Unstructured meshes have also been used in ELLAM simulations [13]. This approach could help resolve the velocity variations in the HHS near the exchange with vascular edges, but at the well-known cost of a less efficient implementation than the structured hexahedral meshes we use. This trade-off could be investigated in a further study.
HHS v a s c u l a r s y s t e m

Performance of the Perfusion-Storage Model
Simulations using the CFDA SE parameters were also used to (a) investigate how the level of detail in the vascular structures influence the resulting outflow concentrations and (b) to assess the computational performance. The datasets used for this comparison are available in Data S1. Figure 7 shows that about 800 leaf nodes in the trees are sufficient to avoid artifacts due to insufficient resolution such as flow heterogeneities. Higher numbers of leaf nodes only change the outflow concentration profile by a minor additional delay, but at the cost of increasing computational effort, see Table 1.
Computational Performance. The simulations can be run on a standard desktop PC, CPU times for the CFDA SE simulations as described above are reported in Table 1 for a Intel Xeon X7550 2.0 GHz CPU. For 800 leaf nodes, the spatially resolved simulation involves the two coupled advection simulations (actually being the largest portion of the workload) in addition to 49 114 instances of the PBPK ordinary differential equation in the HHS. Compared to a four-subcompartment simulation without spatial resolution, the computation time increases only by a factor of 160 241. Simulations were only run single-threaded, first test show a moderate gain in performance by parallelizing the code using OpenMP. However, for the resolution described above, we are far from real-time simulations.
Typical Transit Times. As we are not using the full vascular trees, we need to discuss the influence of incompletely resolved vascular trees on the velocity profile and the time a given metabolite takes to pass through the HHS. For assessing this, let be the average distance between midpoints of terminal edges of the supplying vascular system to the closest midpoint of edges of the draining vascular system, a typical distance the blood needs to pass through the HHS. Moreover, let v avg := 1 Volume(HHS) HHS v be the average Euclidean norm of the velocity, so that we can define a typical transit time for the blood through the HHS. We refer to [14] for an overview of other notions of transit times for metabolites in the liver. The amount of blood flowing between the vascular structures per unit time is independent of the level of detail in the vascular structures, merely the flow through each terminal edge changes. If the perfusion of the liver volume is essentially homogeneous, we obtain t typ ≈ liver volume · porosity total perfusion = Volume(HHS) · ϕ sin Q tot independent of the concrete vascular geometries. In particular, it is independent of the level of detail in the vascular structures, provided the perfusion is sufficiently homogeneous. With Volume(HHS) = 1.077 ml, ϕ sin = 0.104, and Q tot = 2.1 ml/min as above, we obtain t perf = 3.200 s. This estimate indeed matches the values t typ obtained for different vascular resolutions in Table 1.
Let us point out that the time for passing through the vascular structures of course depends on how much detail is contained in them. Thus the total transit time of a concentration depends on the vascular resolution. However, no metabolization takes place inside the vascular structures in our model, so only t typ or t perf matter for the total metabolization.

Software Architecture
In this Section, we give an overview of our implementation of the simulations. Using C++ notation, we present the main structure of the class structure, omitting most of the technical details. In particular, all the technicalities implementing ELLAM for the 1D and 3D advection are left out. The 1D case is described in [5,Section 4], our implementation here is based on the corresponding code available as part of [15].