Physiologically Based Pharmacokinetic Modelling and Simulation to Predict the Plasma Concentration Profile of Doxorubicin

Doxorubicin (DOX) is still an important anticancer agent despite its tricky pharmacokinetics (PK) and toxicity potential. The advent of systems pharmacology enables the construction of PK models able to predict the concentration profiles of drugs and shed light on the underlying mechanisms involved in PK and pharmacodynamics (PD). By utilizing existing published data and by analysing two clinical case studies we attempt to create physiologically based pharmacokinetic (PBPK) models for DOX using widely accepted methodologies. Based on two different approaches on three different key points we derived eight plausible models. The validation of the models provides evidence that is all performing as designed and opens the way for further exploitation by integrating metabolites and pharmacogenomic information.


Introduction
Despite being an "old" drug, Doxorubicin (DOX) still remains an important and valuable therapeutic agent in cancer therapy [1]. Its clinical use, however, is limited due to safety issues. The latter is correlated to the cumulative dose used and is manifested mainly as cardiotoxicity that could potentially lead to congestive heart failure (CHF) [2] causing up to 50% mortality [3]. There are two proposed mechanisms that are correlated to the plasma levels of its two metabolites doxorubicinol (DOXol) and doxorubicinone (DOXone), although there is still a controversy, due to the pleiotropic pharmacological effects of DOX at the molecular level and complex pharmacokinetic (PK) profile [4][5][6]. Both metabolites produce reactive oxygen species (ROSs) intracellularly that trigger cytotoxicity and programmed cell death [4][5][6]. Of special research interest is also the fact that the clinical outcome shown by DOX exhibits a sex-related behaviour, an observation that remains elusive and challenges the scientific community [7].
The pegylated liposomal DOX formulation (Doxil ® , Caelyx ® ) that entered the market in 1995 represents the first FDA-approved nanomedicine that successfully addressed the problematic PK behaviour and specific safety profiles of DOX [8,9]. Specifically, it showed enhanced circulation time of DOX as well as a better tumour accumulation profile [10,11]. Recently, quantitative systems pharmacology (QSP) is in the spotlight of modern pharmacology by providing, through physiology-based pharmacokinetic (PBPK) modelling, mechanistic insights on the PK processes towards achieving better efficacy and safety profiles in the clinical setting [12]. PBPK models are based on in vitro-in vivo correlation (IVIVC) procedures. Striving to be as mechanistic as possible in nature, they are based on the underlying anatomical, physiological, and biochemical characteristics of an organism [13]. In such models, the body is a multicompartment system with every compartment representing a different organ connected to other compartments by blood or lymph circulation, through a system of differential equations describing different phenomena, such as blood flow, cardiac output, organ volumes, glomerular filtration rate etc. [14].
Such capability gives the advantage to PBPK models and leads to better identification of the sources of PK variability allowing to extrapolate to different subpopulations [10]. In this context, precision medicine could be achieved in the clinical setting by connecting PBPK models with pharmacodynamic (PD) prediction models and their capacity for population simulation [i.e., the prediction of the effects of age, gender, comorbidities, genetic polymorphisms, lifestyle factors (e.g., smoking) and more] [13].
Finally, an advantage of such models is their ability to simulate the drug concentration profile on the site of action (e.g., targeted organ or tissue), allowing the refinement of dosage schemes and the achievement of maximum safety and effectiveness profiles [15]. This can be of great value in the case of advanced formulations, such as nanoformulations, as the combination of knowledge about the structure and the function of target organs, with the physicochemical properties of the nanocarriers, the individual parameters of each patient and the drug properties could create favourable conditions for individualized treatment [16].

