Prediction of Survival After Partial Hepatectomy Using a Physiologically Based Pharmacokinetic Model of Indocyanine Green Liver Function Tests

The evaluation of hepatic function and functional capacity of the liver are essential tasks in hepatology as well as in hepatobiliary surgery. Indocyanine green (ICG) is a widely applied test compound that is used in clinical routine to evaluate hepatic function. Important questions for the functional evaluation with ICG in the context of hepatectomy are how liver disease such as cirrhosis alters ICG elimination, and if postoperative survival can be predicted from preoperative ICG measurements. Within this work a physiologically based pharmacokinetic (PBPK) model of ICG was developed and applied to the prediction of the effects of a liver resection under various degrees of cirrhosis. For the parametrization of the computational model and validation of model predictions a database of ICG pharmacokinetic data was established. The model was applied (i) to study the effect of liver cirrhosis and liver resection on ICG pharmacokinetics; and (ii) to evaluate the model-based prediction of postoperative ICG-R15 (retention ratio 15 min after administration) as a measure for postoperative outcome. Key results are the accurate prediction of changes in ICG pharmacokinetics caused by liver cirrhosis and postoperative changes of ICG-elimination after liver resection, as validated with a wide range of data sets. Based on the PBPK model, individual survival after liver resection could be classified, demonstrating its potential value as a clinical tool.

The evaluation of hepatic function and functional capacity of the liver are essential tasks in hepatology as well as in hepatobiliary surgery. Indocyanine green (ICG) is a widely applied test compound that is used in clinical routine to evaluate hepatic function. Important questions for the functional evaluation with ICG in the context of hepatectomy are how liver disease such as cirrhosis alters ICG elimination, and if postoperative survival can be predicted from preoperative ICG measurements. Within this work a physiologically based pharmacokinetic (PBPK) model of ICG was developed and applied to the prediction of the effects of a liver resection under various degrees of cirrhosis. For the parametrization of the computational model and validation of model predictions a database of ICG pharmacokinetic data was established. The model was applied (i) to study the effect of liver cirrhosis and liver resection on ICG pharmacokinetics; and (ii) to evaluate the model-based prediction of postoperative ICG-R15 (retention ratio 15 min after administration) as a measure for postoperative outcome. Key results are the accurate prediction of changes in ICG pharmacokinetics caused by liver cirrhosis and postoperative changes of ICG-elimination after liver resection, as validated with a wide range of data sets. Based on the PBPK model, individual survival after liver resection could be classified, demonstrating its potential value as a clinical tool.

