DigiLoCS: A leap forward in predictive organ-on-chip simulations

Abstract

Digital twins, driven by data and mathematical modelling, have emerged as powerful tools for simulating complex biological systems. In this work, we focus on modelling the clearance on a liver-on-chip as a digital twin that closely mimics the clearance functionality of the human liver. Our approach involves the creation of a compartmental physiological model of the liver using ordinary differential equations (ODEs) to estimate pharmacokinetic (PK) parameters related to on-chip liver clearance. The objectives of this study were twofold: first, to predict human clearance values, and second, to propose a framework for bridging the gap between in vitro findings and their clinical relevance. The methodology integrated quantitative Organ-on-Chip (OoC) and cell-based assay analyses of drug depletion kinetics and is further enhanced by incorporating an OoC-digital twin model to simulate drug depletion kinetics in humans. The in vitro liver clearance for 32 drugs was predicted using a digital-twin model of the liver-on-chip and in vitro to in vivo extrapolation (IVIVE) was assessed using time series PK data. Three ODEs in the model define the drug concentrations in media, interstitium and intracellular compartments based on biological, hardware, and physicochemical information. A key issue in determining liver clearance appears to be the insufficient drug concentration within the intracellular compartment. The digital twin establishes a connection between the hardware chip structure and an advanced mapping of the underlying biology, specifically focusing on the intracellular compartment. Our modelling offers the following benefits: i) better prediction of intrinsic liver clearance of drugs compared to the conventional model and ii)explainability of behaviour based on physiological parameters. Finally, we illustrate the clinical significance of this approach by applying the findings to humans, utilising propranolol as a proof-of-concept example. This study stands out as the biggest cross-organ-on-chip platform investigation to date, systematically analysing and predicting human clearance values using data obtained from various in vitro liver-on-chip systems. Accurate prediction of in vivo clearance from in vitro data is important as inadequate understanding of the clearance of a compound can lead to unexpected and undesirable outcomes in clinical trials, ranging from underdosing to toxicity. Physiologically based pharmacokinetic (PBPK) model estimation of liver clearance is explored. The aim is to develop digital twins capable of determining better predictions of clinical outcomes, ultimately reducing the time, cost, and patient burden associated with drug development. Various hepatic in vitro systems are compared and their effectiveness for predicting human clearance is investigated. The developed tool, DigiLoCs, focuses explicitly on accurately describing complex biological processes within liver-chip systems. ODE-constrained optimisation is applied to estimate the clearance of compounds. DigiLoCs enable differentiation between active biological processes (metabolism) and passive processes (permeability and partitioning) by incorporating detailed information on compound-specific characteristics and hardware-specific data. These findings signify a significant stride towards more accurate and efficient drug development methodologies.

1 Introduction

The drug testing dilemma presents a significant challenge in pharmaceutical development, marked by high costs and a distressing attrition rate in accurately predicting human responses [1, 2]. A pivotal element in preclinical drug development is the accurate estimation of the first-in-human dose and different dosing regimens to keep drug levels within a therapeutic range. This demands precise assessments of hepatic clearance of the drug and human pharmacokinetics [3, 4]. Typically, the gold standard in drug development is the use of simpler in vitro systems to study drug metabolism, including liver microsomes [5] and suspension or plated hepatocytes [6]. The drug depletion data (time-concentration profile) are then analysed to determine the in vitro drug clearance rate. A simple mathematical model has been employed in earlier work, and it considers the in vitro system as a single compartment, the one-compartment PK model [7]. Well-mixing and instantaneous drug distribution is assumed and all biological processes, e.g., permeability and partitioning from cell culture media into intracellular milieu are lumped into drug clearance. This approach also cannot differentiate between compounds actually being metabolised and compounds bound to media proteins or hardware. The so determined in vitro clearance value is then extrapolated to humans (in vitro—in vivo extrapolation) and integrated into human physiologically-based pharmacokinetic (PBPK) models [4, 8] to predict human pharmacokinetics (absorption, distribution, metabolism and excretion (ADME)), before actually testing a new compound in humans. Although this approach is well-established in drug development and easy to use, it also systematically underpredicts human PK [9] by 5–10 fold across studies and compounds.

Microphysiological systems (MPS) and organ-on-chips, as well as 3D organoids, hold great promise to address more complex in vitro ADME, toxicology and pharmacology questions offering miniature, biomimetic systems that replicate key aspects of human organ physiology [10–12]. These technologies create an environment where human cells can grow and interact in an organ-specific context, providing insights into human biology and disease that were previously unattainable in conventional in vitro models or animal studies. OoC and MPS-based systems are not only used in today’s drug development for PK but also for assessing drug efficacy and toxicity [10, 13–15]. While the emulated in vitro biology of MPS and OoCs is getting ever more complex and produces more human-relevant data, these systems still fall short in considerably improving the prediction power of in-human situations, like PK [13, 16]. However, MPS and OoC data are also still analysed using the conventional mathematical model (one-compartment) that does not account for advanced biology. It remains unclear whether the OoC and MPS biology is still not human-relevant enough (and thus not producing human-relevant data) or the conventional mathematical analysis is the cause for the underprediction. Potentially, a digital twin framework that enables the mapping of complex on-chip biology to advanced mathematical models could provide a useful approach to enable OoC and MPS translation to humans and increase the prediction power, but it is currently lacking.

The current study aims to develop a digital twin approach integrating MPS and OoC data within advanced computational models of biology to improve the prediction of clinical clearances. DigiLoCs, our developed digital liver-on-chip simulator, facilitates the accurate description of on-chip complex biology to disentangle biological processes, namely clearance, permeability, and partitioning. The tool comprises and utilises information on complex biological processes (clearance, permeability, partitioning), hardware-specific information from the studied in vitro system, and compound-specific information. By accounting for more multi-dimensional information, the tool enables differentiation between active biological processes, such as metabolism, and passive ones, such as permeability and partitioning of a compound from cell culture media into the cellular environment. This tool offers a significant advancement over conventional approaches, which fail to explicitly consider passive biological processes and conflate them into a singular clearance mechanism. By providing a more detailed understanding of biological processes, our tool has the potential to reveal better insights into liver-chip biology. Drug depletion kinetics of 32 compounds were taken from literature covering commercially available liver-on-chips (CnBio [13, 17], Javelin), and 3D spheroids [18, 19], including fast and slow-cleared compounds. According to these studies, DigiLoCs outperform the conventional prediction approach considerably. The impact of a more accurate description of clinical clearance values on predicting human PK was investigated in a proof-of-concept study using propranolol. The kinetics of propranolol was predicted in humans using the conventional approach, DigiLoCs, and literature approach. The results obtained from DigiLoCs for propranolol in the proof-of-concept study were much closer to the actual observed human values than those of other approaches.