Materials and Methods
To address the issues related to therapeutic peculiarities related to DOX and its metabolites in the body we attempt the development of a PBPK model capable of predicting the PK profile of DOX in the plasma. This first step presented in our work, provides a solid basis for further incorporation of DOX main metabolites [Doxorubicinol (DOXol) and Doxorubicinone (DOXone)] kinetics in the future. With a working model at hand, adjustments will be possible so that the potential application of such a PBPK model in various innovative DOX nano-formulations, e.g., Doxil ® , could be a useful clinical tool in treating cancer patients. The ability of the clinicians to estimate, through such a PBPK model, the proper dose of DOX to achieve the maximum efficacy and safety profiles in individual patients is fundamental within the precision medicine concept.
There are several attempts in the literature trying to develop PBPK models for DOX; however, with different end goals. Dubbelboer et al. created both PBPK and semi-PBPK models for DOX attempting to incorporate intracellular binding as a distribution factor in the latter one [17]. Gustafson et al. created a mouse PBPK DOX model with macromoleculespecific binding as the main factor for distribution and organ-specific metabolism and excretion. These authors believe their model has the potential to predict the magnitude of PK interactions of DOX with other drugs, as well as more efficiently addressing various clinical situations [18]. Hanke et al. developed a PBPK model for a DOX fusion molecule, namely "zoptarelin doxorubicin". To this end, they initially developed a DOX PBPK model (as DOX is the metabolite of the fusion molecule). They also utilized DNA intracellular binding to predict distribution. They had unspecified hepatic clearance and bile elimination [19]. He et al. created a PBPK model to assess DOX disposition at many levels (system, tissue interstitial, cell and subcellular organelles) by analyzing mice data and scaling up for humans thus gaining insights concerning toxicity [20]. In our approach, various relevant models were created and validated using a middle-out approach. In this approach, existing clinical observations are utilized in a reverse translation approach to combine any prior information on the drug or system into the analysis of the clinical observations to project forward beyond the scope of the initial observations. In such methodologies, the models comprise of three different interacting components, the system data, the drug data, and the clinical trial data. System data refers to the properties of the organism for which the model is created (e.g., blood flow or enzyme activity). Drug data refers to the physicochemical data and PK information of the modelled drug (e.g., drug pK a , affinity for enzymes, etc.). Clinical trial data refers to the administration settings (e.g., population age, female percentage, dosage scheme, etc.) [21]. The interaction of these three components is summarized in Figure 1.
organism for which the model is created (e.g., blood flow or enzyme activity). Drug data refers to the physicochemical data and PK information of the modelled drug (e.g., drug pKa, affinity for enzymes, etc.). Clinical trial data refers to the administration settings (e.g., population age, female percentage, dosage scheme, etc.) [21]. The interaction of these three components is summarized in Figure 1.

Clinical Studies Used
To construct a model two types of datasets are needed: (a) a training dataset, utilized in the development of the model and (b) a validation dataset, used for independent validation of the model. For our models, we used the clinical study of Camaggi et al. as the training dataset [22] and the clinical study of Speth et al. as the independent validation dataset [23]. The details of those clinical studies are summarized in Table 1.

Software
The Simcyp (version 19R1) simulator (Simcyp Ltd., Sheffield, UK) was used to simulate the PK profile of DOX. The criteria for assessing the predictive performance of the models were the predicted/observed ratio for the AUC and the Cmax of DOX. A perspective on the qualification and verification of PBPK models, as well as examples of regulatory PBPK submissions, can be found in the work of Shebley et al. [24].

Virtual Population Characteristics (System Data)
As both clinical studies utilized in this work (training and validation dataset) refer to a cancer population, we opted for the cancer population (Sim-Cancer) that is included within the Simcyp simulator. For this special population, many adjustments have been made to better account for the specific changes that are expected to be found in the physiological parameters of such a population.

Clinical Studies Used
To construct a model two types of datasets are needed: (a) a training dataset, utilized in the development of the model and (b) a validation dataset, used for independent validation of the model. For our models, we used the clinical study of Camaggi et al. as the training dataset [22] and the clinical study of Speth et al. as the independent validation dataset [23]. The details of those clinical studies are summarized in Table 1.

Software
The Simcyp (version 19R1) simulator (Simcyp Ltd., Sheffield, UK) was used to simulate the PK profile of DOX. The criteria for assessing the predictive performance of the models were the predicted/observed ratio for the AUC and the C max of DOX. A perspective on the qualification and verification of PBPK models, as well as examples of regulatory PBPK submissions, can be found in the work of Shebley et al. [24].

Virtual Population Characteristics (System Data)
As both clinical studies utilized in this work (training and validation dataset) refer to a cancer population, we opted for the cancer population (Sim-Cancer) that is included within the Simcyp simulator. For this special population, many adjustments have been made to better account for the specific changes that are expected to be found in the physiological parameters of such a population.