INTRODUCTION
Determining liver function is a crucial task in hepatology, e.g., for liver disease diagnosis or evaluation of pre-and postoperative functional capacity of the liver. An accurate assessment is especially relevant in the context of liver surgery as postoperative complications are associated with reduced functional capacity. A comprehensive characterization of the status of a patient and their liver is routinely performed before liver surgery such as partial hepatectomy. This includes among others anthropometric factors (e.g., age, sex, body weight), static liver function tests (e.g., ALT, AST, albumin, bilirubin, INR, prothrombin time), cardiovascular parameters (e.g., cardiac output, blood pressure, hepatic blood flow) and lifestyle factors (e.g., smoking, medication) as well as volumetric information obtained by CT or MRI scans. In addition, quantitative evaluation of liver function is often performed via pharmacokinetic measurements of test compounds specifically metabolized by the liver (dynamical liver function tests) such as methacetin [LiMAx (Rubin et al., 2017)] and MBT (Gorowska-Kowolik et al., 2017), caffeine (Renner et al., 1984), galactose (Bernstein et al., 1960) and indocyanine green (ICG) (Sakka, 2018).
ICG is an inert, anionic, water-soluble, tricarbocyanine dye that is bound to plasma proteins. After intravenous administration ICG is taken up exclusively by the liver and excreted unchanged into the bile. It is not reabsorbed by the intestine and does not undergo enterohepatic circulation (Wheeler et al., 1958). As a result, ICG is an ideal test compound to test hepatic uptake and biliary excretion. Based on the plasma time course of ICG, pharmacokinetic parameters are calculated as a proxy for liver function, the most common parameters being: (i) ICG retention ratio 15 min after administration (ICG-R15) [%]; (ii) ICG plasma disappearance rate (ICG-PDR) [%/min]; (iii) ICG-clearance [ml/min]; and (iv) ICG half-life (ICG-t 1/2 ) [min]. Reduced elimination of ICG by the liver is directly reflected by these parameters (Sakka, 2018).
Liver disease, especially advanced and more severe liver disease, is accompanied by a loss of liver function which can be quantified with dynamical liver function tests. The effects of liver disease on ICG-elimination have been studied extensively, e.g., in different stages of primary biliary cholangitis (PBC) (Vaubourdolle et al., 1991). ICG elimination is reduced in Gilbert's disease (Martin et al., 1976) as well as in patients with hepatic fibrosis and cirrhosis (Gadano et al., 1997). Interestingly, also non-liver diseases can affect ICG parameters, e.g., ICGclearance is significantly reduced in patients with chronic pancreatitis (Andersen et al., 1999).
Liver cirrhosis is an end stage liver disease and highly relevant in the context of hepatobiliary surgery. The most common causes for liver cirrhosis in Europe are alcohol abuse, chronic hepatitis B and/or C virus infection or non-alcoholic fatty liver disease (Hackl et al., 2016). The pathological characteristics of liver cirrhosis include degeneration of hepatocytes as well as a reduction of liver perfusion through increased portal resistance. An additional characteristics is the formation of intrahepatic shunts, which bypass part of the hepatic blood supply around the liver tissue. Consequently, no ICG can be extracted, resulting in further reduced elimination (Schuppan and Afdhal, 2008).
The severity of cirrhosis can be described using the Child-Turcotte-Pugh-Score (CTP) (Child and Turcotte, 1964;Pugh et al., 1973). The CTP is an empiric, qualitative, discontinuous classification of the severity of the "hepatic functional reserve" (Botero and Lucey, 2003). Based on a set of parameters, the CTP assigns a score from 5 to 15 to a cirrhotic patient, where the more severe a patients symptoms are the higher the score is. These parameters are serum bilirubin and serum albumin concentrations, the International Normalized Ratio (INR), amount of ascites in ultrasound and the degree of encephalopathy (Child and Turcotte, 1964;Pugh et al., 1973). Patients are classified depending on their mortality risk as CTP-A (5-6 points, low risk), CTP-B (7-9 points, intermediate risk) or CTP-C (10-15 points, high risk). Differences in ICG-elimination between cirrhotic patients and control subjects has been widely assessed (Caesar et al., 1961;Gilmore et al., 1982;Burns et al., 1991;Figg et al., 1995;Mukherjee et al., 2006;Pind et al., 2016;Møller et al., 2019). A good correlation between CTP-score and ICG-elimination has been reported with ICG-elimination decreasing as CTP-score increases (Figg et al., 1995;Mukherjee et al., 2006;Pind et al., 2016;Møller et al., 2019).
For many primary and secondary liver tumors, partial hepatectomy is the only curative treatment option. Liver resection (partial hepatectomy) describes the surgical removal of a part of the liver. Hepatectomy is an important procedure in general surgery with more than 20,000 liver resections in Germany per year (Filmann et al., 2019). The procedure has been widely performed for the treatment of various liver diseases, such as malignant tumors, benign tumors, calculi in the intrahepatic ducts, hydatid disease, and abscesses (Jin et al., 2013). Despite advances in technology and increasing numbers of highly experienced and specialized clinicians in the field of hepatobiliary surgery, postoperative morbidity and mortality is still a major issue as complex resections are increasingly being performed in older and higher risk patients (Jin et al., 2013). Major hepatectomy in the presence of cirrhosis is considered to be contraindicated due to the high mortality rate. Only selected patients with stable Child A status or an ICG-R15 of less than 10% may potentially be considered for extended liver resection (Kitano and Kim, 1997).
A key challenge in HPB (hepato pancreato biliary) surgery is to predict functional capacity of the future liver remnant and thus reduce morbidity and mortality after extended liver resection. As a result, the decision whether or not a partial hepatectomy can be performed is often based on predictions of postoperative liver function (in addition to the remnant liver volume), which are in turn based on preoperative assessment of liver function (and volume). Understanding how cirrhosis alters liver function as measured via ICG is of high clinical relevance. Elucidating how ICG parameters change with increasing CTP score would be a valuable asset for the functional evaluation of patients with liver disease. Important questions for the evaluation of liver function with ICG in the context of HPB surgery are (i) how liver disease, especially cirrhosis, alter ICG elimination, and (ii) if postoperative survival can be predicted from preoperative ICG measurements. Within this work a physiologically based computational model of ICG pharmacokinetics was developed and applied to study these questions.

Data
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)

Model Implementation
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).

Model Parameterization
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 r i,k = (y i,k − m i,k ( p)) is the residual of time point i in time course k for model prediction m i,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.

Uncertainty Analysis
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).

Indocyanine Green Pharmacokinetic Parameters
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: 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
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 Csupport 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.

