Tacrolimus (Tac) reduces the risk of rejection after solid organ transplantation. However, its toxicities are well known (Kershner and Fitzsimmons, 1996). Consequently, many transplant professionals and pharmacologists have to manage its narrow therapeutic window. Given its highly variable pharmacologic (pharmacokinetic [PK] and pharmacodynamic [PD]) profile (Staatz and Tett, 2004), therapeutic drug monitoring (TDM) is used to individualize Tac dosages and reduce the risks of toxicity and rejection. However, traditional TDM remains a reactive strategy that requires a PK steady state, i.e., approximately 3 days after therapy initiation or dosage change. The delay caused by repetitive dose changes is prohibitive in the early achievement of safe and effective Tac levels (Staatz et al., 2001; Borobia et al., 2009; Richards et al., 2014). Identification of invariable PK biomarkers can help to proactively adjust the dose. However, the identification of useful and relevant biomarkers is only the first step toward therapy individualization. Once the marker is identified, its effect on drug PK variability must be quantified. A population-based PK (popPK) approach can help to model quantitatively the effect of a patient covariate on the drug PK profile in order to simulate the most probable response for a given patient allowing the design of a personalized drug dosage.

Tacrolimus is metabolized in the intestine, in the liver and, to a limited extent, in the kidney by the Cytochromes P450 (CYP) 3A4 and 3A5 enzymes (Staatz and Tett, 2004; Kamdem et al., 2005). It is now universally recognized that a single nucleotide polymorphism (SNP) in the CYP3A5 gene is associated with approximately a 40 to 50% decrease in Tac clearance (Anglicheau et al., 2003; Hesselink et al., 2003; Haufroid et al., 2004, 2006; Macphee et al., 2005; Elens et al., 2007). Inclusion of CYP3A5^∗3/^∗3 loss-of-function (LOF) allelic status for Tac initial dosage calculation achieves therapeutic levels more quickly (Thervet et al., 2010). However, despite the PK improvement it generates, the clinical benefit in terms of outcome of such a pro-active dosage strategy has not been proven yet but some limitations in study designs have been highlighted (van Gelder and Hesselink, 2010). However, even the PK benefit of a proactive dosage based on CYP3A5 genotype solely is controversial (Shuker et al., 2016).

Recently, the CYP3A4^∗22 decrease-of-function (DOF) allele has been suggested as a good candidate to further refine the Tac starting dose (after adjusting for CYP3A5 genotype) (Elens et al., 2011, 2013a,b,c,f; Gijsen et al., 2013; Guy-Viterbo et al., 2014; Hesselink et al., 2014; Kuypers et al., 2014; de Jonge et al., 2015a; Tang et al., 2016). However, to our knowledge, only two PopPK studies have examined the combined effects of the CYP3A4^∗22 and CYP3A5^∗3 SNPs, showing that CYP3A4 DOF can exacerbate the CYP3A5 LOF (Moes et al., 2016; Andreu et al., 2017). The definition of a rationale categorization of the patient into poor (PM), intermediate (IM), and extensive metabolizer (EM) according to these two SNPs has been successfully proposed previously and takes the advantage of being clearly understandable for clinicians and medical staff. The classification is based on the fact that the CYP3A4^∗22 allele is associated with a decrease of CYP3A4 function while CYP3A5^∗3 is linked to a loss of CYP3A5 expression and that both metabolic defects have synergistic effects. Rationally, the PM cluster contains CYP3A5^∗3 homozygotes carrying the CYP3A4^∗22 variant; the IM group contains CYP3A5^∗3 homozygotes but not carrying the CYP3A4^∗22 allele; and EM includes CYP3A5 expressers also not carrying the CYP3A4^∗22 allele.

Apart from these functional SNPs, there are a plethora of satellite genes that could affect the function of CYP3A isoenzymes (Werk and Cascorbi, 2014). We describe two of the most promising SNPs among these genes.

Genetic variation in the Peroxisome proliferator-activated receptor a (PPARA) gene, a nuclear receptor, was discovered as a novel genetic determinant influencing CYP3A4 activity (Klein et al., 2012). The minor allele of the PPARA rs4253728G > A polymorphism has been associated with significantly decreased CYP3A4 expression and activity (Klein et al., 2012; Elens et al., 2013e). This polymorphism might therefore influence the pharmacokinetics of drugs that are primarily metabolized by the CYP3A4 enzyme, such as Tac.

