Establishment and Evaluation of a Parametric Population Pharmacokinetic Model Repository for Ganciclovir and Valganciclovir

Background: Ganciclovir and valganciclovir are used for prophylaxis and treatment of cytomegalovirus infection. However, there is great interindividual variability in ganciclovir’s pharmacokinetics (PK), highlighting the importance of individualized dosing. To facilitate model-informed precision dosing (MIPD), this study aimed to establish a parametric model repository of ganciclovir and valganciclovir by summarizing existing population pharmacokinetic information and analyzing the sources of variability. (2) Methods: A total of four databases were searched for published population PK models. We replicated these models, evaluated the impact of covariates on clearance, calculated the probability of target attainment for each model based on a predetermined dosing regimen, and developed an area under the concentration–time curve (AUC) calculator using maximum a posteriori Bayesian estimation. (3) Results: A total of 16 models, one- or two-compartment models, were included. The most significant covariates were body size (weight and body surface area) and renal function. The results show that 5 mg/kg/12 h of ganciclovir could make the AUC0–24h within 40–80 mg·h/L for 50.03% pediatrics but cause AUC0–24h exceeding the exposure thresholds for toxicity (120 mg·h/L) in 51.24% adults. (4) Conclusions: Dosing regimens of ganciclovir and valganciclovir should be adjusted according to body size and renal function. This model repository has a broad range of potential applications in MIPD.


Introduction
Ganciclovir (GCV) is often used for, but not limited to, the treatment or prophylaxis of cytomegalovirus (CMV) disease in immunocompromised patients, such as those with acquired immunodeficiency syndrome and iatrogenic immunosuppression associated with organ transplantation or chemotherapy of neoplastic disease [1]. GCV's antiviral activity is mediated through its triphosphate, which inhibits viral DNA polymerase and slows DNA elongation [2]. Valganciclovir (VGCV) is an L-valine ester prodrug of GCV formulated as an oral solution or tablet to overcome the low bioavailability of oral GCV (~5%) [3]. VGCV has a higher oral bioavailability of approximately 60% [4]. Plasma protein binding of GCV is minimal, accounting for approximately 1 to 2% over the concentration range of 0.5-51 mg/L [5]. GCV clearance correlates well with renal function as over 90% is eliminated unchanged in urine [6].
Currently, there is no consensus on the pharmacokinetic-pharmacodynamic (PK-PD) endpoints for GCV due to a paucity of clinical studies examining efficacy targets in adults or pediatric populations. Both the trough concentration (C trough ) and the area under the concentration-time curve (AUC) have been used for therapeutic drug monitoring (TDM). Recently, a 24 h AUC (AUC 0-24h ) above 50 mg·h/L has been suggested by Märtson et al. [7] for prophylaxis and 80-120 mg·h/L for treatment of CMV infection.
Several studies have shown that GCV exhibited high interindividual variability (IIV) in different populations [8][9][10][11]. As a result, the use of a one-size-fits-all dose for all patients could result in treatment failure. Individualized dosing and TDM play crucial roles in optimizing GCV efficacy and safety. Population pharmacokinetics (PPK) analyses can identify covariates that influence PK parameters and produce estimates of individual PK parameters through Bayesian forecasting to develop individualized therapy for GCV.
To achieve individualized dosing, a comprehensive understanding of the PK parameters and PPK models of GCV is necessary. However, to our knowledge, there are no published PPK model repositories. Therefore, the overarching aim of our work was to establish a parametric PPK model repository of GCV and VGCV. We would like to highlight that the goal of this study is not to prove that these models are comparable with each other, but rather to showcase the heterogeneity. Specifically, we attempted to (1) identify all published PPK models to the best of our capability through a comprehensive literature search, and to reorganize these models using standardized R codes; (2) assess the performance of different PPK models in this repository by comparing the concentration-time profiles and covariates effects; and (3) illustrate the utility of the model repository with two examples (i.e., Monte Carlo simulation for the probability of target attainment (PTA) and AUC calculator based on maximum a posteriori Bayesian estimation (MAP-BE)). By building this model repository and demonstrating the sources of variability, we aim to facilitate subsequent utilization of these models to improve precision dosing of GCV and VCGV.

