A wide range of heterogeneous data was curated for model building (parameterization) and subsequent model validation (comparison of model predictions to clinical data). An overview of the 29 studies with their respective ICG dosing protocols is provided in Table 1. All data is available via the pharmacokinetics database PK-DB (https://pk-db.com) (Grzegorzewski et al., 2021). PK-DB was used to encode the information on (i) patient characteristics (e.g., age, disease, medication), (ii) applied interventions (e.g., ICG dosing, route of application); (iii) measured ICG time courses in plasma and (iv) ICG pharmacokinetic parameters (e.g., ICG-PDR, ICG-R15, ICG-clearance). The established clinical data base was essential for model building (see section 3.1.3) and model validation (see sections 3.2, 3.3, and 3.4).

The computational model is an ordinary differential equation (ODE) model without delays encoded in the Systems Biology Markup Language (SBML) (Hucka et al., 2019; Keating et al., 2020). It is defined as a set of species (metabolites), compartments (organs and blood compartments), and reactions (processes such as metabolic reactions and blood transport). The model was developed using sbmlutils (König, 2021b), and cy3sbml (König and Rodriguez, 2019) and simulated using sbmlsim (König, 2021a) based on the high-performance SBML simulator libroadrunner (Somogyi et al., 2015).

Values for organ volumes and tissue blood flows were taken from literature (ICRP, 2002; Jones and Rowland-Yeo, 2013). A total of five model parameters were determined by minimizing the residuals between model predictions and clinical data. This optimization-problem was solved using SciPy's least_squares method and differential evolution algorithm (Virtanen et al., 2020). For the objective cost function F depending on the parameters p→ a simple L2-Norm was used consisting of the sum of weighted residuals

where ri,k=(yi,k−mi,k(p→)) is the residual of time point i in time course k for model prediction mi,k(p→) and the corresponding data point y_{i, k}; w_{i, k} is the weighting of the respective data point i in time course k based on the error of the data point and w_k = the weighting factor of time course k. Weighting of time courses was based on the number of subjects per study. The data used for the parameter fit is listed in Table 1. The final parameter set given in Table 2 was determined using 250 runs of the local least squares optimization.

To evaluate the uncertainty of model predictions uncertainty analysis was performed for a subset of simulations. Each model parameter was changed individually by B±25%. From the set of resulting time courses the mean, standard deviation (SD) and minimum and maximum values at each time point were calculated. These uncertainty areas were displayed as shaded areas. Parameters corresponding to physical constants (such as molecular weights) and dosing were not varied in the uncertainty analysis, as well as parameters for conservation conditions such as the fractional blood flow through the lung (must be 1).

Pharmacokinetic parameters of ICG were calculated from the plasma-concentration time courses using non-compartmental methods (Urso et al., 2002). The elimination rate constant (k_el) was calculated by fitting the concentration-decay-curve to an exponential function: c(t)=c(0)·e-kel·t. ICG-PDR is k_el reported in [%/min]. Half-life (t_1/2) was calculated as log(2)/k_el, ICG-clearance as CL = V_d·k_el, with the apparent volume of distribution V_d = D/(AUC_∞·kel). D is the applied dose of ICG and AUC_∞ is the area under the plasma-concentration time curve AUC calculated via the trapezoidal rule, extrapolated until infinity. ICG-R15 = c(15)/c_max was calculated as the ratio between the plasma-concentration after 15 min and the maximum concentration c_max.

Classification models were developed which predict survival after hepatectomy (binary classification: Survivors/Non-Survivors) based on different sets of features (see Table 3). Classification was performed using scikit-learn (Pedregosa et al., 2011) with a C-support vector classifier using a polynomial kernel. For model training and evaluation a dataset of 141 patients with information on survival status, resection rate and preoperative ICG-R15 was used (Wakabayashi et al., 2004; Seyama and Kokudo, 2009). Cross-validation was performed using the ShuffleSplit method with 200 iterations and a train-test-ratio of 75/25%. Based on the confusion matrix the following evaluation metrics were calculated: precision/positive predictive value (PPV), recall, balanced accuracy, F1 score, negative predictive value (NPV), and receiver operator curves (ROC).

For the classification models PBPK1 and PBPK2 the model parameter f_cirrhosis was determined from preoperative ICG-R15 values using

with R15_preop being the clinically measured preoperative ICG-R15 [dimensionless]. The coefficients a = 0.312, b = 1.693, and c = 0.861 were determined by fitting the curve to the predicted dependency of ICG-R15 on f_cirrhosis (see Figure 4C). This provided an estimate of individual liver disease (cirrhosis degree). This estimated parameter allowed in combination with the resection rate to predict individual postoperative ICG-R15 values.

A PBPK model for the prediction of ICG pharmacokinetics was developed which is available in SBML (Hucka et al., 2019; Keating et al., 2020) under CC-BY 4.0 license from https://github.com/matthiaskoenig/icg-model. Within this work version 1.0 of the model was used (König and Köller, 2021). To simulate the whole-body distribution and hepatic elimination of ICG two models were coupled: (i) A whole-body model (Figure 1A) describing the distribution of ICG in the body and to the organs via blood flow. (ii) A liver model (Figure 1B) which describes hepatic uptake of ICG, biliary excretion of ICG and transport of ICG into the feces. ICG pharmacokinetic parameters (i.e., ICG-PDR, R15, clearance, half-life) were calculated from the resulting time course predictions of ICG in the venous plasma.

In a first step the fitted model was evaluated with the data used for model calibration consisting of ICG time courses in healthy subjects (Figure 2). For the simulations infusion protocols and body weights were adjusted as reported in the respective studies (see Table 1 for details). If no body weight was reported 75 kg were assumed. The same unique set of fitted parameters was used for all simulations (see section 3.1.3).

The model predictions for ICG plasma disappearance curves after an ICG bolus are in good agreement with the clinical data (Meijer et al., 1988; Klockowski et al., 1990; Grundmann et al., 1992; Andersen et al., 1999; Kamimori et al., 2000; Niemann et al., 2000) (Figures 2A–E). In addition, more complex infusion protocols as reported in Soons et al. (1991) can also be described (Figure 2F), infusion protocol of three different infusion rates 2.0 → 0.5 → 1.0 mg/min, each for 40 min). Due to the high extraction-ratio of ICG by the liver, the plasma concentration reaches steady state quickly after each change in the infusion rate. Next, simulations of the biliary excretion rate of ICG after bolus administrations of 0.5, 1.0, and 2.0 mg/kg ICG were performed and the results were compared to clinical data (Meijer et al., 1988; Chijiiwa et al., 2000) (Figures 2G,H).