RESULTS
Within this work a PBPK model of ICG pharmacokinetics was developed and applied to study (i) how liver disease, especially cirrhosis, alter ICG elimination, and (ii) if postoperative survival can be predicted from preoperative ICG measurements in the context of partial hepatectomy. 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

LI__f_oatp1b3
Scaling factor of transport protein amount 1 -respective clinical protocols is provided in Table 1. All data is freely available from https://pk-db.com.

Model of ICG Pharmacokinetics
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.

Distribution and Blood Flow
The distribution of ICG on the whole-body level is modeled using a network of blood flows representing the systemic circulation. From the venous blood, ICG is transported through the lung into the arterial blood from where it can reach the liver on two paths: (i) through the hepatic artery and (ii) through the gastrointestinal tract via the portal vein. Because the liver is the only tissue partaking in the elimination of ICG, all other organs (e.g., kidney, adipose tissue, muscle, etc.) were pooled into the rest compartment. Each organ consists of a blood compartment (representing the vessels) and a tissue compartment. ICG transport via blood flow was implemented as irreversible transport. The transport v i from compartment i to the next compartment is determined by the ICG concentration C i in compartment i and a compartment-specific blood flow Q i . Q i is determined by the cardiac output Q CO and a compartment specific fractional tissue blood flow fQ i . Multiple conservation conditions hold in the model to ensure mass and flow balance. First the sum of blood flows from the arterial to the venous compartment must equal the sum of flows in the opposite direction: Q CO = Q lu = Q h + Q re . Flow into an organ must be equal to the flow out of the organ. For example, hepatic venous blood flow must be equal to the sum of hepatic arterial and portal venous blood flow: Q h = Q ha + Q po .

Hepatic Metabolism and Biliary Excretion
The liver model ( Figure 1B) consists of three consecutive transport reactions of ICG. After ICG is taken up in the liver it is excreted into the bile. Both transport reactions are modeled as irreversible Michaelis-Menten-kinetics. From the bile, ICG is transported into the feces modeled via a first order kinetic. All transport kinetics scale with the liver volume V li .

Parameter Fitting
Parameter fitting of the model was performed using a subset of ICG time courses and extraction-ratio measurements (see Table 1). No ICG pharmacokinetic parameters were used in the model fitting. Overall, 5 model parameters were fitted (see Table 2). Two of them determine the import of ICG in the liver, three determine the subsequent excretion in the bile. The agreement between fit data and model predictions improved substantially during parameter fitting and all training data with the exception of three simulations (see biliary excretion in Figure 2H blue and green curves; and the second, slower elimination phase in Figure 6A) could be described very well after parameter fitting.

Modeling Liver Cirrhosis
The reference model, representing a healthy human subject, was adjusted to simulate cirrhosis by including a combination of functional tissue loss (due to scarring and necrosis in cirrhosis) and the formation of intrahepatic shunts, both key hallmarks of cirrhosis ( Figure 1C). The loss of functional liver tissue was controlled via the parameter f tissue_loss ∈ [0, 1) which defines the fraction of parenchymal cell volume lost in the liver due to the disease. For modeling arteriohepatic and portosystemic shunts two additional blood vessels were introduced into the model. They connect the hepatic artery and the portal vein directly to the hepatic vein. As a result, a part of the portal venous and arterial blood bypasses the active liver tissue and is shunted to the hepatic venous blood compartment, so that ICG can not be extracted (corresponding to in silico shunts). The amount of blood that flows through the shunts is determined by the parameter f shunts ∈ [0, 1), which defines the fraction of blood bypassing the liver. The remaining blood (1 − f shunts ) reaches the liver tissue and ICG can be extracted. To simulate various degrees of cirrhosis the parameters f shunts and f tissue_loss were varied in lockstep by coupling them into the parameter f cirrhosis .

Modeling Hepatectomy
The developed model allows to predict changes in ICG pharmacokinetic parameters after partial hepatectomy ( Figure 1D). In silico liver resections were simulated by reducing the fractional liver volume FVli by up to 90% (corresponding to a resection rate of 90%). The absolute liver volume is determined with the body weight BW via FVli · BW. All liver resections were simulated under varying degrees of cirrhosis as described in section 3.1.4.