To the best knowledge of the authors, this is the first and biggest study so far, comparing head-to-head the performance of different hepatic in vitro systems to predict human clearance and demonstrating the impact OoC and MPS systems can have in the drug development process, enhanced through the modelling and prediction features of DigiLoCs.

2 Methods

In this section, the following are described: i) data used in the study for predicting human clearance, ii) DigiLocs, digital twin for liver-on-chip, iii) mathematical model, parameter estimation and sensitivity analysis for DigiLocs, and iv) translation to humans and prediction of human pharmacokinetics.

2.1 Data

In this work, published data on pharmacokinetics (metabolism) or toxicology studies of 32 drugs are used (See Table 1) to predict human pharmacokinetics.

[Figure omitted. See PDF.]

2.2 Mathematical model

A typical single-compartment model is described as follows. Let C(t) be the drug concentration in the chip at time t, V the volume of the chip, and CL_c the clearance parameter. Then we havewith initial value C(t = 0) = C₀. The solution to this ODE is(1)

Taking the logarithm on both sides,(2)which is equivalent to regression on log-transformed kinetic data.

A digital twin of liver-on-chip with three compartments that incorporate much more information on parameters related to both on-chip characteristics and drug-specific properties was developed. The three-compartment model considering media, interstitium and intracellular compartment is described as follows: Let C_m(t), C_i(t), C_c(t) and V_m, V_i, V_c be the concentration of the drug at time t, and volume of the media, interstitium, and intracellular compartment respectively. CL_c is the clearance parameter. Then we have(3) (4) (5)

The parameters here are defined as follows,(6)where K_{int_med} is the partition coefficient for the transfer of drug between the media and interstitium, P_endothelial is the permeability coefficient of the drug between the endothelial layer, SA_{med_int_liver} is the surface area of the interstitium, K_{water_cell} is the partition coefficient for water exchange or movement within the interstitium, K_{water_cell} is the partition coefficient for water exchange or movement within the intracellular compartment, PA_{cell_int} is the product of the surface area and permeability coefficient of the cellular membrane in the intracellular compartment, PA_{int_cell} is the product of the surface area and permeability coefficient of the cellular membrane in the interstitium compartment and F_unbound is the fraction unbound (media, reference value). All the relevant variables and parameters with descriptions are given in Table 2.

[Figure omitted. See PDF.]

We can write this system of linear ODEs in matrix form,(7)

The general solution for the system of ODEs at time t is:(8)where X is the matrix of eigenvectors of A, Λ is the diagonal matrix with λ₁, λ₂, λ₃ (eigenvalues of A) as diagonal entries, and C₀ is the initial value of variables at time 0.

The objective function for optimisation is as follows:where is the observed data point and is the computed value at time i respectively and n is the number of observed data points.

2.3 DigiLoCs: Digital twins for cell-based liver assays

DigiLoCs is a software tool, developed within this work that describes the on-chip complex biology more accurately in the context of use to predict clinical clearance values. The software comprises (Fig 1):

* modelling of complex biological processes (clearance, permeability, partitioning),

* hardware-specific information from the studied in vitro system and

* compound-specific information.

[Figure omitted. See PDF.]

Contrasting conventional approach, the DT approach uses biological, hardware, and physicochemical information to map the biological processes on liver-chip more accurately to in silico, thereby maximising the information leveraged. This results in the disentanglement of active (metabolism) and passive (permeability, partitioning) processes. Created with biorender.com.

The tool differentiates between active biological processes, such as metabolism, and passive ones, like permeability and partitioning of a compound from cell culture media into the cellular milieu. This contrasts with conventional approaches, where passive biological processes are not considered especially and lumped together into a single process, i.e., clearance.

Liver-on-chip technology provides a more physiologically relevant environment compared to traditional cell cultures or animal models, enhancing the simulation of drug responses using mathematical models. Hence, a more accurate mathematical description is needed. The three primary compartments of the liver chip considered in the model are media, interstitium, and intracellular space, which serve as dynamic environments where drugs are distributed, metabolised, and interact with hepatic cells. This compartmentalisation is based on concepts applied in human whole-body PBPK modelling.

The software tool is developed in the open-source programming language R [21] and seamlessly communicates with PK-Sim (https://www.open-systems-pharmacology.org/) via in-house developed functions. For more information, see esqlabsR package (https://github.com/esqLABS/esqlabsR). All analysis and plotting were also done in R. The proposed workflow does not interfere with existing wet lab Standard Operating Procedures (SOPs) for performing biological experiments and does not add an extra considerable burden to the user. DigiLoCs uses existing biological data, and its performance may be improved by measuring cell-associated compound concentrations in addition to the compound media depletion time course, which would add a minor extra step in the lab SOP. This, however, is negligible given the improvement in performance power and the confidence in the prediction.

2.3.1 Implementation of hardware specifications.

DigiLoCs map the chip architecture to a compartmental model to describe the time-dependent distribution of a compound on-chip. The compartment models use time-dependent ordinary differential equations (ODEs) and assume well-mixing within compartments. These are generally accepted to describe the distribution of exogenous and endogenous compounds and molecules.

A physical chamber separated by a membrane or connected by flow to another chamber is represented by a compartment in the software. Serial compartments are connected via media flow rates (typically in ml/min) between the compartments describing mass transport via gradients and diffusion and normalised by the volume of the originating compartment.

2.3.2 Implementation of biological specifications.

The biology (more precisely, the cell type exerting the biological function under investigation; here: metabolism) is mapped by two additional compartments representing the interstitial and intracellular space of the investigated biology (Fig 1). Transport rates from the cell culture media into the interstitial and intracellular milieu are described by two core processes:

* passive diffusion (driven by concentration gradients between compartments)

* thermodynamic equilibrium (distribution of compound between two phases)

These processes are then described by two main parameters in the computational model: permeability (how fast is a compound taken up?) and partitioning (how much of the compound is taken up by cells?). Lastly, the metabolism rate is reallocated to the intracellular compartment and corrected by the unbound fraction of the compound in the intracellular compartment. Similarly, in pharmacokinetics, the distribution of compounds is often described using a compartmental framework, which involves dividing a system into distinct compartments and modelling the processes that govern the movement of compounds between them.

2.3.3 Implementation of compound-specific information.

The following physicochemical properties of the investigated compounds are used in the software: i) lipophilicity (logP), ii) molecular weight (MW) and iii)fraction unbound (fu); to calculate up to six dependent downstream parameters (listed below). These parameters describe the partitioning from the main media compartment into the interstitial space and between the water fraction and both interstitial and intracellular space. Additionally, permeability across the endothelial barrier and between interstitial and intracellular spaces are calculated.