P450 oxidoreductase (POR) is a membrane-bound protein, which is responsible for the transfer of electrons from NADPH to microsomal type II cytochrome P450 enzymes. Liver-specific POR-knockout mice are phenotypically normal but accumulate lipids in the liver and show considerably decreased hepatic drug metabolism (Henderson et al., 2003). Numerous POR missense mutations in humans have been discovered and linked to anarchic steroidogenesis, ambiguous genitalia, and Antley–Bixler syndrome (Huang et al., 2008). In the general population, the 1508C > T SNP (rs1057868; POR^∗28) is the most common variant with a reported minor allelic frequency (MAF) of 30% in the white population. POR^∗28 encodes the amino acid variant A503V, which has been associated with differential CYP450 activity (de Keyser et al., 2013; Elens et al., 2013d, 2014). For instance, CYP3A5 expressers carrying one or two POR^∗28 alleles have shown a 45 % lower midazolam metabolite conversion and higher Tac dose compared with CYP3A5 expressers without POR^∗28 (Elens et al., 2013d).

The first objective of our study is to use a population-based PK approach to simultaneously evaluate the relevance of genotypic and non-genotypic covariates formerly identified as influencing Tac PK. The second objective is to translate our findings into rationale initial dosage recommendations for clinicians that maximize the probability of achieving desired Tac concentrations after the initial dose.

Baseline characteristics of the patients and genotype frequencies are reported in Table 1. The genotype distributions were in accordance with the Hardy-Weinberg principle and with the frequencies reported¹. In total, considering the CYP3A clustering strategy described earlier (Elens et al., 2011), there were 5 patients classified as poor metabolizers (PM = CYP3A5^∗3 homozygotes carrying the CYP3A4^∗22 variant), 36 as Intermediate metabolizers (IM = CYP3A5^∗3 homozygotes not carrying the CYP3A4^∗22) and 18 as extensive metabolizers (EM = CYP3A5 expressers not carrying the CYP3A4^∗22 allele).

In classical candidate-gene association studies, the effect of CYP3A4^∗22 on Tac PK is well accepted (Elens et al., 2011, 2013a,b,c,f; Gijsen et al., 2013; Guy-Viterbo et al., 2014; Hesselink et al., 2014; de Jonge et al., 2015a; Tang et al., 2016; Andreu et al., 2017). Contrasting with these observations, a number of previous studies have failed to highlight the benefit of introducing CYP3A4^∗22 in modeling Tac inter-individual variability through popPK-approaches, with a few exceptions (Shi et al., 2011; Zuo et al., 2013; Zhang et al., 2015; Moes et al., 2016; Andreu et al., 2017). This is not surprising as the majority of studies were performed in Asian populations where CYP3A4^∗22 is absent, as it is for individuals of African origins. We clearly showed here that PK prediction can be improved by inclusion of patient CYP3A4^∗22 allelic status, particularly via a previously described CYP3A cluster classification that takes into account both CYP3A4^∗22 and CYP3A5^∗3 alleles. Even if it is well accepted that CYP3A5^∗3 is the primary factor explaining Tac PK metabolic defect, we showed here that the CYP3A4^∗22 PK influence is additive. However, even if the amplitude of the CYP3A4^∗22 effect might be comparable to that of CYP3A5^∗3, the CYP3A4^∗22 influence is not as statistically significant probably because of the wide PK variability observed among the CYP3A4^∗22 carriers. As a consequence, it may also explain why it is not always identified as a significant covariate when considering the few studies reported to date (Shi et al., 2011; Zuo et al., 2013; Zhang et al., 2015; Moes et al., 2016). Another possible explanation is the fact that CYP3A5^∗3 completely blunts the CYP3A5 activity whereas for CYP3A4^∗22, some isoenzyme activity remains. An alternative hypothesis is that this lack of statistical reproducibility is due to the lower allelic frequency of CYP3A4^∗22 compared to CYP3A5^∗3 on the one hand (Bigdeli et al., 2014), and to the fact that CYP3A4 activity is more variable than that of CYP3A5 on the other. Furthermore, the parametric (or semi-parametric) modeling strategy used in previous studies is probably less efficient in detecting inter-individual variability. Indeed, our non-parametric approach benefits from using multiple support points for iterative processing of the data and, as such, each patient is considered as having its own PK parameters distribution and does not rely on the supposition that PK parameters are normally distributed in the general population. This allows better individual prediction and increases the ability to detect differences between individuals. Finally, many confounders can potentially impact on the CYP3A4^∗22/CYP3A5^∗3 effect. In the present study, patients were still hospitalized and environmental influencing factors were potentially better controlled than in ambulatory studies. One can also consider the fact that patients were in the very early period after transplantation where steroids that are known to induce CYP3A activity are still at a quite high dosage. Consequently, steroid induction can have a different impact on CYP3A activity depending on the genetic profile and can boost the difference between the different CYP3A clusters making the effect of CYP3A4^∗22 even more significant (de Jonge et al., 2015b). Consequently, in previous studies, the effect of CYP3A4^∗22 might still be clinically true but just hidden because of the study design and/or modeling method.