Development of DOX PBPK Model (Drug Data)
To construct the models, we need to input the physicochemical properties of DOX, values for DOX renal, metabolic, and biliary clearance and to select an appropriate model for DOX distribution. The physicochemical properties of DOX were either calculated, found by literature review or online database utilization. The pharmacokinetic parameters were calculated based on the data from the training dataset. For the calculations, two patients of the original case study were excluded from further analysis: (a) Patient 69 because of hepatic metastases, extrahepatic obstruction, and percutaneous biliary drainage and (b) Patient 72 as the hepatic clearance exceeded the hepatic blood flow that was calculated. The calculated  Table 2. P.E. tool: parameter estimation tool (i.e., a Simcyp simulator tool); HEP: intrinsic metabolic clearance calculated per 10 6 hepatocytes; HLM: intrinsic metabolic clearance calculated per mg of microsomal protein; HLC: intrinsic metabolic clearance calculated per mg of cytosolic protein; V ss : Volume of distribution in steady state; V sac : volume of single adjusting compartment (see Supplementary Materials Section S12 for details); Q sac : single adjusting compartment blood flow (see Supplementary Materials Section S12 for details); K p Scalar: A scaling value for the calculated K p values of each tissue in a full PBPK model (see Supplementary Materials Section S12 for details).

Calculating DOX Renal Clearance
The Simcyp simulator requires renal clearance to be inputted as the renal clearance of a healthy 20-30 y.o. male (refCL R ). We calculated the refCL R based on the patient's GFR and the GFR values expected for 20-30 y.o. males as shown in Equation (1).
where refCL R,i is the reference renal clearance for each patient i, CL R,i is the renal clearance of each patient i, GFR ref is the reference GFR for a 20-30 y.o. male and eGFR i is the expected GFR for each patient i.
As GFR values per patient were not mentioned in the work of Camaggi et al. used as a training dataset [22], we calculated those values based on two methods: Method A and Method B.

Method A
In this method, we used the approach of Davies and Shock [25] for calculating the expected GFR (eGFR) for each patient i by utilizing expected GFR values per age group (Equation (2)).
where eGFR i is the expected GFR for each patient i, GFR g,i is the expected GFR for the group that each patient i belonged based on age, and BSA i is the is body surface area of each patient i, a value that was recorded in the training dataset clinical study.

Method B
In this method, we used the approach of Wright et al. [26] and calculated the expected GFR values based on the expected serum creatinine levels of each patient i (Equation (3)).
where eGFR i is the expected GFR, age i is the age and SCr i is the serum creatinine levels of each patient i. The parameter SEX i takes two distinct values, 0 for males and 1 for females based on the gender of each patient i. However, as the gender of each patient i was not included in the study, an average GFR value was calculated by averaging the GFR values assuming both male and female gender (Equation (4)).
where eGFR i,m is the estimated GFR values for patient i if considered male and eGFR i,f if considered female.
The calculated values of renal clearance of DOX for a 20-30 y.o. healthy male (refCL R ) are summarized in Table 2. For the values of each individual patient for each method, see Supplementary Materials Section S3. It is noteworthy that the calculated renal clearance for DOX exceeds the calculated GFR values, thus implicating the involvement of active processes in renal excretion, a fact that is discussed below in the model limitations.

Calculating DOX Hepatic Clearance
The Simcyp simulator can calculate the hepatic clearance of a drug by scaling using different in vitro systems based on the following method: Hepatic metabolic clearance can be predicted using either hepatocyte, cytosolic fraction, or microsomal fraction in vitro systems. Biliary excretion can be calculated using the hepatocytes in vitro system. Thus, to be able to predict the hepatic clearance for a simulated patient, we must calculate the hepatic intrinsic clearance first and then we must calculate the intrinsic metabolic clearance and the intrinsic biliary excretion and correct them by the appropriate scaling factors as seen in Table 3.
To calculate the hepatic intrinsic clearance, we can use Equation (5) which is based on the well-stirred liver model.
where CL H,b,i is the hepatic blood clearance for each patient i, Q H,i is the hepatic blood flow of each patient i, fu b is the unbound fraction of the drug in blood and CLu int,H,b,i is the intrinsic hepatic clearance of the unbound drug for each patient i. By rearranging Equation (5), we obtain Equation (6) and thus the intrinsic clearance of the unbound drug can be calculated.
From the above, it is obvious that there are three parameters that need to be calculated in order to calculate the hepatic intrinsic unbound clearance of a drug, namely the (1) hepatic blood flow of each patient i, Q H,i, , (2) the unbound fraction of a drug in blood, f u,b and (3) the hepatic blood clearance for each patient i, CL H,b,i .