Model Validation for Healthy Controls
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 ICGclearance 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(Martin et al., , 1976Meijer et al., 1988) (Figures 3E-G).
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 ICGclearance 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 ICGt 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  (Martin et al., 1975(Martin et al., , 1976Meijer et al., 1988).

FIGURE 4 | Dependency of ICG pharmacokinetic parameters on cirrhosis. Simulations of 4 specific cirrhosis degrees are indicated by stars. (A-D)
Dependency of ICG pharmacokinetic parameters on the degree of shunting (green), degree of tissue loss (yellow) and degree of cirrhosis (black-blue-red). (E) Correlation between ICG-R20 and ICG-t 1/2 in cirrhotic and control subjects (Cherrick et al., 1960). (F) Correlation between ICG-R20 and ICG-kel in cirrhotic and control subjects (Caesar et al., 1961). (G) Correlation between ICG-PDR after an ICG dose of 0.5 and 5.0 mg/kg in control subjects and subjects with various liver diseases (Leevy et al., 1967). (H) Correlation between ICG-clearance after a bolus administration and during a constant infusion of ICG in cirrhotic subjects (Burns et al., 1991). multiple datasets in which ICG pharmacokinetic parameters were reported in patient subgroups of different CTP-Scores (Figg et al., 1995;Møller et al., 1998Møller et al., , 2019Herold 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  Venous concentration after a bolus ICG administration in a healthy subject and a cirrhotic patient (f cirrhosis was set to 0.7 corresponding to moderate cirrhosis) (Burns et al., 1991). (C) Venous concentration during a constant ICG infusion in healthy and cirrhotic subjects (Caesar et al., 1961). (D,E) Hepatic venous and arterial ICG concentration and ICG extraction ratio in a cirrhotic patient (f cirrhosis was set to 0.54 corresponding to mild-moderate cirrhosis) (Keiding et al., 1993). (F-H) ICG extraction ratio in cirrhotic and control subjects during a constant ICG infusion (Caesar et al., 1961;Leevy et al., 1962;Gadano et al., 1997). 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(Møller et al., , 2019Herold 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.

Model Validation for Hepatectomy
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.  (Stockmann et al., 2009;Thomas et al., 2015). (F) Correlation between the measured postoperative ICG-kel and the estimated postoperative ICG-kel (product of preoperative ICG-kel and the future liver remnant) (Okochi et al., 2002;Ohwada et al., 2006;Sunagawa et al., 2021). (G) Correlation between postoperative ICG-PDR and intraoperative ICG-PDR during trial clamping (Thomas et al., 2015). (H) Correlation between postoperative and preoperative ICG-PDR (Thomas et al., 2015). Simulations for (E-H) were performed for varying resection rates in healthy subjects and three different degrees of cirrhosis. In summary, the model allows to systematically predict the changes of ICG pharmacokinetic parameters in HPB surgery under various degrees of liver disease (cirrhosis).

Prediction of Post-hepatectomy Survival
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 databased approaches such as preoperative ICG-R15 and calculated postoperative ICG-R15.

Model Building and Validation
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 hepatopancreatico-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.

Model Assumptions and Limitations
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.
FIGURE 8 | Classification of survivors/non-survivors after hepatectomy. All classification models are set up to predict non-survival rather than survival. Classification data (n = 141) from Seyama and Kokudo (2009) and Wakabayashi et al. (2004). (A) PBPK model based prediction of postoperative ICG-R15 in survivors and non-survivors. f cirrhosis was estimated for every subject based on preoperative ICG-R15 and post-operative ICG-R15 calculated using the respective resection rate (used in PBPK1). Corresponding measured preoperative ICG-R15 as inset (used in Data1A). (B) Dependency of the evaluation metrics of the cross-validated classification model PBPK1 on the predicted postoperative ICG-R15 cutoff (mean ± SD). PPV, positive predictive value; NPV, negative predictive value. (C) ROC curve for the prediction of non-survival after hepatectomy using the complete dataset with cross-validation (mean ± SD) for classification models Data1A, Data1B, and PBPK1. (D) Predicted postoperative ICG-R15 and survival status depending on the resection rate and f cirrhosis (used in PBPK2). f cirrhosis was estimated for every subject individually based on preoperative ICG-R15 and post-operative ICG-R15 calculated using the respective resection rate. Shaded blue areas correspond to respective predicted postoperative ICG-R15. (E) Decision boundary of the two-dimensional classification model PBPK2 based on the resection rate and f cirrhosis using the complete dataset. White area: predicted survivor; blue area: predicted non-survivor. (F) ROC curve for the prediction of non-survival after hepatectomy using the complete dataset with cross-validation (mean ± SD) for classification models Data2 and PBPK2.
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.

Prediction of Survival
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.

Previous Research
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).

Future Orientations
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 abovementioned additional clinical data.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://pk-db.com.

AUTHOR CONTRIBUTIONS
AK and MK designed the study, developed the computational model, built the classifiers, implemented and performed the analysis, and wrote the initial draft of the manuscript. JG provided support with PK-DB and data curation. All authors discussed the results and revised the manuscript critically.