Importantly, our study is in total agreement with the conclusions of the study of Andreu et al. (2017) even if the research strategies were dissimilar. Indeed, some important differences in both the design, as well as in the population modeling method render our study different and more robust than the Spanish study. The first difference resides in the study design. In the discovery cohort used to build their model, Andreu et al. (2017), included only 7 patients that were intensively sampled and, as such, providing a complete PK course early after transplantation (day 7). The rest of the samples (98% of the patients) were through levels collected at five different time points with only one collected in the very early post-transplant phase. In our study, a complete PK profile was available for all the patients in the early post-transplant phase. The second main difference resides in the population modeling method. Indeed, in the study of Andreu et al. (2017), Tac concentration-time data were analyzed using a parametric population PK approach with NON-MEM, which assumes that the estimated parameters are normally distributed in the population. However, this assumption might not be true especially if under-represented polymorphic alleles and minority clusters are present in the sample. In non-parametric statistics, no assumptions are made about the underlying distribution of the PK parameters and each patient can serve as a support point for the model-building iterative process and the estimation of PK parameters. As such, instead of obtaining only single-point parameter estimates for the population, one gets multiple estimates, up to one for each subject studied. Consequently, the model comes the closest to the collection of each subject’s exactly known parameter values. Other strengths of the non-parametric approaches include mathematical consistency, good statistical efficiency, and good asymptotic convergence (Jelliffe et al., 2000).

PPARα SNPs have been associated with differences in exposure and/or metabolite formation of drugs metabolized through CYP3A. More particularly, PPARα rs4253728G > A SNP has been associated with the risk of developing Post-Transplantation Diabetes Mellitus in patients treated with Tac (Elens et al., 2013e). However, our previous investigation failed to explain this increased risk through a PK difference. Here, data suggest that PPARα might have an influence on Tac PK, but in a recessive manner. However, in our cohort, only 4 patients were homozygous for the variant allele, among whom 2 were CYP3A PM. This obviously renders the statistical power very low and might explain the fact that it was not retained in the final model. Moreover, the effect of PPARα is thought to be exerted through an indirect effect on CYP3A4 activity. As a consequence, its effect in CYP3A4^∗22 carriers might be lowered and potentially confounded. This information might also partly clarify the fact that PPARα SNP was significant only when CYP3A cluster was not included in the model.

By simulations of multiple dosing scenarios across the different CYP3A clusters, we provide here clear dosage recommendations with well-defined deliverables. With the table presented in this paper, the clinician can use our predictions directly for a given patient. Our model was proved to be efficient to predict the Tac though blood concentration obtained after the first dose administered directly after transplantation. Given our PTA prediction table, the results suggest a starting dose around 0.1 mg/kg bodyweight b.i.d. for PM, 0.13 mg/kg bodyweight b.i.d. for IM and 0.2 mg/kg bodyweight b.i.d. for EM. However, by comparing PTA analysis with observed data, even if we can see that predictions were quite accurate for IM and EM, they were less precise for the PM cluster where our model seems to slightly underestimate the defect caused by CYP3A4^∗22. Consequently, in line with what has been proposed earlier (Haufroid et al., 2006; Thervet et al., 2010; Birdwell et al., 2015; De Meyer et al., 2016) and because of the recent shifts toward lower Tac target ranges (Ekberg et al., 2007), we suggest revising our above advices for PM with a dose of 0.07 mg/kg bodyweight b.i.d. These new guidelines are reasonable and in accordance to the original suggestions of Haufroid et al. (2006), whose guidelines have been successfully tested in a randomized clinical trial (Thervet et al., 2010) and further translated in clear recommendations by the CPIC (Birdwell et al., 2015). With the present analysis, we add a slight nuance to their proposal by considering the DOF caused by the CYP3A4^∗22 allele. Besides, our innovative classification implies different PM/IM/EM proportions in each group explaining also the subtle modifications we propose here. However, whereas some studies have identified a relationship between Tac exposure and the risk of acute rejection, this has not been a universal finding (Bouamar et al., 2013). This observation clearly questions the clinical relevance of dosage guidelines based on the probability of trough concentration achievement. It has been speculated that currently applied targets saturate the Tac response and that the concentration-effect relationship is reaching its maximum at lower concentrations (Bouamar et al., 2013; Storset et al., 2015). Nonetheless, with the recent trend toward lower Tac target ranges (Ekberg et al., 2007), the need for prediction tools to avoid underexposure will probably increase. Moreover, Tac is known for its exposure-dependent diabetogenic as well as nephrotoxic effects, which reinforces the relevance of such a tool enabling to avoid too high drug exposure.