Calculating the Hepatic Blood Flow for Each Patient i
Blood is supplied to the liver via two paths: (a) through the hepatic artery and (b) through the portal vein. In the virtual population we selected for our models (i.e., cancer population), the liver receives a predefined percentage of the cardiac output ( f CO L ) via the above-mentioned methods, different for each gender. As the clinical study used as the training dataset did not include gender, the average percentage for the two genders was used, as shown in Table 4. In the cancer population used for the models, the formula for calculating the cardiac output for each patient i is a function of body surface area (BSA) and age as shown in Equation (7).
where CO i is the cardiac output of each patient i, BSA i is the body surface area of each patient i and age i is the age of each patient i. Thus, the hepatic blood flow for each patient i (Q H,i ) can be calculated using Equation (8).
For the values of the cardiac output calculated for each individual patient, see Supplementary Materials Section S4. f CO L represents the percentage of CO for the liver and is calculated in Table 4.

Calculating Unbound Fraction in Blood
We know that for the unbound fraction of a drug in blood, Equation (9) applies.
As we mentioned earlier, the unbound fraction of DOX in blood was 0.25 and the B:P ratio (R B:P ) was calculated to be 0.87 (see Supplementary Materials Section S5). Thus, for DOX, f u,b was calculated to be 0.2175 as shown here:

Calculating Hepatic Blood Clearance
Hepatic blood clearance (CL H,b ) can be calculated for each patient i by utilizing Equation (10): where CL H,i is the hepatic plasma clearance for each patient i and R B:P is the ratio of DOX concentration in blood vs. plasma. Thus CL H,B was calculated for each patient i. For the values of each individual patient, see Supplementary Materials Section S6.

Calculating Intrinsic Hepatic Clearance
Using the values for the hepatic blood flow for each patient i, the unbound fraction of DOX in blood and the values of hepatic blood clearance calculated for each patient i in Equation (6), we calculated CLu int,H,b for each patient i.

Separating Hepatic Clearance to Hepatic Metabolic Clearance and Biliary Excretion
Based on the literature [27], approximately 40% of DOX is excreted in the bile as unchanged drug while 5-12% of DOX and its metabolites are excreted in urine. In the training datasets, the patient had an average plasma clearance (CL P ) of 51.75 L/h, an average renal clearance (CL R ) of 5.59 L/h, and thus an average fraction excreted in urine (f e ) of 11.07%. (For detailed calculations see Supplementary Materials Section S7). Based on the above, the fractions that are eliminated via different routes used for the construction of the model are presented in Table 5.

Way of Elimination Percentage
Thus, the biliary excretion can be calculated to be 44.94% of the hepatic clearance ( f CL,H,bile ) by using Equation (11) and the values from Table 5.
where f bile is the percentage of DOX excreted in bile and f e is the percentage of DOX excreted in urine as unchanged drug. By assuming that the ratio of biliary clearance to hepatic clearance is the same as intrinsic biliary clearance to intrinsic hepatic clearance, we can calculate the intrinsic biliary clearance and intrinsic metabolic clearance by the intrinsic hepatic clearance we calculated before for each patient i. For the values of each patient, see Supplementary Materials Section S8.

Correction of Intrinsic Biliary Clearance for In Vitro System Scaling Factors
Intrinsic biliary clearance must be corrected for Scaling Factor 2 (i.e., liver weight-LW) and Scaling Factor 1 (in this case, hepatocellularity per gram liver-HPGL). The correction is carried out using Equation (12).
where CLu int(Bile),i the biliary excretion per million hepatocytes for each patient i, CLu int,bil,b,i the intrinsic biliary excretion for each patient i and HPGL i the hepatocellularity per gram liver for each patient i. The individual values calculated and information about calculating HPGL for each patient can be found in Supplementary Materials Section S9. The average value of unbound intrinsic metabolic clearance per million hepatocytes is summarized in Table 2.