Finally, simulations of constant infusions were performed and compared to reported arterial and hepatic vein time courses of ICG (Leevy et al., 1962) and ICG extraction ratios (Leevy et al., 1962; Grainger et al., 1983) (Figures 2I–L).

Overall, the model shows the ability to accurately predict ICG time courses for venous and arterial plasma concentrations, for hepatic vein concentrations, the biliary excretion rate and extraction ratios when compared to clinical data. Especially plasma time courses of ICG after ICG bolus and ICG infusion are very well-predicted by the model, even for varying administration protocols (dosing and infusion rates).

Discrepancies between model predictions and data can be observed in case of the ICG extraction ratios, where data points are outside of the uncertainty bounds of the model simulation in Figures 2J–L, and the biliary excretion in Figure 2H at 0.5 mg/kg (blue) and 1.0 mg/kg (green), whereas biliary excretion at 2.0 mg/kg (red) is in good agreement with the model.

In a next step, a systematic analysis of the dose dependency of ICG pharmacokinetic parameters was performed (Figure 3). A dose-dependency of the ICG parameters can only be observed if the ICG dose exceeds 100 mg (much higher then the typically applied doses of 20–35 mg), resulting in a reduction in ICG-clearance and ICG-PDR as well as an increase of ICG-R15 and ICG-t_1/2 (Figures 3A–D). The model predictions could be validated with clinical data (Martin et al., 1975, 1976; Meijer et al., 1988) (Figures 3E–G).