1. Partitioning: i) K_{int_med} ii) K_{water_cell} iii) K_{water_cell}

2. Permeability: i) P_endothelial ii) PA_{cell_int} iii) PA_{int_cell}

Here int refers to interstitial space, med refers to media, water refers to water exchange fraction, and cell refers to intracellular space. These parameters are calculated based on well-established and documented equations implemented in PK-Sim software [3], which is a comprehensive software tool for whole-body PBPK modelling. It enables rapid access to all relevant anatomical and physiological parameters for humans and common laboratory animals contained in the integrated database for model building and parameterisation. The same partition coefficient calculation methods as implemented in PK-Sim are also readily available and can be investigated:

* PK-Sim standard

* Poulin and Theil

* Rodgers and Rowland

* Schmitt

* Berezhkovskiy

Further, only the unbound fraction of a compound can be taken up by cells and be metabolised by cells. The unbound fraction in the cell culture media is typically informed by biological experiments. However, the intracellular unbound fraction is not often available or measured. Thus, two established QSAR models (quantitative structure-activity relationship) are implemented in the software to predict the unbound intracellular fraction of the investigated compound as a function of its physicochemical properties [22, 23].

2.4 Parameter estimation

Parameter estimation aims to find unknown parameters in a computational model and is estimated using experimental data collected from well-defined and standard conditions. By minimising the distance of theoretical function values and experimentally known data, the set of parameters in the model can be estimated. The parameters which are not directly measurable can be estimated using least squares or any other fitting methods to analyse the model quantitatively. Nominal parameter values are obtained from PK-Sim software.

Parameter estimation in DigiLoCs is a two-step process. Firstly, a customised cost function is implemented. This cost function calculates the weighted difference (ssq) between the model simulation (pred) from a specific compartment and the corresponding observed data (obs) for each time point according to the equation(9)

Common parameter estimation methods include maximum likelihood estimation and Nelder Mead optimisation. Nelder Mead, a non-linear optimisation method, is used to find the minima of the objective function in this work. Additionally, the partition coefficient between the intracellular (IC) and the main media compartment is estimated using the area under the simulated time-concentration profile of the IC and interstitial (IST) compartment and corrected for by the QSAR-predicted cellular unbound fraction (fu_cell) and the unbound fraction in the media (measured, fu_media):(10)

Kp_uu,obs is calculated from literature [24], where media and intracellular concentrations in hepatocytes were determined. Initially investigated for suspension hepatocytes in a 2D setting, the authors provide a scaling factor (∼4.9) to apply to human hepatocytes. Further, the ionisation state of the investigated compound (-1, 0, 1) results in a different partitioning. Otherwise, a range of possible partition coefficients are investigated. This is an additional anchor point for estimating the cost function value and links the simulated intracellular and main media compartment concentrations. Eventually, both differences are squared and summed up, resulting in the final sum of residuals. Based on this, a compound-specific scaling factor (SF) is calculated and used to scale the predicted human clearance:(11)

Specifically for the liver use cases, on-chip liver clearance and surface area between the main media compartment and the cell layer are estimated. It is possible to estimate other parameters, such as pre-calculated permeability or partition coefficient values.

2.4.1 Implementation of software.

Methodologically, DigiLoCs is implemented in the open-source programming environment R with its own package. A library of two common chip architectures and two cell types with six different chip-specific settings have already been implemented.

1. One chamber, no media flow

2. Two chamber, recirculating flow

3. Organ-on-chip (hepatocytes)

1. a. CnBio

2. b. Hurel 1

3. c. Hurel 2

4. d. Dynamic42

5. e. Javelin

The most straightforward system is a single, perfused microfluidic chamber containing one kind of cultured cell (e.g., hepatocytes) that exhibits functions of one tissue type linked to channels for fluid transport. In more complex designs, two or more microchannels are connected by porous membranes, lined on opposite sides by different cell types, to recreate interfaces between different tissues. The CnBio and Javelin chip settings are shown in Fig 2. For detailed information on the processes and modelling, readers are directed to Bhatia et al. [25]. These building blocks can be interchangeably used and connected, similar to the building blocks in PK-Sim. While the package provides a step-by-step guide to generate and run a simulation, the code communicates seamlessly with a generic PK-Sim model to determine partitioning and permeability values as described above, which are used in the simulation.

[Figure omitted. See PDF.]

Q_mix is the mixing flow rate in mL/min. Created with biorender.com.

2.5 Sensitivity analysis

Both local and global sensitivity analyses are used to quantify the impact of input parameters on the output variables. This involves varying certain input parameters and observing changes in the output variable intracellular concentration. The local sensitivity is estimated by changing one input parameter at a time while other parameters are held constant. The study provides insights into the sensitivity of various parameters and how they affect the output of the model.

The input parameters K_{int_med}, P_endothelial, SA_{med_int_liver}, K_{water_int}, K_{water_cell}, PA_{cell_int}, PA_{int_cell}, fu and CL_c are varied to evaluate the sensitivity of the output variable intracellular concentration (C_c). First, the local sensitivity is estimated by changing one input parameter by 10% at a time while the other parameters are held constant, and the changes in output variable C_c are compared. Eq 12 shows the local sensitivity index for C_c with respect to the varying model parameter (P_i), which is approximated by a small perturbation ΔP_i,(12)where C_c(P) is the model prediction of the intracellular concentration for parameter set P. The local sensitivity index is normalised to eliminate the effect of units,(13)

Global sensitivity analysis evaluates the effect of potential interactions of the input parameters in an output variable. The Sobol sensitivity analysis of the SALib package in Python is used to perform the global sensitivity analysis. Input parameters are sampled using the Saltelli sampler. The lower and upper bound of the parameters are set as 0.1-fold and 10-fold of the baseline parameter values, respectively. The first-order and total-order indices are estimated using the Sobol sensitivity analysis.

2.6 Translation to humans