Calculating Intrinsic Metabolic Clearance
As was the case with intrinsic biliary excretion, again in the case of intrinsic metabolic clearance, we must correct the values for both Scaling Factor 2 (i.e., LW) and Scaling Factor 1 (i.e., HPGL, MPPGL or CPPGL). Here, one can follow two approaches: either assign all metabolic clearance to be predicted using HPGL (Method C) or by considering the metabolic pathway of DOX, try to distribute it to different in vitro systems based on the location of the actual metabolizing enzymes (Method D). The latter will be useful as it can be utilised to simulate the formation of the different DOX metabolites in the future.
The following equations are used for the calculation of each in vitro system: CLu int(Met),HLC,i = CLu int,met,b,i LW × CPPGL i The values for correcting the entire intrinsic metabolic clearance for the three in vitro systems are summarized in Table 2. Values for MPPGL and HPGL come from the literature [28,29] while the Simcyp simulator calculates values for CPPGL. Details of calculating individual HPGL, MPPGL and CPPGL values as well as the values for each patient are in Supplementary Materials Section S9. Individual patient values for biliary excretion are in Supplementary Materials Section S10. The corrected values for intrinsic metabolic clearance for each patient are in Supplementary Materials Section S11. Table 6 shows the final values of intrinsic metabolic clearance for each in vitro system when considering the expected relative contribution of each based on the DOX metabolic pathway. This distribution was performed by considering the following facts [27]:

•
The enzymes that participate in the primary metabolic pathway, which is the twoelectron reduction, mainly aldoketoreductases (AKRs) and carbonylreductases (CBRs) are located in the cytoplasm and are expected to be found in the cytoplasmic fraction after centrifugation.

•
The enzymes that participate in the secondary metabolic pathway, which is the oneelectron reduction, are located mainly in mitochondria and sarcoplasmic reticulum and are expected to be found in the microsomal fraction after centrifugation.

•
The enzymes that participate in the minor metabolic pathway, which is the deglycosylation, are not specified, contribute only 1-2% of the total metabolism and thus their contribution can be attributed per 10 6 hepatocytes.

Selecting a Distribution Model for DOX
To simulate DOX distribution, either the minimal PBPK (mPBPK) model or the full PBPK (fPBPK) models can be used. Both are models provided by the Simcyp simulator to predict the distribution of drugs. In the case of the mPBPK, the V ss value for DOX was user-inputted based on the values of the training dataset. In the case of the fPBPK model, V ss was calculated in silico. The V ss values are summarized in Table 2. Thus, there are two methods for calculating DOX distribution, either using the mPBPK model (Method E) or the fPBPK model (Method F).
By using the Simcyp parameter estimation (PE) tool, the appropriate parameters were calculated (i.e., the V sac volume (L) and V sac blood flow (L/h) in the case of mPBPK and the K p scalar in the case of fPBPK). The calculated values are summarized in Table 7. For more information on distribution models utilized by Simcyp, see Supplementary Materials Section S12.

Generating DOX Models
Summarizing all the above, we have two different approaches on three key calculations thus resulting in eight possible DOX models. The models are presented in Table 7.  Table 6).

Virtual Patient Demographics for the Development of the Model Based on the Training Dataset (Clinical Settings Data)
The age of the virtual patients was set from 42 to 72 years to match the age of the patients in the training dataset. The proportion of females to males was set to 0.5 as there was no mention of gender in the clinical study. All simulations run for 10 groups of 10 persons each to study the population PK. The virtual study ran for 168 h. Patients received 60 mg/m 2 of DOX via IV bolus infusion over 2 min at time 0 in accordance with the conditions of the clinical study. Finally, we opted to record 10,000 virtual plasma samples over the course of 168 h. Figure 2A represents the predicted mean, 95th percentile and 5th percentile concentration versus time course of DOX for model 8, presented as a sample, for the virtual patient population. Figure 2B shows the mean concentration vs. time values for all 8 models. For the plots of the remaining seven models, see Supplementary Materials Section S13. patients in the training dataset. The proportion of females to males was set to 0.5 as there was no mention of gender in the clinical study. All simulations run for 10 groups of 10 persons each to study the population PK. The virtual study ran for 168 h. Patients received 60 mg/m 2 of DOX via IV bolus infusion over 2 min at time 0 in accordance with the conditions of the clinical study. Finally, we opted to record 10,000 virtual plasma samples over the course of 168 h. Figure 2A represents the predicted mean, 95th percentile and 5th percentile concentration versus time course of DOX for model 8, presented as a sample, for the virtual patient population. Figure 2B shows the mean concentration vs. time values for all 8 models. For the plots of the remaining seven models, see Supplementary Materials Section S13.  Table 8 summarizes the AUC and Cmax for the different models. The observed mean value for was reported by the authors in their clinical study and was calculated  Table 8 summarizes the AUC and C max for the different models. The observed mean value for AUC 0-168 was reported by the authors in their clinical study and was calculated after the exclusion by us of two patients (see Section 2.4). The predicted population mean AUC 0-168 for each model was calculated by the Simcyp simulator. The observed mean C max value was calculated for time 0 h based on the triexponential equations the authors of the clinical study suggested for DOX resulting from their measurements. The predicted mean C max of the population for each model was calculated by the simulator.  Table 6). For the origin of observed values, see Section 2.7.