To simulate changes of ICG pharmacokinetics in cirrhosis, hepatic tissue loss and shunts were included in the model as described above. The same unique set of fitted parameters used for the validation in healthy subjects was used here. First, a systematic analysis of the effect of intrahepatic shunts (f_shunts), functional tissue loss (f_{tissue_loss}) and the combination of both (f_cirrhosis) on ICG pharmacokinetic parameters was performed (Figures 4A–D). All three parameters were varied from 0 (no effect, healthy control) to 0.9 (severe effect). ICG-clearance and ICG-PDR decrease with increasing f_cirrhosis whereas ICG-R15 and ICG-t_1/2 increase. The loss of a fraction of functional liver tissue appears to have a smaller effect on ICG pharmacokinetic parameters than shunting of an equal fraction of blood past the liver. When f_shunts and f_{tissue_loss} are combined to f_cirrhosis their effect on ICG pharmacokinetic parameters is additive. For ICG-clearance and ICG-PDR the effects of both parameters combine to an almost linear dependency on f_cirrhosis (Figures 4A,B).

By varying the f_cirrhosis parameter from 0 to 0.9 different degrees of cirrhosis were simulated and the nonlinear relation between ICG-R20 and ICG-kel as well as ICG-R20 and ICG-t_1/2 could be predicted (Figures 4E,F). As seen in the systematic analysis (Figures 4A–D) ICG-t_1/2 and ICG-R20 increase with cirrhosis whereas ICG-kel decreases. The correlation between the ICG pharmacokinetic parameters is predicted accurately by the model when compared to a clinical dataset that lacks information about the severity of liver cirrhosis of its patients (Cherrick et al., 1960; Caesar et al., 1961). Next the ICG-PDR in cirrhotic patients and control subjects after different doses of ICG (0.5 and 5.0 mg/kg ICG) was compared to the model predictions. The clinical data shows higher ICG-PDR values after an ICG dose of 0.5 mg/kg than after an ICG dose of 5.0 mg/kg (Leevy et al., 1967) (Figure 4G). ICG-clearances after a bolus administration and during a constant infusion show good positive correlation in cirrhotic patients (Burns et al., 1991). This correlation is predicted accurately by the model (Figure 4H).

Having evaluated and validated the effect of f_cirrhosis on the model prediction of ICG parameters, we were interested how the model f_cirrhosis parameter compares to the in vivo estimation of cirrhosis degree via the CTP-score (Figure 5). As described in the introduction, the CTP-score is a semi-quantitative scoring system that describes the severity of liver cirrhosis. An important step to apply the developed PBPK model in a clinical setting, is the ability to adjust the model individually to the respective status of liver disease in a patient. Therefore, the relationship between the f_cirrhosis parameter and the CTP-Score was evaluated using multiple datasets in which ICG pharmacokinetic parameters were reported in patient subgroups of different CTP-Scores (Figg et al., 1995; Møller et al., 1998, 2019; Herold et al., 2001).

The clinical results of the ICG pharmacokinetic parameters in different CTP-classes were mapped onto their respective systematic scan (Figures 5A–D). The resulting f_cirrhosis values were then compared between the patient groups. Additional individual data is shown (Figg et al., 1995) for validation.