Drug-related parameters extracted from OoC or any other in vitro studies can be scaled to predict clinical parameters using in vitro-in vivo translation (IVIVT) [13, 17]. The typical value of unbound intrinsic clearance CL_int(u) determined for each drug from the pharmacokinetic analysis of the in vitro depletion data is scaled up to a human liver equivalent unbound intrinsic clearance CL_int(u),(H) using(14)where HC is the human hepatocellularity of 120 million cells / g of liver, LW is the average human liver weight of 25.7g / kg of body weight [17] and fu_inc is the unbound fraction of drug in the incubation medium. The hepatic clearance (referring to whole blood concentrations) is then predicted (CL_H,pred) using the Well-Stirred (WS) model:(15)where Q_H is the average hepatic blood flow of 20.7 mL/min/kg of body weight and fu_b is the fraction of the drug unbound in blood. The fraction unbound in the blood was calculated for each compound from the known fraction unbound in the media (fu_p) and blood-to-plasma ratio (Rbp) according to the equation fu_p = fu_p/Rbp or directly used, if available from the literature. The predicted hepatic clearance (CL_H,pred) values were then compared to observed hepatic clearance (CL_H,obs) values (referring to whole blood concentrations). Following, the ratio of predicted (either via conventional or digital twin approach) and observed clinical clearance values for all investigated compounds was calculated. Using these ratios, a density distribution function was computed (function geom_density from ggplot [26]) for visual purposes only.

2.7 Prediction of human pharmacokinetics

Initially, a PBPK model is developed using qualified installations of the software PK-Sim that ensures the software has been installed correctly and thoroughly validated according to established procedures. A whole-body PBPK model includes an explicit representation of the organs most relevant to the uptake, distribution, excretion, and metabolism of the drug. These typically include the heart, lungs, brain, stomach, spleen, pancreas, intestine, liver, kidney, gonads, thymus, adipose tissue, muscles, bones, and skin. More information can be found in S1 File.

The tissues are interconnected by arterial and venous blood compartments, and each is characterised by an associated blood flow rate, volume, tissue partition coefficient, and permeability. If applicable, R (Distribution 4.0) and RStudio (Version 1.2.5) are used in the analysis for preprocessing and post-processing of data and model outputs [27]. The analytical approach is based on the principles set out in the EMA, FDA, and/or OECD guidelines for reporting on PBK M&S [11]. The developed PBPK model is used to describe the human kinetics of propranolol. Key kinetic parameters are informed by either clinical data, literature values or on-chip predictions. The translational workflow that integrates organ-on-chip results to predict human pharmacokinetics is shown in Fig 3.

[Figure omitted. See PDF.]

The digital twins of the humanised organ-on-chip systems, together with chip-specific information and physicochemical information, are developed in R. Created with biorender.com.

3 Results

We report the following results in this section. First, results from a simulation model based on a digital twin for selected compounds; second, sensitivity analyses; third, prediction of human clearance values; and finally, translation to human PK using propranolol as a proof-of-concept study. The liver clearance and surface area of the chip are estimated after fitting the drug kinetic data. The Poulin and Theil method of partition coefficient calculation was used due to its superior fit to observed drug kinetics, which outperformed alternative methods.

3.1 Simulating on-chip compound depletion

The digital twins for the investigated in vitro liver systems were successfully implemented in R and used to simulate the on-chip kinetics. After parameter estimation, the resulting model simulations describing the observed compound depletion data were visually inspected. The final parameter values can be found in Tables A and B in S1 File. Additionally, the squared sum of residuals was evaluated and deemed acceptable as it was less than < 0.01, which was the case for all simulations (data not shown). An example of on-chip kinetics is presented in Fig 4.

[Figure omitted. See PDF.]

Observed data are shown in blue dots, where the data for diclofenac, midazolam, and oxazepam are from Docci et al. [13], while propranolol is from Tsamandouras et al. [17]. The red, violet and green curve plots the drug concentration in the intracellular, medium and interstitium compartments, respectively: IC = intracellular, Ist = interstitium.

As can be seen, the digital twin approach (violet line) captures the on-chip kinetics (blue dots) very well. Simultaneously, the intracellular (IC) kinetics are plotted (red lines), clearly highlighting the difference in compound uptake and, thus, clearance rates. The remaining figures are presented in S1 File (See Figs A-D).

3.2 Sensitivity analysis

The sensitivity analysis, both local and global, was conducted to quantify the sensitivity of model output intracellular concentration with input parameters. The analyses were performed for various parameters, and the results indicated that the output is more sensitive to parameters such as the permeability coefficient of the endothelial layer, surface area of the liver sinusoids, and clearance. These parameters were estimated or calculated from clinical data/ experimental results. Clearance (CL_c) is identified as the most sensitive parameter with respect to intracellular concentration. The results imply that accurate values of these sensitive parameters are crucial for the model’s accuracy.

The normalized local sensitivity indices (Fig 5a) and the first-order and total-order global sensitivity indices (Fig 5b) for intracellular concentration across the input parameter set are shown. The results from both local and global sensitivity analyses show that the output is more sensitive to parameters P_endothelial, SA_{pls_int_liver}, fu and CL_c. SA_{pls_int_liver} and CL_c were estimated and nominal values were used for all other parameters. SA results imply that we need correct values of the constants P_endothelial, fu as they are more sensitive.

[Figure omitted. See PDF.]

(a) Blue bars indicate that the output and the input changes in the same direction, and the red bar indicates that the output decreases when the input increases. (b) The blue and orange bars represent first-order and total-order indices, respectively.

3.3 Predicting human clearance

The on-chip estimated clearance values were translated to total human clearance according to Eq 15. More detailed information is available in S1 File (see Tables A and B). Likewise, from the investigated studies (Table 1), in vitro unbound clearance values were available and scaled to human equivalents.

The conventional approach is a single-compartment model explained in Section 2.2 using a single ODE. The ratio of clinical observed human clearance values and either predicted human clearances using conventional mathematical modelling or the digital twin approach were estimated and converted into a density function for easier graphical visualisation as described in Section 2.6. As can be seen in Fig 6, the digital twin approach (DigiLoCs) outperforms the conventional approach considerably. The centre of the distribution is around 1, indicating a non-biased prediction of clinical clearance values, while the width of the distribution is very small. Quantitatively, the ratio for the digital twin approach overall compounds is 1.04 ± 0.31, with a coefficient of variation of 30%. In contrast, the conventional approach (red curve) majorly under-predicts the clearance values while maintaining a broad distribution and thereby adding to uncertainty in the prediction (0.56 ± 0.44, CV = 79.3%). The correlation plot between the observed and herein predicted clinical clearance values highlights the improved prediction performance of the DigiLoCs approach on an individual drug level. As can be seen in Fig 7 (example graph for CnBio Liver-on-Chip data), most of the compounds fall within the 1.5-fold line (Average fold error, AFE = 0.965). Similar correlation plots are presented in the S1 File (See Figs E-G) for the other in vitro systems.