Patient Demographics for the Development of the Model Based on the Validation Dataset
As mentioned earlier, the clinical study of Speth et al. was used as the validation dataset. In their work, Speth et al. studied the pharmacokinetic behaviour of DOX, by analysing its concentration both in plasma and at a cellular level [23]. Eighteen patients with leukaemia participated in the study. Their age ranged from 17 to 67 y.o. and there were 8 females and 10 males. All patients had normal renal and hepatic function.

Administration and Sample Retrieval
The patients were administered Vincristine on day 2 (dose of 1 mg/m 2 ), and Cytarabine each day for days 1 to 7 (dose 200 mg/m 2 ). In sixteen patients, DOX was administered on days 1, 2 and 3 (dose of 30 mg/m 2 ) and seven patients received DOX as an IV bolus injection, five as a 4 h IV infusion, four as an 8 h IV infusion. In two patients, DOX was administered on day 1 as a 72 h IV infusion. Blood samples were taken from 5 to 240 min after administration from at least 2 patients per therapeutic scheme. For the rest of the patients in each therapeutic scheme, blood samples were taken when DOX is expected to reach its maximum and minimum concentration. After centrifuging for 10 min, plasma was kept at −20 • C. Two of the therapeutic schemes were selected for independent validation of the models. Specifically, the IV bolus and the 8 h IV infusion of DOX.

Analytical Method and Pharmacokinetic Analysis
The samples were analysed via high-performance liquid chromatography (HPLC). The sensitivity of the analytical methodology was 1 ng/mL. DOX plasma concentrations were described by biexponential equations. C max values (at 5 min after the third administration) and AUC for each therapeutic scheme are summarized in Table 9.

Clinical Settings of Virtual Patients for Validation Dataset
Again, the simulated cancer population was selected for the simulation, as it better represents the clinical study population of Speth et al., used as a validation dataset. As the Simcyp simulator requires a minimum age of 20 for this specific population, the age of the virtual patients was set from 20 to 67 years. The female analogy was set at 44.5% in accordance with the study. Ten groups of 10 patients each were simulated. The study duration was set at 168 h. The dose was set at 30 mg/m 2 , based on the selected therapeutics schemes. Figure 3A represents the predicted mean, 95th percentile and 5th percentile concentration versus time course of DOX for model 8, presented as a sample, for the virtual patient population for the IV bolus injection trial (administration of 30 mg/m 2 on days 1, 2 and 3). Figure 3B shows the mean concentration vs. time values for all 8 models. For the plots of the remaining seven models, see Supplementary Materials Section S14.  Figure 4A represents the predicted mean, 95th percentile and 5th percentile concentration versus time course of DOX for model 8, presented as a sample, for the virtual patient population for the DOX IV infusion trial (administration over 8 h of 30 mg/m 2 on days 1, 2 and 3). Figure 4B shows the mean concentration vs. time values for all 8 models.  Figure 4A represents the predicted mean, 95th percentile and 5th percentile concentration versus time course of DOX for model 8, presented as a sample, for the virtual patient population for the DOX IV infusion trial (administration over 8 h of 30 mg/m 2 on days 1, 2 and 3). Figure 4B shows the mean concentration vs. time values for all 8 models. For the plots of the remaining seven models, see Supplementary Materials Section S15.

Observed Concentration Values for the Validation Dataset
Tables 10 and 11 summarize the AUC and Cmax values (observed vs. predicted) for the eight different models for the two validation trials (IV bolus and IV infusion accordingly).