The resulting mapping between f_cirrhosis and the CTP-classes shows a good positive correlation (Figure 5E). The f_cirrhosis values for the controls groups are close to 0, increasing with the CTP-class. The relation appears nonlinear, as f_cirrhosis shows little difference between CTP-class B and C. The mappings of the CTP-class to f_cirrhosis for the different ICG parameters each give very similar results. From the mapping of all four pharmacokinetic parameters a mean value of f_cirrhosis was calculated for each CTP-class (Control: f_cirrhosis = 0.0; Mild cirrhosis: 0.41; Moderate cirrhosis: 0.70; Severe cirrhosis: 0.82). The resulting values were used in all simulations of control, mild, moderate and severe cirrhosis, as well as in the above described dose dependency analysis (Figures 3A–D).

Only a single study reported the numerical CTP-score of the patient groups in combination with ICG-clearance (Figg et al., 1995). All other studies instead used the CTP-classes (A, B, C) (Møller et al., 1998, 2019; Herold et al., 2001). With a dataset of individually reported CTP-scores in combination with ICG pharmacokinetic parameters of cirrhotic patients, it would be possible to calculate the relationship of the CTP-score on the f_cirrhosis parameter more accurately. Such an improved mapping would allow to adjust the model via the f_cirrhosis parameter individually based on the respective severity of liver disease/cirrhosis of the patient reported as CTP-score.

After establishing the CTP mapping the model was further validated via several comparisons with clinical data of ICG time courses in cirrhotic and control subjects (Figure 6).

Assuming moderate cirrhosis (f_cirrhosis = 0.7), the model prediction of an ICG time course in a cirrhotic patient agrees well with the clinical data (Burns et al., 1991) (Figures 6A,B). The main alteration compared to the healthy control is the slower disappearance rate resulting in higher ICG plasma concentrations. The same effect is observed in steady state via a constant ICG infusion (Figure 6C). Using the f_cirrhosis values from the CTP mapping above, the steady state concentrations are predicted in agreement with the clinical data (Caesar et al., 1961). Figures 6D,E shows the relation between the hepatic venous and arterial ICG concentrations and the extraction ratio in a cirrhotic subject. Here, f_cirrhosis was set to 0.54 which allowed to predict arterial and hepatic vein concentration as well as ICG extraction ratio. Finally, in Figures 6F–H the ICG extraction ratio predicted for controls and three different cirrhosis degrees was compared to clinical data (Caesar et al., 1961; Leevy et al., 1962; Gadano et al., 1997). The extraction ratio in cirrhotic subjects is reduced compared to healthy controls, as predicted by the model.

After validating the model predictions of ICG pharmacokinetics in liver cirrhosis, the model was applied to liver surgery. To analyze the effect of partial hepatectomy on ICG elimination the change in ICG pharmacokinetic parameters as a function of the resection rate was simulated (Figures 7A–D). The scan was performed for healthy controls as well as three different degrees of cirrhosis.

ICG-clearance and ICG-PDR are highest in the preoperative liver (resection rate = 0) and decrease with increasing resection rate whereas ICG-t_1/2 and ICG-R15 are lowest in the healthy liver and increase with increasing resection rate. The effect of varying the degree of cirrhosis is in accordance with the results shown in Figures 4A–D. Importantly, increasing resection rate and increasing degree of cirrhosis affect ICG pharmacokinetic parameters in the same manner. The dependencies of ICG-clearance, ICG-PDR, ICG-t_1/2 and ICG-R15 on the resection rate are fairly linear up to 50–60% resection, and become much more non-linear for higher resection rates.

For model validation the predictions were compared to clinical data of subjects undergoing hepatectomy. For these simulations the resection rate was varied from 0 to 0.9. First, the relative change of ICG-PDR after partial hepatectomy as a function of the resection rate was simulated (Figure 7E). The model predicts a nonlinear dependency of change in ICG-PDR on the remnant liver volume independent of the degree of cirrhosis. This prediction is in good agreement with the clinical data (Stockmann et al., 2009; Thomas et al., 2015). Furthermore, the correlation between measured postoperative ICG-kel and estimated remnant ICG-kel (ICG-kel · fractional liver remnant) was simulated under various degrees of cirrhosis (Figure 7F). A good correlation can be observed. The model predictions were compared to three different data sets (Okochi et al., 2002; Ohwada et al., 2006; Sunagawa et al., 2021) and are in good agreement with them. In addition, all data sets are in good agreement with each other. The simulated correlation line is independent of the cirrhosis degree, but with increasing cirrhosis ICG-kel decreases. A large variability can be observed in the experimental data, but as our simulations indicate is most likely not due to the underlying liver disease (cirrhosis).