Our study has, however, some limitations such as the potential confounding effect of co-medications interfering with ABCB1 function. However, as we did not find any influence of ABCB1 SNPs on the Tacrolimus PK in univariate analysis, it is most likely that these ABCB1 inhibitors will not substantially affect our model, especially given the low dosage of these potentially interacting co-medications. Furthermore, considering these factors would have increased the number of covariates to test and this comes against the general principle of parsimony. Indeed, multiple statistical testing would have amplified the chance of spurious associations leading to over-parametrization of the model. Besides, we did test the potential impact of these co-medications in univariate analysis and no significant associations were found. One second surprising finding is the fact that hematocrit was not retained as a significant covariate in the final model whereas most of previous Tac popPK studies reported a significant effect (Benkali et al., 2009; Woillard et al., 2011; de Jonge et al., 2012; Asberg et al., 2013; Han et al., 2013; Storset et al., 2014, 2015; Andreu et al., 2017). This lack of association can potentially arise from the design of our study. Indeed, given that our patients were still hospitalized and closely monitored, the variability in hematocrit values was not substantial (CV = 15%) and even reduced in PM (CV = 13.7%) and IM (13.6%) after CYP3A genotype stratification, providing an explanation on why the effect of hematocrit is no longer observed in a multivariate context when CYP3A genotype is introduced the model. Also, concerning the lack of influence of the age of the patient, as reported in Table 1, our population was 51.9 years old on average and it has been described that the age-related PK changes were essentially observed between ages 40 and 50 but that bioavailability was constant at lower and higher relative values in younger and older patients, respectively (Storset et al., 2014). This might explain why we did not find any association between age and Tac PK. Similarly, it has been observed that gender differences in the PK of CYP3A substrates seem to be more pronounced at younger ages compared with in the elderly (Cotreau et al., 2005), providing an explanation why gender was not associated with Tac PK in our cohort.

Finally, it is obvious that our recommendations should be validated through a randomized clinical trial and are open to future amendments with the potential discovery of new biomarkers. For instance, different studies highlighted the importance of P450 oxidoreductase SNPs to explain differences in CYP3A-driven metabolism and in particular the POR^∗28 allele. In the present study, we failed to replicate this observation but this could be due to an insufficient statistical power and/or imperfect study design. In particular, the analysis of POR^∗28 is complex as it not only depends on the CYP450 activity alteration but also on the inter-protein cooperation which relies on the substrate size and the CYP450 isoform implied.

We developed here a practical tool to predict Tac exposure after renal transplantation taking into consideration the two patient’s CYP3A4 and CYP3A5 genotypes and their linked predicted metabolic phenotypes. We also showed that our model is accurate in predicting the exposure subsequent to the very first Tac dose, indicating that it can be used even when steady state is not yet reached and thus, produces additional information to TDM that can only be initiated profitably when PK steady state is guaranteed, unless a Bayesian approach is used (Asberg et al., 2013). In conclusion, based on our simulations, we predict a different starting dose for each CYP3A genotype profile. Therefore, we recommend new starting doses of 0.07 mg/kg bid for PM, 0.13 mg/kg bid for IM and 0.2 mg/kg bid for EM. Subsequently, after therapy initiation, this tool would probably benefit the clinician if used in a Bayesian adaptive control system (Storset et al., 2015).

MM, AC, and VH designed the research study; JBW, MN, VH, and LE performed the experiments; JBW, MM, MN, AC, RHvS, TvG, NL, DAH, PM, VH, and LE analyzed the results; JBW, MM, MN, AC, RHvS, TvG, NL, DAH, PM, VH, and LE wrote the manuscript; all authors read and approved the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.