[Figure omitted. See PDF.]

In total, a set of 32 compounds across three different in vitro liver systems have been investigated. The x-axis presents the ratio of predicted/observed clinical clearance values using either the DigiLoCs or the conventional approach, and the y-axis shows the frequency of the ratio observed.

[Figure omitted. See PDF.]

The solid line shows the line of unity, while the dotted line is 1.5-fold, and the dashed line has a 3-fold deviation.

3.4 Translation to human PK

The impact of accurately predicting human clearance values based on in vitro cell-based assays on predicting human PK was assessed using propranolol as a proof-of-concept case study. First, a human PBPK model describing the human kinetics of propranolol was implemented in PK-Sim and qualified with clinical observations. Next, the predicted human clearance value using either the conventional modelling approach or based on the same on-chip kinetic data was implemented in the human PBPK model simulating the kinetics after a single oral dose. Further, a population of n = 1000 patients was simulated to account for inter-patient variability. As shown in Fig 5, the implemented human PBPK model describes observed clinical data well (using clinical clearance values). When substituting only the clearance value with the conventional or the digital twin-based values, the impact on simulating human PK becomes apparent, while the conventional approach would overpredict (i.e., the on-chip clearance is underpredicted) the human PK (3-fold C_max, up to 6-fold overprediction of AUC). Moreover, this approach would actually simulate non-negligible concentrations of propranolol left over after 24 h. For repeated daily dosing, this would result in accumulation of propranolol in this hypothetical setting, which would have immediate implications for potential toxicity or efficacy considerations. On the other hand, the digital-twin based approach still slightly overpredicts the AUC and C_max, however only by 1.5-fold and captures the terminal phase correctly.

4 Discussion

The aim of this work is to improve the current prediction of human clearance values and to present a framework for translating in vitro findings to relevant clinical situations. The presented integrated translational approach combined quantitative OoC and cell-based assay compound depletion kinetics with an OoC-digital twin to simulate drug kinetics in humans.

Initial investigations revealed the potential to describe clinical clearance values more appropriately than is currently possible with the conventional approach. This simpler approach lumps biological processes together into a single process—clearance—and uses only minimal information available, e.g., only the cell number and media volume. While biological systems have evolved rapidly in the last decade, especially in the field of organ-on-chip and microphysiological systems, the applied mathematical models to analyse the quantitative complex biological data have been the same for decades (the early concept of clearance was introduced by Möllers in 1928, while Well-stirred model was introduced in 1971).

In contrast, the digital twin approach for the organ-on-chip and 3D spheroids comprises three building blocks described here: biological, hardware, and physicochemical information. The distinction between active and passive processes is achieved by an explicit description of uptake, distribution, and metabolism involved in the biological processes. Further, the digital twin links the architecture of the hardware chip with an advanced mapping of the underlying biology (intracellular compartment). The on-chip kinetics for 32 compounds (six compounds were removed from the initial set due to missing information) was well described, highlighting the drug-specific effects on cellular uptake and hence metabolism. Note that this analysis used the same biological information as used in the conventional approach and that no additional biological experiments were needed or performed to improve the outcome of the digital twin approach.

The predictive power of organ-on-chip and 3D spheroids over conventional approaches was revealed when the depletion data was analysed with the digital twins (Fig 6). Not only was the systematic underprediction issue resolved, but the uncertainty in prediction was also reduced by a factor of 3 (comparing CVs).

Lastly, we aimed to demonstrate the clinical impact of this approach by translating the results to humans using propranolol as a proof-of-concept example. Here, the head-to-head comparison clearly demonstrated the superior power of both quantitative biological data from OoCs and digital twins over conventional approaches in predicting human PK more appropriately (Fig 8). Although only one compound was used to demonstrate clinical impact, the workflow and process are easily applicable to other compounds. To the best of our knowledge, this study is the biggest comprehensive report to systematically assess the predictive power of organ-on-chip in the context of use of liver clearance.

[Figure omitted. See PDF.]

Pink dots are clinical observations (digitised from Borgström et al. [29]), while the blue solid line represents the mean of the patient population using the clinical observed clearance value. When using the clearance value (red line) based on conventional approaches, the area under the curve is 6-fold overpredicted. In contrast, using the digital twin-based clearance, the AUC is only 1.5-fold overpredicted, also simulating the right kinetics at 24 h (black curve). Shaded areas represent ± 1 SD.

The mathematical algorithm to determine liver clearance depends on the time-concentration profiles and, if available, on intracellular or cell-associated compound concentrations. The algorithm minimises a cost function by identifying a clearance value such that model prediction and data observation match. The cost function takes both the PK profile from cell culture media and the cell-associated levels into account, which is not the case for the conventional approach. Further, binding of the compound to plastic/hardware of the chips, to proteins contained in the cell culture media, or any other intracellular lipids can be accounted for to accurately determine liver clearance. DigiLocs, further, also does not depend on a scaling factor, which overcomes the systematic underprediction of conventional approaches (5–10 fold on average across multiple studies).

So far, limited information is available from the literature or in-house measurements on the observed partitioning of compounds into the intracellular or cell-associated milieu of hepatocytes. If that data becomes available, it may be incorporated compared to the adjustments made in the software to match the clinical clearance values. If the predicted and observed Kp_uu values match, the digital twin approach truly improves the prediction. If there is a discrepancy between these values, the fitting process can be re-run including the observed Kp_uu value. This would inform the maximum capacity of the system to metabolise the compound. If this final rate is still lower than the observed clinical clearance value, two options are possible to understand and improve the prediction:

1. Calculate a correction factor, which is compound-specific and chip-specific and not generic like in the conventional approaches.

2. Investigate other model-specific parameters to optimise, e.g., permeability or partitioning.

Although initially developed for hepatic clearance, the mathematical model can be employed for toxicity or efficacy-related questions depending on the context of use. In such a setting, time-concentration profiles will be simulated and linked to other measured biomarkers (e.g., ATP (adenosine triphosphate), TEER (barrier integrity)) to determine IC50 or EC50 values, and parameters to assess toxicity and efficacy, respectively. Likewise, the same integration of complex biological processes, hardware-, and drug-specific information can be used to model other cell and chip types, e.g., a blood-brain-barrier-chip, which is used to determine the permeability of compounds across the barrier.