Search Strategy
In an attempt to identify all parametric PPK models, a literature search was performed on the PubMed, Embase, Scopus and Web of Science databases from inception to 28 May 2023, according to the Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) reporting guideline [12]. Search terms included those related to the drug of interest (ganciclovir, BW-759, valganciclovir, Cytovene, and Valcyte) and terms specific to PPK models included population pharmacokinetic, NONMEM, MONOLIX, and nlmixr. The latter terms were derived from a publication by Li et al. [13]. Two authors (W.Y. and M.G.) conducted the literature search and study selection independently using EndNote (version 20.0.0; Thomson Scientific, Box Hill, Victoria, Australia). A third senior investigator (Q.H.) was consulted to resolve any discrepancies between the two authors. The complete search strategies for each database, and inclusion and exclusion criteria are available in Supplementary Materials.

Information Extraction
We extracted the following information from the included studies: (1) study characteristics including subject demographics, sampling strategies, dosage regimens, and quantitative methods; (2) PPK modeling characteristics, including software and algorithms used, PK parameters and related formulas, between-subject variability (BSV, which was recorded as the coefficient of variation (CV), and %CV= omega 2 * 100%), residual unex-Pharmaceutics 2023, 15, 1801 3 of 20 plained variability (RUV), validation methods; and (3) covariate information, including the list of all covariates that were tested, selection criteria and the subset of covariates that had significant effects on the PK parameters.

Quality Control of the PPK Model Repository
Quality control (QC) procedures were undertaken to screen and rectify any issues related to the establishment of the model repository. We created 4 age-stratified cohorts of virtual patients (neonates, infants, children and adults), these typical virtual patients were designed to reflect the target population of each model as accurately as possible (See Table 1 for details). Although the age covariate had been converted into a categorical variable, this should not impede the validity of the virtual patient cohort as the representative values of each cohort were derived from the median of observed data. The corresponding steadystate concentration-time curves for each virtual patients were obtained by simulations (each simulation was performed on 1000 typical virtual patients) with rxode2 package (version 2.0.12) of R (version 4.2.2). The principle that guided the QC procedures was that, given the models were appropriate to describe the corresponding PPK characteristics, the simulated concentration-time curves of the same typical virtual patients generated by different models should be comparable. This was defined as the 95% confidential interval (CI) of the geometric mean of time to maximum concentration (T max ), maximum concentration (C max ) and other PK parameters of 1000 virtual patients in a study should fall within the range of 50-200% of the geometric mean of PK parameters of all the same virtual patients. When significant differences in the simulated concentration-time curves were detected (i.e., the C max was outside of this range), potential model reproduction errors were first examined and excluded, and then the covariates included in the models were compared to identify possible causes that affected the PK behaviors. The PK parameters used for similarity comparison, including C max , T max , and the half-life time (t 1/2 ) were derived from non-compartment analysis (NCA) analysis of drug concentrations in virtual patients by using IQnca (Version 1.3.0) and Rmisc (Version 1.5.1) package.
All virtual patients received 5 mg/kg GCV by intravenous infusion over 1 h every 12 h (q12h). For VGCV, neonates, infants, and children received 10 mg/kg q12h while adult received 900 mg q12h. For pediatrics, their dosing regimens followed the commonly used dosing regimens of the included studies; for adults, the dosing regimen was based on the third international consensus guidelines on the management of cytomegalovirus in solid-organ transplantation [14]. The R codes of model repository establishment were provided in Supplementary Materials.