Thomas et al. (2015) found significant correlation between post-hepatectomy ICG-PDR and intraoperative ICG-PDR measured under trial clamping of those parts of the liver that were to be removed. This was simulated by changing hepatic blood flow and liver volume in separate simulations but in the same intervals. This was performed for a healthy liver as well as three different degrees of cirrhosis (Figure 7G). The predictions agree well with the clinical data and show that reducing hepatic blood flow (clamping of liver volumes which will be resected) has a very similar effect on ICG elimination as actually removing the respective liver volume via hepatectomy.

Finally, the correlation between preoperative and postoperative ICG-PDR for different resection rates and cirrhosis degrees was simulated and compared to clinical data (Figure 7H). ICG-PDR is reduced in cirrhosis preoperatively as well as postoperatively. The model prediction agrees with the clinical data (Thomas et al., 2015).

Overall the predictions of liver resection in severely cirrhotic liver is not in good agreement with the clinical data. This reflects the fact that no resections are performed in severely cirrhotic liver due to high risk of postoperative complications. As a consequence, most of the liver resection are performed in mild to moderate cirrhosis. The model allows to perform these risky hepatectomies in silico and predict there effect.

In summary, the model allows to systematically predict the changes of ICG pharmacokinetic parameters in HPB surgery under various degrees of liver disease (cirrhosis).

An interesting application of the presented PBPK model is the prediction of postoperative outcome for patients undergoing hepatectomy. Preoperative ICG-R15 and the planned resection rate are key parameters included in the decision process whether a patient is eligible to receive liver resection surgery.

As shown above, the presented PBPK model accurately predicts ICG-R15 in liver cirrhosis as well as the changes in ICG-R15 following hepatectomy. As such, we were interested how a classification model based on the PBPK model prediction of postoperative ICG-R15 compares to classification approaches only using clinical data (preoperative ICG-R15, resection rate and calculated postoperative ICG-R15).

Five alternative classification models were developed to predict survival after partial hepatectomy using a dataset of 141 patients (109 survivors and 32 non-survivors) (Wakabayashi et al., 2004; Seyama and Kokudo, 2009). Three of the classification models were solely based on clinical data: (i) Data1A–using preoperative ICG-R15, (ii) Data1B–using the calculated postoperative ICG-R15 by multiplying the future liver remnant (1-resection rate) with the preoperative ICG-R15; and (iii) Data2–using both the resection rate and the preoperative ICG-R15. In addition two classification models were developed using PBPK model predictions as input: (iv) PBPK1–using the predicted postoperative ICG-R15. Hereby, the model parameter f_cirrhosis was estimated for every subject based on preoperative ICG-R15, and postoperative ICG-R15 was predicted using f_cirrhosis and the corresponding resection rate; and (v) PBPK2– using the resection rate and the estimated f_cirrhosis model parameter. An overview of the classification results of these five models is provided in Table 3, Figure 8, and Supplementary Figure 1.

Both PBPK-based classifiers (PBPK1, PBPK2) as well as the Data2 classifier outperform the data-based classifiers using a single feature (Data1A, Data1B) in predicting survival after partial hepatectomy.

When comparing the classification models using a single feature (Data1A, Data1B, PBPK1) the physiological-based predicted postoperative ICG-R15 (PBPK1) clearly outperforms the preoperative (Data1A) as well as calculated postoperative ICG-R15 (Data1B).

Figure 8A shows the postoperative ICG-R15 in survivors and non-survivors predicted by the model as well as corresponding preoperative ICG-R15 (inset).