Eventually, we envision DigiLoCs to support the pharmaceutical decision-making process by reducing animal testing and ultimately streamlining the drug development process.

5 Conclusion

The development of digital twins for organ-on-chips, reported here, incorporating systems of differential equations-based models and leveraging published data, holds great potential to enhance our understanding of drug behaviour and clinical outcomes. The in vitro liver clearance for 32 drugs was predicted using DigiLoCs and a proof-of-concept (translation to human pharmacokinetics) study on propranolol was done. DigiLoCs are envisioned to serve as a decision-support tool for pharmaceutical research, aiding in estimating first-in-human doses, evaluating human pharmacokinetics, and importantly, diminishing reliance on animal experimentation, thereby fostering more efficient, expedited, and sustainable drug development processes. Our approach is generalisable across various physiological contexts and not limited to liver metabolism but may be extended to other organs as well, such as gut metabolism and barrier models such as the brain or placenta.

Supporting information

S1 File. SUPPLEMENTARY INFORMATION FOR DigiLoCS: A leap forward in predictive organ-on-chip simulations.

https://doi.org/10.1371/journal.pone.0314083.s001

References

1. 1. Franzen N, van Harten WH, Retèl VP, Loskill P, van den Eijnden-van Raaij J, IJzerman M. Impact of organ-on-a-chip technology on pharmaceutical R&D costs. Drug Discovery Today. 2019;24(9):1720–1724. pmid:31185290

* View Article

* PubMed/NCBI

* Google Scholar

2. 2. Denayer T, Stöhr T, Roy MV. Animal models in translational medicine: validation and prediction. European Journal of Molecular & Clinical Medicine. 2014;2(1):5.

* View Article

* Google Scholar

3. 3. Thelen K, Coboeken K, Willmann S, Burghaus R, Dressman JB, Lippert J. Evolution of a detailed physiological model to simulate the gastrointestinal transit and absorption process in humans, part 1: oral solutions. Journal of Pharmaceutical Sciences. 2011;100(12):5324–5345. pmid:21993815

* View Article

* PubMed/NCBI

* Google Scholar

4. 4. Jones H, Chen Y, Gibson C, Heimbach T, Parrott N, Peters S, et al. Physiologically based pharmacokinetic modelling in drug discovery and development: a pharmaceutical industry perspective. Clinical Pharmacology & Therapeutics. 2015;97(3):247–262. pmid:25670209

* View Article

* PubMed/NCBI

* Google Scholar

5. 5. Obach RS. Prediction of human clearance of twenty-nine drugs from hepatic microsomal intrinsic clearance data: an examination of in vitro half-life approach and nonspecific binding to microsomes. Drug Metabolism and Disposition. 1999;27(11):1350–1359. pmid:10534321

* View Article

* PubMed/NCBI

* Google Scholar

6. 6. Brown HS, Griffin M, Houston JB. Evaluation of cryopreserved human hepatocytes as an alternative in vitro system to microsomes for the prediction of metabolic clearance. Drug Metabolism and Disposition. 2006;35(2):293–301. pmid:17132764

* View Article

* PubMed/NCBI

* Google Scholar

7. 7. Reddy MB, Mccarley KD, Bunge AL. Physiologically relevant one-compartment pharmacokinetic models for skin. 2. comparison of models when combined with a systemic pharmacokinetic model. Journal of Pharmaceutical Sciences. 1998;87(4):482–490. pmid:9548902

* View Article

* PubMed/NCBI

* Google Scholar

8. 8. Murata Y, Neuhoff S, Rostami-Hodjegan A, Takita H, Al-Majdoub ZM, Ogungbenro K. In vitro to in vivo extrapolation linked to physiologically based pharmacokinetic models for assessing the brain drug disposition. The AAPS Journal. 2022;24(1).

9. 9. Hallifax D, Foster JA, Houston JB. Prediction of human metabolic clearance from in vitro systems: Retrospective analysis and prospective view. Pharmaceutical Research. 2010;27(10):2150–2161. pmid:20661765

* View Article

* PubMed/NCBI

* Google Scholar

10. 10. Shroff T, Aina K, Maass C, Cipriano M, Lambrecht J, Tacke F, et al. Studying metabolism with multi-organ chips: new tools for disease modelling, pharmacokinetics and pharmacodynamics. Open Biology. 2022;12(3). pmid:35232251

* View Article

* PubMed/NCBI

* Google Scholar

11. 11. European Medicines Agency. Guideline on the investigation of drug interactions; 2012. https://www.ema.europa.eu/en/documents/scientific-guideline/guideline-investigation-drug-interactions-revision-1_en.pdf.

12. 12. Maass C, Stokes CL, Griffith LG, Cirit M. Multi-functional scaling methodology for translational pharmacokinetic and pharmacodynamic applications using integrated microphysiological systems (MPS). Integrative Biology. 2017;9(4):290–302. pmid:28267162

* View Article

* PubMed/NCBI

* Google Scholar

13. 13. Docci L, Milani N, Ramp T, Romeo AA, Godoy P, Franyuti DO, et al. Exploration and application of a liver-on-a-chip device in combination with modelling and simulation for quantitative drug metabolism studies. Lab on a Chip. 2022;22(6):1187–1205. pmid:35107462

* View Article

* PubMed/NCBI

* Google Scholar

14. 14. Fowler S, Chen WLK, Duignan DB, Gupta A, Hariparsad N, Kenny JR, et al. Microphysiological systems for ADME-related applications: current status and recommendations for system development and characterization. Lab on a Chip. 2020;20(3):446–467. pmid:31932816

* View Article

* PubMed/NCBI

* Google Scholar

15. 15. Maass C, Sorensen NB, Himmelfarb J, Kelly EJ, Stokes CL, Cirit M. Translational assessment of drug‐induced proximal tubule injury using a kidney microphysiological system. CPT: Pharmacometrics & Systems Pharmacology. 2019;8(5):316–325. pmid:30869201

* View Article

* PubMed/NCBI

* Google Scholar

16. 16. Herland A, Maoz BM, Das D, Somayaji MR, Prantil-Baun R, Novak R, et al. Quantitative prediction of human pharmacokinetic responses to drugs via fluidically coupled vascularized organ chips. Nature Biomedical Engineering. 2020;4(4):421–436. pmid:31988459