Effect of Covariates on Clearance Variation
Clearance (CL) is a crucial parameter for AUC, and AUC plays a central role in the individualized dosing of GCV. Thus, the comparison of the effects of different covariates on CL was necessary. To explore whether the covariate effect on CL was clinically meaningful and to understand between-study differences in covariate impacts on CL, we used a forest plot to comprehensively compare covariate effects across studies. Weight, estimated glomerular filtration rate (eGFR) and creatinine clearance (CLcr) were scaled to the same  a uniform covariate value was set as the reference (refer to Table S1 for details). For other continuous covariates (serum creatinine, SCR and body surface area, BSA) that were only identified in one study, the minimum and maximum values of that covariate were used to calculate CL and the reference CL value was calculated using the median covariate value of each study.
For binary covariates, such as the presence of critical illness (defined as "1" for critically ill patients and "0" for others), the common condition would be treated as the reference (COV i = 0). The uncommon condition would be treated as the test (COV i = 1). CL i = CL common + CL diff * COV i . The range of CL i would be [CL common , CL commpn + CL diff ] (if CL diff > 0), or [CL common + CL diff , CL common ] (if CL diff < 0).
Then, the effect range of identified covariate on CL was calculated by the following formula: The minimum CL or the maximum CL Reference CL × 100% We considered covariates effects outside of the 80-125% boundary as clinically significant, based on the standard used in bioequivalence studies [13,15]. Detailed R codes can be found in the Supplementary Materials.

Application of the PPK Model Repository
In order to demonstrate the practical application of the model repository in MIPD, we will provide two illustrative examples.

Monte Carlo Simulation for the Probability of Target Attainment
We calculated the PTA of simulated concentration-time curves to evaluate the commonly used dosing regimens of GCV and VGCV. The trapezoidal method was used to calculate the steady-state AUC 0-24h for 1000 virtual patients. The prophylaxis target used was 40-80 mg·h/L [16], while AUC 0-24h of 80-120 mg·h/L was considered treatment target [7]. AUC exceeding 120 mg·h/L posed toxic risk.

AUC Calculator Based on Maximum a Posteriori Method
When using GCV and VGCV in a real-world setting, clinicians often need to calculate the AUC to determine whether the patient has achieved the appropriate level of exposure. To highlight the model repository's convenience in MIPD, we developed an AUC 0-24h calculator based on the maximum a posteriori-Bayesian Estimation (MAP-BE) using R shiny (version 1.7.4). Four main components made up the core function of the calculator: (1) PPK model parameter information; (2) model defined by rxode2 package; (3) the objective function; (4) function used for MAP estimation [17].