The PBPK model predicted postoperative ICG-R15 is able to distinguish better between survivors and non-survivors than the preoperative ICG-R15 as can be seen by the clearer separation of the histograms (Figure 8A) and the respective ROC curves (Figure 8C). Both preoperative ICG-R15 as well as calculated postoperative ICG-R15 are not very useful for the prediction of survival after partial hepatectomy, whereas predicted postoperative ICG-R15 using PBPK1 is a very good measure to predict survival after partial hepatectomy.

To determine possible cutoffs for predicted postoperative ICG-R15 based on the PBPK1 classifier the dependency of evaluation metrics on the cutoff was analyzed (Figure 8B). Balanced accuracy as well as f1-score have a maximum at around 35%. The negative predictive value (NPV) as well as recall is 1 up to a predicted postoperative ICG-R15 < 20%.

Figure 8D depicts how the predicted postoperative ICG-R15 depends on the resection rate and f_cirrhosis. The data confirms that a cutoff value slightly below 40% would correctly predict most of the non-survivors with a stricter cutoff of 20% avoiding any death after partial hepatectomy. Similar analysis of the data-based single feature classification models failed to find a significant optimum of evaluation metrics for either preoperative ICG-R15 or calculated postoperative ICG-R15 (see Supplementary Figure 1).

The two-feature classification models (PBPK2, Data2) show good performance in the survival prediction comparable to PBPK1 (Figures 8E,F). Whereas, the one-dimensional PBPK1 classifier provides a simple interpretation and cutoff value, the two dimensional classifiers are more difficult to interpret.

In summary, we developed a single-feature classification model based on a physiological-based model of ICG elimination (PBPK1) which allows to predict post-hepatectomy survival solely based on preoperative ICG-R15 input. Importantly, this computational model-based approach clearly outperforms data-based approaches such as preoperative ICG-R15 and calculated postoperative ICG-R15.

In summary, a PBPK model for ICG based liver function evaluation was developed, validated, and applied to the prediction of postoperative outcome after liver surgery, i.e., survival after partial hepatectomy. The model takes into account physiological factors such as the degree of cirrhosis and the estimated functional liver remnant, which allowed an accurate prediction of postoperative liver function in agreement with clinical data. As such, the model has proven its potential of becoming a valuable clinical tool for the planning of hepato-pancreatico-biliary surgery.

The physiologically-based modeling approach allowed us to predict ICG pharmacokinetics data from 29 studies (see Table 1) using only a small set of parameters and processes. The model accurately predicts changes in ICG pharmacokinetic parameters in a wide range of conditions including varying degrees of cirrhosis. Additionally, in silico hepatectomies with underlying cirrhosis are in good agreement with clinical data. As an important note, all clinical data besides the time courses in healthy subjects used for model calibration was used for model validation.

Important model assumptions were (i) the pooling of all tissues not involved in the metabolization of ICG into the rest compartment, (ii) the use of ordinary differential equations without delays, and (iii) the focus on the first phase of ICG elimination relevant for the calculation of ICG parameters. By pooling all other tissues in a single compartment, tissues with different perfusion rates were combined. Due to transit times in the vessel trees of each organ a delay between the flow supplying the organ and draining from the organ occurs. A possible solution to address this issue is to use delayed differential equations which explicitly set these delays in the model. The model did not include any delays, but instead the cardiac output was set at the lower end of the physiological range to correctly describe the observed hepatic extraction ratios. Using an ODE model without delays allows for much faster simulation times as well as exchangeable and reusability of the model in other contexts which do not support delayed differential equation models. The model focuses on the accurate description of the first 15–30 min of ICG elimination relevant for the calculation of the ICG parameters, which is well-described by a monophasic disappearance pattern. In reality, a biphasic disappearance pattern of ICG with a second, slower phase of plasma elimination can be observed (Stoeckel et al., 1980; Meijer et al., 1988; Burns et al., 1991). This phase is not relevant for the calculation of pharmacokinetic parameters and the presented model focuses on the first phase of ICG elimination.