* View Article

* PubMed/NCBI

* Google Scholar

17. 17. Tsamandouras N, Kostrzewski T, Stokes CL, Griffith LG, Hughes DJ, Cirit M. Quantitative assessment of population variability in hepatic drug metabolism using a perfused three-dimensional human liver microphysiological system. Journal of Pharmacology and Experimental Therapeutics. 2016;360(1):95–105. pmid:27760784

* View Article

* PubMed/NCBI

* Google Scholar

18. 18. Kanebratt KP, Janefeldt A, Vilén L, Vildhede A, Samuelsson K, Milton L, et al. Primary human hepatocyte spheroid model as a 3D in vitro platform for metabolism studies. Journal of Pharmaceutical Sciences. 2021;110(1):422–431. pmid:33122050

* View Article

* PubMed/NCBI

* Google Scholar

19. 19. Bonn B, Svanberg P, Janefeldt A, Hultman I, Grime K. Determination of human hepatocyte intrinsic clearance for slowly metabolized compounds: comparison of a primary hepatocyte/stromal cell co-culture with plated primary hepatocytes and HepaRG. Drug Metabolism and Disposition. 2016;44(4):527–533. pmid:26851239

* View Article

* PubMed/NCBI

* Google Scholar

20. 20. Rajan SAP, Sherfey J, Ohri S, Nichols L, Smith JT, Parekh P, et al. A novel milli-fluidic liver tissue chip with continuous recirculation for predictive pharmacokinetics applications. The AAPS Journal. 2023;25(6). pmid:37891356

* View Article

* PubMed/NCBI

* Google Scholar

21. 21. R Core Team. R: A language and environment for statistical computing; 2021. Available from: https://www.R-project.org/.

22. 22. Poulin P, Haddad S. Hepatocyte composition-based model as a mechanistic tool for predicting the cell suspension: aqueous phase partition coefficient of drugs in in vitro metabolic studies. Journal of Pharmaceutical Sciences. 2013;102(8):2806–2818. pmid:23670739

* View Article

* PubMed/NCBI

* Google Scholar

23. 23. Austin RP, Barton P, Mohmed S, Riley RJ. The binding of drugs to hepatocytes and its relationship to physiochemical properties. Drug Metabolism and Disposition. 2004;33(3):419–425. pmid:15616151

* View Article

* PubMed/NCBI

* Google Scholar

24. 24. Mateus A, Matsson P, Artursson P. Rapid measurement of intracellular unbound drug concentrations. Molecular Pharmaceutics. 2013;10(6):2467–2478. pmid:23631740

* View Article

* PubMed/NCBI

* Google Scholar

25. 25. Bhatia SN, Ingber DE. Microfluidic organs-on-chips. Nat Biotechnol. 2014;32(8):760–772. pmid:25093883

* View Article

* PubMed/NCBI

* Google Scholar

26. 26. Wickham H. Ggplot2. 2nd ed. Use R!. Cham, Switzerland: Springer International Publishing; 2016.

27. 27. Kuepfer L, Niederalt C, Wendl T, Schlender J, Willmann S, Lippert J, et al. Applied concepts in PBPK modelling: how to build a PBPK/PD model. CPT: Pharmacometrics & Systems Pharmacology. 2016;5(10):516–531. pmid:27653238

* View Article

* PubMed/NCBI

* Google Scholar

28. 28. Rodrigues AV, Alexandre-Pires G, Valério-Bolas A, Santos-Mateus D, Rafael-Fernandes M, Pereira MA, et al. 3D-hepatocyte culture applied to parasitology: immune activation of canine hepatic spheroids exposed to leishmania infantum. Biomedicines. 2020;8(12):628. pmid:33352885

* View Article

* PubMed/NCBI

* Google Scholar

29. 29. Borgström L, Johansson CG, Larsson H, Lenander R. Pharmacokinetics of propranolol. Journal of Pharmacokinetics and Biopharmaceutics. 1981;9(4):419–429. pmid:7310641

* View Article

* PubMed/NCBI

* Google Scholar

Citation: Aravindakshan MR, Mandal C, Pothen A, Schaller S, Maass C (2025) DigiLoCS: A leap forward in predictive organ-on-chip simulations. PLoS ONE 20(1): e0314083. https://doi.org/10.1371/journal.pone.0314083

About the Authors:

Manoja Rajalakshmi Aravindakshan

Roles: Data curation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

Affiliation: Department of Computer Science and Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal, India

ORICD: https://orcid.org/0000-0003-2861-4453

Chittaranjan Mandal

Roles: Methodology, Supervision, Writing – review & editing

Affiliation: Department of Computer Science and Engineering, Indian Institute of Technology Kharagpur, Kharagpur, West Bengal, India

Alex Pothen

Roles: Supervision, Writing – review & editing

Affiliation: Department of Computer Science, Purdue University, West Lafayette, Indiana, United States

Stephan Schaller

Roles: Writing – review & editing

Affiliation: ESQlabs Gmbh, Saterland, Germany

Christian Maass

Roles: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Supervision, Visualization, Writing – original draft, Writing – review & editing

E-mail: [email protected]

Affiliations: ESQlabs Gmbh, Saterland, Germany, MPSlabs, ESQlabs Gmbh, Saterland, Germany

ORICD: https://orcid.org/0000-0001-7618-0421

[/RAW_REF_TEXT]

[/RAW_REF_TEXT]

[/RAW_REF_TEXT]

[/RAW_REF_TEXT]

References

1. Franzen N, van Harten WH, Retèl VP, Loskill P, van den Eijnden-van Raaij J, IJzerman M. Impact of organ-on-a-chip technology on pharmaceutical R&D costs. Drug Discovery Today. 2019;24(9):1720–1724. pmid:31185290

2. Denayer T, Stöhr T, Roy MV. Animal models in translational medicine: validation and prediction. European Journal of Molecular & Clinical Medicine. 2014;2(1):5.

3. Thelen K, Coboeken K, Willmann S, Burghaus R, Dressman JB, Lippert J. Evolution of a detailed physiological model to simulate the gastrointestinal transit and absorption process in humans, part 1: oral solutions. Journal of Pharmaceutical Sciences. 2011;100(12):5324–5345. pmid:21993815

4. Jones H, Chen Y, Gibson C, Heimbach T, Parrott N, Peters S, et al. Physiologically based pharmacokinetic modelling in drug discovery and development: a pharmaceutical industry perspective. Clinical Pharmacology & Therapeutics. 2015;97(3):247–262. pmid:25670209