Observed Concentration Values for the Validation Dataset
Tables 10 and 11 summarize the AUC and C max values (observed vs. predicted) for the eight different models for the two validation trials (IV bolus and IV infusion accordingly).
The observed mean values for AUC 0-120 were reported by the authors in their clinical study. The predicted population mean for AUC 0-120 for each model was calculated by the simulator. This applies to both validation approaches (IV bolus and IV infusion).
In the multiple IV bolus administration, the authors of the clinical study used as the validation dataset reported for DOX, a t 1/2 a (distribution half-life) of 4 ± 2 min. The injection was administered over 1 min. However, the first sampling time was 5 min. Based on these facts, we believe that if the biexponential equations the authors suggested based on their observations are used for the calculation of C max values, then they would be underestimated. To this end, and since as noted in the study, the authors measured the concentration of DOX for one patient every 30 s after the administration and found the C max value to be 9.98 mg/L at 90 s, we elected to use this value as the observed C max for the multiple IV bolus administration, as we consider it to be more accurate. The predicted C max values were calculated by the simulator. In the multiple IV infusion administration, the mean observed C max value is directly reported by the authors. The predicted C max values were calculated by the simulator. * This term refers to DOX distribution model. CL met : represents the metabolic clearance for each model calculated either by 10 6 hepatocytes (HEP) or using our custom distribution (DIST) of metabolic clearance on in vitro systems (see Table 6).

Discussion
To accomplish the goals, the development of all DOX models followed the middleout approach as presented above. The data regarding DOX physicochemical properties were collected from the literature or calculated, and two independent clinical studies were selected. The first served as the training dataset (clinical study by Camaggi et al.) and the second as the validation dataset (clinical study by Speth et al.) for the models. In particular, the training dataset was used for the calculation of DOX clearance. DOX distribution was calculated using data from the training dataset as well as utilizing the Simcyp parameter estimation tool (P.E. tool). The validation dataset was used to validate the models' performance. There were three key points where different approaches could be used, thus leading to eight different possible models each of which was validated against the validation dataset. The validation was performed comparing the observed versus the predicted values for AUC and C max of DOX.  Table 8, the best performance based on C max was observed on models 3 and 7 [−9.3% and −9.4% (relative difference between predicted and observed values), respectively] followed by that of models 1 and 5 (+11.2% and −16.3%, respectively). It seems that models based on the minimal PBPK distribution model better predict C max , while models using the full PBPK model seem to overestimate C max by approximately 60%. This is an intriguing result; however, by considering that: (a) in the clinical case study used as the training dataset the authors measured the first concentration value at 15 min, (b) the observed C max value was calculated based on the triexponential equations the authors calculated and (c) the first half-life of DOX is approximately 5 min, then one could assume that the authors probably had underestimated C max .
The best-performing models when using AUC 0-168 as the criterion seem to be models 4, 2, 7 and 5 (+0.07%, +0.16%, −0.24% and +0.56%, respectively), without significant differences for the rest. The results were satisfying, ranging from −1.79% to +1.9%. Overall, as expected, the generated DOX models fit the data of the training dataset.

Model Performance Based on IV Bolus Validation
The first validation of the models was performed using the IV bolus methodology described in the validation dataset clinical study. As seen in Table 10, when using C max as the criterion, the best performing models are 2, 4, 6 and 8 (−30.9%, −30.9%, −31.0% and −31.0%, respectively). In fact, it seems that the most accurate models are those using full PBPK model for distribution. It must be noted, however, that the observed C max value corresponds to one patient whose plasma DOX concentration was measured at 90 s after administration, while the predicted values correspond to approximately 60 s after administration. For that reason, we elected to utilise the AUC predicted/observed ratio for the selection of the most appropriate model.
Utilizing AUC 0-120 as the criterion, models 8, 7, 5 and 6 (+9.59%, +9.64%, +9.90% and +10.33%, respectively) seem to be most accurate. Among them, the best performance was by models using the approach by Wright et al. when predicting renal clearance (models 8 and 7 versus models 5 and 6). Considering all the above, it seems that the best modelling approach would be to use the full PBPK model for distribution (Method F), combined with the approach based on the work by Wright et al. when calculating renal clearance (Method B).