In addition, some discrepancies between model predictions and clinical data existed for the biliary excretion and the ICG extraction ratio. The focus of the model was on accurate prediction of plasma ICG time courses and ICG parameters. The use of a simple biliary excretion model was not able to capture all the observed dynamics.

Importantly, all model simulations were performed with an identical model parameter set not accounting for inter-individual or inter-study differences or biases.

An important result of this work is that the presented model is able to predict clinical data ranging from ICG timecourses after various dosing regimes to changes in ICG parameters after partial hepatectomy. Although the model predictions appear to be in very good agreement with the data, this agreement was not analyzed systematically but was limited to a graphical interpretation. A more thorough statistical analysis would be required in the future including parameter identifiability and uncertainty analysis.

One main outcome of this study is a single-feature classification model based on a physiological-based model of ICG elimination (PBPK1) which allows to predict post-hepatectomy survival solely based on preoperative ICG-R15 and resection rate. Further, it has become apparent that preoperative measurements of ICG pharmacokinetics alone are not sufficient for predicting postoperative outcome.

The developed classification models demonstrated the potential of using PBPK predicted postoperative ICG-R15 values in the clinical decision process. Whereas, the PBPK1 classifier provides a simple cutoff based on the individual model prediction (Figures 8A,B), PBPK2 provides the dependency of predicted postoperative ICG-R15 on resection rate and cirrhosis degree (Figure 8D), both key factors for survival after liver resection.

Comparing different approaches of predicting postoperative outcome after partial hepatectomy showed the importance of taking resection rate into account. The data-based classifier combining resection rate with preoperative ICG (Data2) allowed to achieve comparable classification results to the PBPK-based classifiers. In contrast, the classification models Data1A and Data1B failed to achieve good prediction results.

This work is the first published PBPK model of ICG. The use of model-predicted ICG pharmacokinetic parameters to predict postoperative survival after partial hepatectomy using a classification model is a novel approach. Two-compartment models have been used among other things to further the understanding of ICG kinetics (Stoeckel et al., 1980) or for the development of an at the time new method of determining the ICG-clearance and hepatic extraction ratio (Grainger et al., 1983). The use of a PBPK model facilitated studying more complex clinical questions such as the effect of partial hepatectomy on ICG elimination in patients with underlying liver cirrhosis of different severity. Furthermore, for the potential clinical application of the presented model, individualizing it to a patient's physiology via the model's physiological parameters would be essential, which is a strong argument for the PBPK modeling approach.

Clinical data supporting our results has been reported by Haegele et al. (2016) who showed that an ICG-R15 >20% on postoperative day 1 predicted poor postoperative outcome, which agrees well with the results shown in Figure 8B. The cutoff of ICG-R15 <20% allows to identify low-risk patients that are unlikely to have poor postoperative outcome after partial hepatectomy. This was confirmed by the high negative and low positive predictive values (>80 and 30%, respectively), suggesting that ICG-R15 is especially useful for the identification of low-risk patients. A recommendation was that subjects with ICG-R15 20–40% should undergo a more careful evaluation of the treatment options and additional information should be taken into consideration. In agreement with our results, Haegele et al. could show a significant improvement of the prediction of postoperative survival when using postoperative ICG-R15 (AUC_roc = 0.893) compared to preoperative ICG-R15 (AUC_roc = 0.719). As a side note, the experimental cutoff of >20% by Haegele et al. was determined in a Western population providing evidence that our classification approach could be generally applied despite developed on data from Japanese subjects.