Identification of the Included Studies
598 studies were found in the initial search across various databases. A total of 16 studies were included in the final analysis, as shown in Figure 1. No additional records were identified from other sources.
BSV was described by an exponential model in all the included studies. RUV was described by a proportional model in ten studies [8,10,11,18,19,22,24,25,28,29], exponential model in three studies [20,23,26], additive model in one study [9], and combined proportional and additive model in two study. Two studies reported interoccasion variability (IOV). Facchin et al. [28] found an IOV in CL of 14.4%, 77.2% in peripheral volume of distribution (Vp), and 111.4% in absorption rate constant (ka). Perrottet et al. [22] reported an IOV in CL of 12%. Two studies published before 2000 did not report model evaluation results [18,19], but after careful inspection, the performance of these two models was comparable to others. Other studies were evaluated by internal validation; three of them also underwent external validation [11,26,27]. GOF, VPC, normalized prediction distribution error (NPDE) and bootstraps were often used as internal validation methods. Modeling strateg are summarized in Ta two that used Monoli interaction (FOCE-I) w GCV PPK as a two-co cluded with a one-com the pediatric populati gested that the one-co was sparse. In additio likely to use the one-c described as a first-or netics parameters rep (Lalagkas et al. [27] a the conversion equat dosage for inclusion difference in molecul reported bioavailabili ꝏ BSV was describ described by a propo tial model in three stu portional and additive   [18,19], but aft parable to others. Oth underwent external v error (NPDE) and boo

QC result
Similarity comparison was used to ensure accuracy of model repository construction. After excluding construction errors, the 95%CI of the PK parameters' geometric mean was mainly distributed within 70-150% of the geometric means of same virtual patients.
PK behaviors of GCV in the same pediatric typical virtual patients from different PPK models were comparable (Supplementary Material Figure S1). However, t 1/2 of Horvatits et al. [24] was larger than others because its original data came from critically ill patients receiving continuous venovenous hemodiafiltration (CVVHDF). It should be noted that the final PPK model had controlled for the extracorporeal therapy given to the patients. After oral administration of VGCV, the PK behaviors of GCV in the same typical virtual patients from most PPK models were also comparable (Supplementary Material Figure S2) except that virtual infants of Vezina et al. [25] showed lower t 1/2 values, and this could be attributed to a limited 3 out of 95 subjects between 0 and 24 months, thus this model could not describe the process of GCV in typical virtual infants very well.
After QC and subsequent verification, the model repository quality was deemed satisfactory.

Comparison of GCV and VGCV PK Profiles
Simulated GCV concentration-time profiles were displayed in Figure 2. The PK profiles across pediatrics were comparable because the weight-based dosing regimen had largely resolved PK differences between pediatrics, indicating that body weight could significantly influence GCV's PK. Adults showed higher serum concentrations than pediatrics at a similar dose of 5 mg/kg. Furthermore, simulated concentrations based on the models established by Krens et al. [10] and Horvatits et al. [24] showed high variability when the dosing regimen of 5 mg/kg q12h was used. Figure 3 shows the concentration-time profiles of VGCV. In pediatric groups, the simulated profiles were comparable between published models. Adults had higher concentrations than pediatrics.
Simulated GCV concentration-time profiles were displayed in Figure 2. The PK profiles across pediatrics were comparable because the weight-based dosing regimen had largely resolved PK differences between pediatrics, indicating that body weight could significantly influence GCV's PK. Adults showed higher serum concentrations than pediatrics at a similar dose of 5 mg/kg. Furthermore, simulated concentrations based on the models established by Krens et al. [10] and Horvatits et al. [24] showed high variability when the dosing regimen of 5 mg/kg q12h was used.    All patients were assumed to be males. For neonates, infants and children, VGCV monotherapy was given at a dose of 10 mg/kg q12h, while 900 mg q12h was for adults.

Covariate Screening and Covariate Effect
All tested covariates that had an effect on CL, distribution volume of central compartment (Vc), intercompartment clearance (Q) and distribution volume of the peripheral compartment (Vp) are summarized in Table 4. The stepwise method typically employed for covariate screening included forward inclusion and backward elimination. No covariates were investigated in Horvatits et al. [24] due to the limited number of included sub- All patients were assumed to be males. For neonates, infants and children, VGCV monotherapy was given at a dose of 10 mg/kg q12h, while 900 mg q12h was for adults.

Covariate Screening and Covariate Effect
All tested covariates that had an effect on CL, distribution volume of central compartment (Vc), intercompartment clearance (Q) and distribution volume of the peripheral compartment (Vp) are summarized in Table 4. The stepwise method typically employed for covariate screening included forward inclusion and backward elimination. No covariates were investigated in Horvatits et al. [24] due to the limited number of included subjects (n = 9). The most influential covariates were body weight and renal function indicators, such as eGFR and CLcr. The impact of each covariate on CL is shown in Figure 4. Body size, including weight and BSA, were evaluated and included as significant covariates in 9 (56.3%) studies. Compared to the reference value, seven out of nine demonstrated significant impact of body size on CL with greater than 20% change under the normal range of body size [9,11,18,20,25,27,29]. Furthermore, the effect of renal function, such as CLcr, eGFR or SCR, was reported in 13 (81.3%) studies. Given the normal range of renal function, all of them demonstrated greater than 20% change in CL when compared to the reference value. Renal function affects renal clearance of GCV, which is the primary eliminate pathway of GCV. Only two studies investigated the influence of sex on CL (limited impact with less than a 20% difference) [22,28]. Yuen et al. [18] and Nguyen et al. [11] also investigated the effect of transplants, CMV-shedding, and critical illness on CL, where only transplant and CMV-shedding showed a significant influence. Body size and sex were reported to be the significant covariates on Vc. The only covariate identified on Q and Vp was weight.     Figure 5 depicts PTA for commonly used dosing regimens based on published PPK models, with each model characterizing a specific population. Results indicated that 5 mg/kg q12h GCV dosing in adults carried a risk of overexposure (AUC 0-24h above 120 mg·h/L) with three out of six models suggesting that half of the subjects were overexposed [10,22,24]. According to the PTA of all PPK models involving adults (Table S2), 51.24% adults' AUC 0-24h could exceed 120 mg·h/L with 5 mg/kg/12 h GCV dose regimen. In contrast, this dosing regimen showed good prophylaxis in most pediatric studies, 46.44% pediatrics could reach the prophylaxis target.  Figure 5 depicts PTA for commonly used dosing regimens based on published PPK models, with each model characterizing a specific population. Results indicated that 5 mg/kg q12h GCV dosing in adults carried a risk of overexposure (AUC0-24h above 120 mg·h/L) with three out of six models suggesting that half of the subjects were overexposed [10,22,24]. According to the PTA of all PPK models involving adults (Table S2), 51.24% adults' AUC0-24h could exceed 120 mg·h/L with 5 mg/kg/12 h GCV dose regimen. In contrast, this dosing regimen showed good prophylaxis in most pediatric studies, 46.44% pediatrics could reach the prophylaxis target.
Further, 10 mg/kg/12 h VGCV was also effective in achieving prophylaxis in pediatrics, while 900 mg/12 h VGCV could lead to overexposure in adults.

AUC calculator based on MAP-BE
As an example, we used the model established by Franck et al. [9] to develop a MAP-BE-based calculator for calculating the AUC 0-24h post-administration of GCV and VGCV. It requires the patient's dosing and sampling information, weight, and creatinine clearance (CrCL, mL/min/1.73 m 2 ) to calculate AUC 0-24h after the first dose. Its results were consistent with NONMEM (see the Supplementary Materials). This AUC calculator demo is available online. (https://ganciclovir.shinyapps.io/Example_ GCVAUCCalculator/). Researchers can create similar calculators by modifying the model code (Supplementary Materials).

Discussion
To our best knowledge, this study was the first to build and share a parametric PPK model repository of GCV and VGCV. The repository is characterized by simulations of concentration-time profiles and covariate effect evaluation, where we demonstrated the potential of the constituent models in estimating the AUC and PTA of GCV and VGCV. The work provides evidence for individualized dosing based on the patient's weight and renal function, as the commonly used dosing regimens tend to miss treatment targets in pediatrics.

Simulated Concentration-Time Profiles of Ganciclovir
The simulated concentration-time profiles (Figure 2) of most studies demonstrated a comparable curve in pediatrics when a dose of 5 mg/kg GCV was administered, suggesting that weight-based dosing is critical. Three studies' (Horvatits et al. [24], Krens et al. [10], and Li et al. [29]) simulated profiles that were different from the others. We found that all three studies included critically ill patients in whom PK can be complicated by hemodynamic instability with varied volumes of distribution and fluctuating renal function. For the study of Horvatits et al., the model was also limited by the small number of patients (9 patients), with no covariates successfully included, and the IIV not clearly explained.
There were three other possible reasons why Li et al.'s PK profiles were different from others: (1) only 11.5% of the patients were children over six years old, so the 10-year-old virtual children did not match with the studied population; (2) the very sparse sampling strategy (138 samples from 104 subjects) might also limit the model performance; (3) the Gao formula [30] used in this study is more applicable to calculate eGFR of children with moderate renal failure.

Simulated Concentration-Time Profiles of Valganciclovir
The simulated concentration-time profiles of VGCV ( Figure 3) were similar across pediatric studies except for Zhao et al. [23]. The patient population in Zhao et al. was unique as they were children who received a kidney transplant and mycophenolate mofetil as the immunosuppressant. It is known that mycophenolate mofetil could reduce the renal clearance of GCV, thus resulting in higher plasma concentrations in the simulation [31]. In Perrottet et al. [22], patients were stratified into three subgroups according to the type of organ transplant received. Patients who received a heart transplant had significantly higher exposure to GCV than the other two subgroups (liver and lung). The use of tacrolimus (heart subgroup) could have reduced renal blood flow to a greater extent than cyclosporine (liver and lung subgroups), resulting in the reduced renal elimination of GCV and the subsequent larger exposure to the drug [32]. Under the commonly used dosing regimen, the exposure level of VGCV in adults was similar to that of GCV at different ages. However, the exposure in pediatrics was significantly lower than GCV, indicating that dose optimization of VGCV for pediatrics is urgently needed.

Covariates Effects on Estimated PK Parameters
Eight of the 16 included studies identified body weight as a covariate on CL. Our results suggested that as age increased, the influence of weight on CL decreased. In neonates, the effect of weight on CL was marked with the clearance of those at the highest weight (5 kg) being at least 3.34 fold the CL in neonates weighing 1 kg. In contrast, the difference in clearance between the heaviest adult patient weighing 100 kg and the lightest patient weighing 40 kg was at most only 2.43 fold. This phenomenon may occur because pediatric patients are in a state of continuous growth and development, so the fluctuations are larger than in adults (the percentage of weight gain decreases with age). This is consistent with the notion that weight is often considered a significant covariate in pediatric population studies but is rarely included in adult studies.
Thirteen of the sixteen included studies identified renal function as a covariate on CL. GCV is mainly excreted in the urine through glomerular filtration and active tubular secretion. As such, the patient's renal function would have a greater impact on the PK behavior of GCV. The forest plot suggested that regardless of the renal function indicator (i.e., eGFR, CLcr, SCR), the covariate effect on CL was clinically significant (outside the range of 80%-125%). This indicates that testing renal function's influence is important when constructing the PPK model of GCV.
Two studies identified gender as a significant covariate, but the subsequent effect evaluation suggested that gender was not a clinically meaningful covariate on CL. The gender effect could manifest in other covariates, such as a discrepancy in body size.

Envisioning the Application of Model Repository
This study constructed a model repository of GCV and VGCV population pharmacokinetic models. Our research has the potential to significantly decrease the amount of time required for conducting literature searches. Leveraging our repository, other researchers can quickly perform external evaluations on their data to identify a suitable model or to explore the factors that influence the predictive ability of the model in different clinical settings. In addition, researchers can perform model averaging to reduce the impact of uncertainty in a single model and to identify the most appropriate predictions for individual patients, thus simplifying the process of precision dosing and reducing the burden of model validation.
New dosing regimens can be evaluated comprehensively using the models in the repository and the provided codes. Our work could significantly improve the efficiency of dosing regimens evaluation and the AUC calculator based on MAP-BE could be useful for determining AUC. The source code related to the PPK model could be modified to develop calculators that are based on new models. In addition, an AUC calculator can be built for each model in the repository, and medical staff can make decisions considering the AUC results calculated by multiple models.
The model repository can also be docked with the "PopED" package, and the published model in the repository can be selected as the prior information to optimize the experimental design. The package also supports the use of "rxode2" to build a preset model. The "nlmixr" package, as a free and open-source package, supports nonlinear mixed effects analysis in R, and can be seamlessly connected with the relevant code of the model repository, improving the efficiency of modeling.

Limitations
Our study has some limitations. Firstly, only the parametric PPK models were included in this study, and the non-parametric PPK models were excluded because the parameters of the non-parametric models were hard to bridge to parametric models. Secondly, in our simulation, we used several typical virtual patients in order to facilitate comparisons across studies. Hence, the simulation results may not be representative enough for the distribution of covariates in clinical practice. In other words, the typical patients may not represent the populations studied. Thirdly, the construction of the model repository was done manually in this study, and the next step could be to automatically generate a model repository using artificial intelligence technology [33].

Conclusions
The model repository of parametric population pharmacokinetic models for GCV and VGCV is useful for promoting MIPD. Optimization of the GCV and VGCV dosing regimen should consider the patient's renal function. Our AUC calculator may provide a useful tool for clinicians to perform TDM during their routine work.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/pharmaceutics15071801/s1, Table S1: The uniform covariate range and reference covariate value; Table S2: Total probability of target attainment of pediatrics and adults; Figure S1: Similarity comparison (%) of simulated pharmacokinetic profiles after intravenous infusion of ganciclovir; Figure   Data Availability Statement: The data generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.