Model Performance Based on IV Infusion Validation
The second validation of the models was performed using the IV infusion methodology described in the validation dataset clinical study. As seen in Table 11, when using C max as the criterion, the best-performing models are 2, 4, 6 and 8 (−59.4%, −59.5%, −59.7% and −59.7%, respectively), thus also indicating that the models using the full PBPK model for distribution tend to better predict the C max ; however, the rest of the models show a similarly good prediction for C max (values ranging from −62.0% to −65.4%).
On the other hand, using AUC 0-120 as the criterion, the best-performing models seem to be 7 and 5 followed by 8 and 1 (24.71%, 25.08%, 25.67% and 25.91% respectively). Thus three (5,7,8) out of four best-performing models use Method B for the calculation of renal clearance, also three (1,5,7) use the mPBPK model for distribution and two of the best performing models calculate hepatic metabolic clearance by each method (C and D).
It is noteworthy, however, that the differences between observed and predicted values are relatively small and well within the two-fold ratio for acceptance set by the industry as a performance standard [30][31][32].

Limitations of the Models Based on Procedure
During the procedure followed for the development of the different models, there were three fundamental limitations that were knowingly ignored at this point but will be thoroughly investigated in the future.
The first limitation was that as can be observed from the training dataset data the renal clearance of the drug exceeds the expected GFR of the patients thus indicating the participation of active procedures in renal excretion. Despite that, we calculated the renal clearance of a 20-30 y.o. male solely based on the expected GFR rate.
The second limitation is about the calculated blood to plasma ratio of DOX. The value used refers to rats. Thus, it is possible that it slightly differs in humans. The future development of our model will consider this factor as well.
The third and last limitation is the distribution of metabolic clearance into its different components. Based on qualitative descriptions, we attempted to quantify the relative percentage of each path. This will also be further investigated as we expand our models.

Conclusions
Based on the data presented in this work, all generated DOX models perform quite well according to the two-fold standard that is widely considered the acceptable measure for the performance of PBPK models.
However, depending on the relevant application where the model might provide insights (e.g., whether the prediction of plasma or specific tissue concentration is in question), one should choose the appropriate model. In the case that DOX nanoformulation pharmacokinetics is under consideration, the selection of a model utilizing full PBPK for distribution is probably most appropriate. The latter also applies if the purpose of the model is the prediction of DOX metabolites or DOX toxicokinetics as more accuracy is needed in the prediction of concentration at the tissue level. To further investigate the above-mentioned problematic PK behaviour of DOX (relating to gender among other pharmacogenomic factors), more clinical data are needed to refine and validate the model for also predicting DOX and its metabolites in various tissues. Overall, we believe that the developed models form a solid basis for further development of even more informative models expanding to new formulations and pharmacogenomic investigations for DOX.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/pharmaceutics14030541/s1. Section S1: Pharmacokinetic data of DOX from the training dataset. Section S2: Plasma Protein Binding. Section S3: Calculating Renal Clearance for each patient. Section S4: Calculating liver blood flow for each patient. Section S5: DOX blood to plasma ratio. Section S6: Hepatic Blood clearance (CLH.B) of each patient. Section S7: Calculating mean fraction excreted in urine for the patients. Section S8: Intrinsic Hepatic Clearance (CLu int,H,b ), Intrinsic Biliary Excretion (CLu bile,b ) and Intrinsic Metabolic Clearance (CLu int,met,b ) for each patient. Section S9: Calculating Liver Weight, HPGL, MPPGL and CPPGL for each patient. Section S10: Biliary excretion per 10 6 of patient hepatic cells. Section S11: Calculating the corrected intrinsic metabolic clearance of each patient based on the three different in vitro systems. Section S12: Explaining Simcyp distribution models. Section S13: Figures of DOX population concentration vs time for the 8 models based on the Training dataset. Section S14: Figures of DOX population concentration vs time for the 8 models based on the Validation dataset multiple IV Bolus administration. Section S15: Figures of DOX population concentration vs time for the 8 models based on the Validation dataset multiple IV Infusion administration. Table S1: Pharmacokinetic data of DOX from the training dataset.  Table S7: Intrinsic Hepatic Clearance (CLu int,H,b ), Intrinsic Biliary Excretion (CLu bile,b ) and Intrinsic Metabolic Clearance (CLu int,met,b ) for each patient. Table S8: Calculated Liver Weight, HPGL, MPPGL and CPPGL for each patient. Table S9: Calculated biliary excretion per 10 6 of patient hepatic cells. Table S10: Calculated corrected intrinsic metabolic clearance of each patient based on the three different in vitro systems. Figures S1-S7