5. Obach RS. Prediction of human clearance of twenty-nine drugs from hepatic microsomal intrinsic clearance data: an examination of in vitro half-life approach and nonspecific binding to microsomes. Drug Metabolism and Disposition. 1999;27(11):1350–1359. pmid:10534321

6. Brown HS, Griffin M, Houston JB. Evaluation of cryopreserved human hepatocytes as an alternative in vitro system to microsomes for the prediction of metabolic clearance. Drug Metabolism and Disposition. 2006;35(2):293–301. pmid:17132764

7. Reddy MB, Mccarley KD, Bunge AL. Physiologically relevant one-compartment pharmacokinetic models for skin. 2. comparison of models when combined with a systemic pharmacokinetic model. Journal of Pharmaceutical Sciences. 1998;87(4):482–490. pmid:9548902

8. Murata Y, Neuhoff S, Rostami-Hodjegan A, Takita H, Al-Majdoub ZM, Ogungbenro K. In vitro to in vivo extrapolation linked to physiologically based pharmacokinetic models for assessing the brain drug disposition. The AAPS Journal. 2022;24(1).

9. Hallifax D, Foster JA, Houston JB. Prediction of human metabolic clearance from in vitro systems: Retrospective analysis and prospective view. Pharmaceutical Research. 2010;27(10):2150–2161. pmid:20661765

10. Shroff T, Aina K, Maass C, Cipriano M, Lambrecht J, Tacke F, et al. Studying metabolism with multi-organ chips: new tools for disease modelling, pharmacokinetics and pharmacodynamics. Open Biology. 2022;12(3). pmid:35232251

11. European Medicines Agency. Guideline on the investigation of drug interactions; 2012. https://www.ema.europa.eu/en/documents/scientific-guideline/guideline-investigation-drug-interactions-revision-1_en.pdf.

12. Maass C, Stokes CL, Griffith LG, Cirit M. Multi-functional scaling methodology for translational pharmacokinetic and pharmacodynamic applications using integrated microphysiological systems (MPS). Integrative Biology. 2017;9(4):290–302. pmid:28267162

13. Docci L, Milani N, Ramp T, Romeo AA, Godoy P, Franyuti DO, et al. Exploration and application of a liver-on-a-chip device in combination with modelling and simulation for quantitative drug metabolism studies. Lab on a Chip. 2022;22(6):1187–1205. pmid:35107462

14. Fowler S, Chen WLK, Duignan DB, Gupta A, Hariparsad N, Kenny JR, et al. Microphysiological systems for ADME-related applications: current status and recommendations for system development and characterization. Lab on a Chip. 2020;20(3):446–467. pmid:31932816

15. Maass C, Sorensen NB, Himmelfarb J, Kelly EJ, Stokes CL, Cirit M. Translational assessment of drug‐induced proximal tubule injury using a kidney microphysiological system. CPT: Pharmacometrics & Systems Pharmacology. 2019;8(5):316–325. pmid:30869201

16. Herland A, Maoz BM, Das D, Somayaji MR, Prantil-Baun R, Novak R, et al. Quantitative prediction of human pharmacokinetic responses to drugs via fluidically coupled vascularized organ chips. Nature Biomedical Engineering. 2020;4(4):421–436. pmid:31988459

17. Tsamandouras N, Kostrzewski T, Stokes CL, Griffith LG, Hughes DJ, Cirit M. Quantitative assessment of population variability in hepatic drug metabolism using a perfused three-dimensional human liver microphysiological system. Journal of Pharmacology and Experimental Therapeutics. 2016;360(1):95–105. pmid:27760784

18. Kanebratt KP, Janefeldt A, Vilén L, Vildhede A, Samuelsson K, Milton L, et al. Primary human hepatocyte spheroid model as a 3D in vitro platform for metabolism studies. Journal of Pharmaceutical Sciences. 2021;110(1):422–431. pmid:33122050

19. Bonn B, Svanberg P, Janefeldt A, Hultman I, Grime K. Determination of human hepatocyte intrinsic clearance for slowly metabolized compounds: comparison of a primary hepatocyte/stromal cell co-culture with plated primary hepatocytes and HepaRG. Drug Metabolism and Disposition. 2016;44(4):527–533. pmid:26851239

20. Rajan SAP, Sherfey J, Ohri S, Nichols L, Smith JT, Parekh P, et al. A novel milli-fluidic liver tissue chip with continuous recirculation for predictive pharmacokinetics applications. The AAPS Journal. 2023;25(6). pmid:37891356

21. R Core Team. R: A language and environment for statistical computing; 2021. Available from: https://www.R-project.org/.

22. Poulin P, Haddad S. Hepatocyte composition-based model as a mechanistic tool for predicting the cell suspension: aqueous phase partition coefficient of drugs in in vitro metabolic studies. Journal of Pharmaceutical Sciences. 2013;102(8):2806–2818. pmid:23670739

23. Austin RP, Barton P, Mohmed S, Riley RJ. The binding of drugs to hepatocytes and its relationship to physiochemical properties. Drug Metabolism and Disposition. 2004;33(3):419–425. pmid:15616151

24. Mateus A, Matsson P, Artursson P. Rapid measurement of intracellular unbound drug concentrations. Molecular Pharmaceutics. 2013;10(6):2467–2478. pmid:23631740

25. Bhatia SN, Ingber DE. Microfluidic organs-on-chips. Nat Biotechnol. 2014;32(8):760–772. pmid:25093883

26. Wickham H. Ggplot2. 2nd ed. Use R!. Cham, Switzerland: Springer International Publishing; 2016.

27. Kuepfer L, Niederalt C, Wendl T, Schlender J, Willmann S, Lippert J, et al. Applied concepts in PBPK modelling: how to build a PBPK/PD model. CPT: Pharmacometrics & Systems Pharmacology. 2016;5(10):516–531. pmid:27653238

28. Rodrigues AV, Alexandre-Pires G, Valério-Bolas A, Santos-Mateus D, Rafael-Fernandes M, Pereira MA, et al. 3D-hepatocyte culture applied to parasitology: immune activation of canine hepatic spheroids exposed to leishmania infantum. Biomedicines. 2020;8(12):628. pmid:33352885

29. Borgström L, Johansson CG, Larsson H, Lenander R. Pharmacokinetics of propranolol. Journal of Pharmacokinetics and Biopharmaceutics. 1981;9(4):419–429. pmid:7310641