Multiple studies showed only moderate performance in predicting post-hepatectomy outcome using preoperative ICG measurements. Gu et al. (2020) found that preoperative ICG-R15 achieved an AUC_roc of 0.657 (95% CI 0.576–0.739) and 0.640 (95% CI 0.445–0.836) for the prediction of post-hepatectomy liver failure and 90-day mortality, respectively. Wong et al. (2013) failed to achieve any significant prediction of postoperative severe morbidity using preoperative ICG-R15 (AUC_roc = 0.51, 95% CI 0.38–0.72). Wang et al. (2018) found that preoperative ICG-R15 surpassed both CTP-score and MELD for the prediction of severe post-hepatectomy liver failure, but with moderate AUC_roc = 0.724 (95% CI 0.654–0.787). In summary, this data provides a strong argument for our approach of using predicted postoperative ICG-R15 values via a PBPK model for predicting postoperative survival, which allowed to improve the discriminatory power (see Figure 8A).

The presented PBPK model predicts postoperative ICG-R15 immediately after partial hepatectomy, corresponding to clinical measurements on the first postoperative day (POD1). Using postoperative ICG-R15 on POD1 is supported by Haegele et al. (2016) who showed that an ICG clearance test was able to predict poor postoperative outcome as early as POD1.

A very interesting result is that our model accurately predicted intraoperative ICG measurements (during trial clamping) (Figure 7G), in which the blood flow to the area to be resected was clamped. These intraoperative ICG measurements are in very good agreement with post-operative ICG measurements (Thomas et al., 2015). Such an approach would allow to measure the expected postoperative ICG-R15 intraoperatively. Additional clamping data showed that the postoperative hospital stay was significantly longer for intraoperatively clamped ICG-R15 >20% (27.5 ± 14.1 days) compared to <20% (17.9 ± 9.2 days). Due to the good correlation between intraoperative clamped and postoperative ICG measurements this data provides additional support for the proposed lower cutoff of predicted postoperative ICG-R15 <20%. It is likely that postoperative ICG-R15 is not only a good predictor for survival but for postoperative complications in general as indicated by the data from Akita et al. (2008). Additional support comes from Haegele et al. which showed significantly reduced liver dysfunction (3.6 vs. 42.9%, p = 0.001), reduced severe morbidity (16.1 vs. 42.9%, p = 0.016) and hospitalization (7 vs. 11 days, p = 0.019) for POD1 ICG-R15 <20% vs. POD1 ICG-R15 >20% (Haegele et al., 2016).

Due to the high mortality rate, extended liver resection in the presence of cirrhosis is considered to be contraindicated. Recommendations are often that only selected patients with Child's A status or preoperative ICG-R15 of less than 10% undergo major hepatectomy (Kitano and Kim, 1997). As can be seen in Figure 8A (inset), even such a strict cutoff can still result in mortality after hepatectomy. We suggest using a combination of resection rate and individual liver damage (cirrhosis degree) estimated from preoperative ICG-R15 instead of relying on preoperative ICG-R15 alone. As shown by the PBPK2 classifier (and indirectly by the PBPK1 classifier) this approach allows a better evaluation of postoperative risk.

Overall, the clinical data shows large variability in ICG pharmacokinetic measurements, mostly due to inter-individual differences (e.g., Figure 7F). Possible explanations are differences in blood flow, plasma proteins or protein amount or activity of the ICG transporter. An important next step would be a systematic analysis of these possible causes of variability and account for these confounding factors in the model.

Importantly, due to the physiologically based modeling approach our predictions could easily be further individualized with the availability of respective data. A personalized risk prediction based on an individualized PBPK model could include (i) general information such as age, sex, and ethnicity; (ii) physiological information such as body weight, body fat percentage, cardiovascular parameters and organ volumes; (iii) information regarding the liver specifically such as liver perfusion, liver volume and quantification of ICG protein amounts as well as assessment of liver disease such as degree of cirrhosis would be of high relevance. Such an individualization could substantially improve the models prediction of postoperative liver function and outcome in patients undergoing partial hepatectomy. Going forward, an important next step will be to evaluate the model in the clinical context using a high quality dataset reporting individual ICG time courses in combination with a subset of the above-mentioned additional clinical data.