Semiphysiologically Based Pharmacokinetic Model of Leflunomide Disposition in Rheumatoid Arthritis Patients.

A semiphysiologically based pharmacokinetic (semi-PBPK) population model was used to evaluate the influence of enterohepatic recycling and protein binding, as well as the effect of genetic variability in CYP1A2, CYP2C19, and ABCG2, on the large interindividual variability of teriflunomide (active metabolite) concentrations following leflunomide administration in rheumatoid arthritis (RA) patients. The model was developed with total and free teriflunomide concentrations determined in RA patients taking leflunomide, as well as mean teriflunomide concentrations following the administration of leflunomide or teriflunomide extracted from the literature. Once developed, the 15-compartment model was able to predict total and free teriflunomide concentrations and was used to screen demographic and genotypic covariates, of which only fat-free mass and liver function (ALT) improved prediction. This approach effectively evaluated the effects of multiple covariates on both total and free teriflunomide concentrations, which have only been explored previously through simplistic one-compartment models for total teriflunomide.

Study Highlights

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC? ☑ Teriflunomide, the active metabolite of leflunomide, has highly variable pharmacokinetics due in part to multiple pharmacogenetic markers. Teriflunomide is highly protein-bound and free concentrations have been hypothesized to be linked to response. • WHAT QUESTION DID THIS STUDY ADDRESS? ☑ Total and free teriflunomide concentrations following leflunomide administration were characterized in this semi-PBPK model, which can analyze the multiple pharmacogenetic variants that have been linked to altered concentrations. • WHAT THIS STUDY ADDS TO OUR KNOWLEDGE ☑ For the first time, the physiological processes that control the conversion of leflunomide to teriflunomide, as well as teriflunomide enterohepatic recycling and binding effects have been included within a semi-PBPK population model. • HOW THIS MIGHT CHANGE CLINICAL PHARMACOLOGY AND THERAPEUTICS ☑ This model can describe teriflunomide disposition following either leflunomide or teriflunomide administration and thus appears to be an important step in quantifying the differences between the compounds. Given this is the first population model to simultaneously describe free and total teriflunomide concentrations, it may be used within pharmacokinetic-pharmacodynamic models to assess the concentration–effect relationship.

Population pharmacokinetic models have traditionally used empirical mammillary compartmental models of drug disposition,1,2 but these mammillary models are not ideal in their description of hepatic first-pass effects and plasma protein binding.3,4 These potentially important issues can be addressed through physiologically based pharmacokinetic models (PBPK), whereby physiology is mimicked by compartments (representing organs) linked by blood perfusion.2,4 However, complete-PBPK models rely heavily on knowledge of tissue concentrations, and as such application to clinical datasets has been minimal, raising doubts as to their necessity in most cases.2 Semi-PBPK may be more useful, as the complexity is guided by the level of available information and the purpose of the model.2

Leflunomide is commonly used in the treatment of rheumatoid arthritis (RA) and psoriatic arthritis, with new therapeutic uses under investigation.5–7 Activity is mediated via the active metabolite teriflunomide, itself now registered for multiple sclerosis (MS).8 Oral leflunomide is almost completely metabolized during its first pass through the liver and is undetectable in plasma. Approximately 70% of administered leflunomide is converted to active teriflunomide, predominately mediated by the CYP450 enzymes 1A2, 2C19, and 3A4,9,10 and carriage of the loss-of-function CYP2C19*2 allele has been associated with lower teriflunomide plasma concentrations.11 Additionally, teriflunomide production may be mediated by rapid degradation of leflunomide within the blood.10 Teriflunomide's volume of distribution is approximately 11L, consistent with a high degree of plasma protein binding (>99%), and the clearance is 31 mL/h, resulting in an approximate 2-week half-life (mean ∼15.7 days).8,9 Enterohepatic recycling, mediated in part by the ABCG2 transporter, contributes significantly to the long half-life, and genetic variability in ABCG2 (specifically the 421C>A single nucleotide polymorphism) has been shown to influence teriflunomide concentrations.12 When recycling is inhibited by the administration of the bile acid binding resin cholestyramine, the half-life of teriflunomide is reduced to 1–2 days.9 Although there is a strong relationship between total and free (unbound) teriflunomide concentrations, substantial variability remains,13,14 leading to speculation that unbound concentrations may better reflect teriflunomide toxicity and efficacy.

A wide range in total teriflunomide steady-state plasma concentrations has been observed after administration of oral leflunomide11,15,16 (3–150 mg/L for 10 and 20 mg doses) and oral teriflunomide (0.1–235 mg/L9,17 for 7 and 14 mg doses). The physiological processes discussed above may contribute to this variability, and such variability could be quantitated and described using a modeling approach. Previous studies have modeled teriflunomide pharmacokinetics following leflunomide administration using one-compartment models assuming a first-order input of teriflunomide. Given the rapid conversion of leflunomide to teriflunomide, such an approach has been sufficient to investigate the effects of demographic factors and genetic variability in CYP2C19 and CYP1A2 on teriflunomide concentrations. However, it does not account for the first-pass metabolism of leflunomide to teriflunomide mediated by the CYP450 enzymes. We aimed to develop a semi-PBPK model that describes total and free teriflunomide concentrations following leflunomide or teriflunomide dosing. A structural model was developed from published mean data, and then expanded into a population model using data from patients treated with leflunomide attending our rheumatology outpatient clinic. Covariate analysis, sensitivity analysis, and simulations were used to examine the significance of factors that may account for the large interindividual concentration variability following leflunomide dosing.

METHODS

General model development

The model was developed in a stepwise manner. The initial structural model described key physiological properties of leflunomide and teriflunomide disposition based on mean literature data. This model was then refined and converted to a population model using total and free concentration determined from a cohort of 69 RA patients who were taking leflunomide. The population model was used for covariate and sensitivity analysis (Supplementary Methods).

Software

Modeling was performed using a Dell Power Edge R910 server with 4x10 core Xeon 2.26 Ghz processors and 256 GB of RAM running Windows Server 2008 R2 Enterprise with a 64-bit operating system. DigitizeIt (Boremann, Germany) was used to extract concentration data. The R Software v. 2.15.2 (R Core Team 2012) was used for data manipulation and graphical output using the ggplot2 (Wickham, 2009) and the doBy (Hojsgaard et al. 2012) packages. Model development and nonlinear mixed effect modeling of the population was performed using NONMEM v. VII Level 2.0 (Beal et al., 2009) with the Wings for NONMEM (v. 720) interface (http://wfn.sourceforge.net/) and the G95 Fortran compiler.

Mean total teriflunomide concentration data (literature data)

Mean teriflunomide concentrations following the administration of oral leflunomide in healthy young males, and oral/intravenous (IV) teriflunomide in healthy young males, and MS patients were extracted from publicly available reports to the US Food and Drug Administration (FDA) and Therapeutic Goods Administration (TGA) by Hoechst-Marion-Roussel (1998) and Sanofi-Aventis (2011, 2012).10,18,26 Mean teriflunomide data were available from one intravenous teriflunomide study, eight studies of repeat oral teriflunomide in healthy volunteers, and eight studies in MS patients. Three studies included cholestyramine “washout.”9,17 Mean teriflunomide concentrations were available from three studies following the administration of oral leflunomide in healthy volunteers.18

Total and free teriflunomide concentration data (population data)

Participants were age 18 years, diagnosed with treatment-naïve RA according to revised American College of Rheumatology Criteria,19 and were enrolled in an RA inception cohort at the Royal Adelaide Hospital (RAH) in 2000–2013. All participants were treated according to a previously published treat-to-target protocol. Leflunomide was added as a second-line therapy in those who failed to respond to maximal tolerated doses of "triple therapy" (methotrexate, sulfasalazine, and hydroxychloroquine).20 Total and free teriflunomide plasma concentrations were determined using a liquid chromatography, tandem mass spectrometry (LC-MS/MS) assay adapted from the method previously published by Rakhila et al.13 (Supplementary Methods). Participant blood samples, demographics, and genetic data were collected as previously described13,14,21–23 (Supplementary Methods), with ethics approval obtained from the Research Ethics Committee of the RAH and the University of South Australia.

Population model development

The population model was developed using the combined literature and population datasets. The model was coded using the ADVAN5 subroutine of NONMEM and fit using the FOCE-I method. Initial parameter values were based on the literature (Table 1). Importance sampling was used to assess parameter identifiability and sensitivity. Models to describe residual error of the total and free teriflunomide concentrations were explored, including additive, proportional, and combined. The model was developed in a stepwise fashion in which population mean parameters were estimated, followed by the addition of parameters describing between subject variability. Population parameter variability was assumed to follow either a log-normal or normal distribution.

Table 1

Parameter values for structural model, base model, and covariate model

	Description	Initial literature estimate	Ref.	Base model	SE%	Covariate model	SE%
OBJ		1627.814		−235.703		−249.277
Population Parameters
QCO (L/day/70 kg FFM)	Cardio output	8640	26, 27	Fixed		Fixed
QH (L/day/70 kg FFM)	Hepatic blood flow	2160	26, 27	Fixed		Fixed
QHP (L/day/70 kg FFM)	Portal hepatic blood flow	1620 [= QH × 0.75]	26, 28	Fixed		Fixed
QHA (L/day/70 kg FFM)	Arterial hepatic blood flow	540 [= QH × 0.25]	26, 28	Fixed		Fixed
QBODY (L/day/70 kg FFM)	Body blood flow	6480	26, 27	Fixed		Fixed
VLUNG (L/70 kg FFM)	Volume of lung	1	26, 27	Fixed		Fixed
VLIV (L/70 kg FFM)	Volume of liver	1.5	26, 27	Fixed		Fixed
VPORT (L/70 kg FFM)	Volume of portal blood	0.5	22, 24	Fixed		Fixed
CPROT (mg/L)	Concentration of unoccupied protein sites	4500	29	Fixed		Fixed
VBODYLEF (L/70 kg FFM)	Leflunomide volume of distribution	7.5	9, 18	Fixed		Fixed
KALEF (day−1)	Leflunomide absorption rate	24	11	Fixed		Fixed
RATEBBM (day−1)	Leflunomide blood born degradation rate constant	23.31	10	Fixed		Fixed
KA (day−1)	Teriflunomide absorption rate	39.6	11	Fixed		Fixed
KOFF (day−1)	Teriflunomide dissociation rate constant	69120	3	Fixed		Fixed
VBILE (L/70 kg FFM)	Volume of the bile	0.5	30, 31, 32	Fixed		Fixed
BILEFLOW (L/day/70 kg FFM)	Bile flow for turnover of bile volume over 24 hours	0.5	30, 31, 32	Fixed		Fixed
aQREABS (L/day/70 kg FFM)	Reabsorbed bile flow	0.475 [= BILEFLOW × 0.95]	30, 31, 32	Fixed		Fixed
aQE (L/day/70 kg FFM)	Excreted bile flow	0.025 [= BILEFLOW × 0.05]	30, 31, 32	Fixed		Fixed
CLM (L/day/70 kg FFM)		566.4	9, 10	Fixed		Fixed
MFCLL (θ1)		0.7		0.705	11.7	0.65	20
CLL (L/day/70 kg FFM)	Leflunomide to teriflunomide clearance parameter	396.5 [=CLM × MFCLL]	9, 10	399.3
CLLM (L/day/70 kg FFM)	Leflunomide to other metabolites clearance parameter	169.9 [=CLM × (1-MFCLL)]	9, 10	167.1
LGTFU (θ2)		−6.14	13	−6.15	1.6	−6.15	1.4
FU	Teriflunomide fraction unbound	0.00215	13	0.00213
CLINT (θ3) (L/day/70 kg FFM)	Teriflunomide clearance parameter	36	9	179	15.9	176	20.5
QBILE (θ4) (L/day/70 kg FFM)	Teriflunomide clearance into bile parameter	1140	9	1890	24.8	1870	35.5
VBODY (θ5) (L/70 kg FFM)	Teriflunomide volume of distribution	7.5	9, 17	8.03	14.8	7.95	11.4
COVALT(θ6)						-0.669	27.7
Population Parameter Variability
FU η (η1) (ω1; SD)				0.198	41.2	0.206	44.6
η1 shrinkage (%)				45.9		44.1
CLINT (η2) (ω2; %CV)				106.3	25.8	108.6	31.3
η2 shrinkage (%)				7.7		8.22
VBODY (η3) (ω3; %CV)				56.2	30.7	56.9	30.5
η3 shrinkage (%)				43.8		43.9
Residual Variability
Total teriflunomide proportional error
(σ1; %CV)		65.1		30.4	24.8	29.7	24.8
EPS1 shrinkage (%)				16.9		17.9
Free teriflunomide proportional error
(σ2; %CV)		66.1		39.7	19.6	39.2	19.7
EPS2 shrinkage (%)				12.0		12.2
Study data proportional error
(σ3; %CV)		339.1		64.3	24.0	64.0	24
EPS3 shrinkage (%)				6.6		6.4

When cholestyramine is administered, assume all bile contents to be eliminated; i.e., QREABS = 0, QE = BILEFLOW.

Selection criteria for the final model was guided by mechanistic plausibility, the lowest value of the objective function (OBJ) with a statistically significant (P < 0.05) improvement of model fit using the Log-Likelihood test, and the precision of parameter estimates (%RSE<50% for fixed and random effects) and visual inspection of standard diagnostic plots. The influence of total body weight, fat-free mass (FFM), and normal fat mass (NFM) were considered via allometric scaling where clearance was scaled to an exponent of 0.75 and volume to 1.24 The base model was evaluated using a visual predictive check (VPC), whereby the median, 5th and 95th percentile of the observed data were compared against the median, 5th and 95th percentiles and their empirical 95% confidence interval (CI) of 1,000 simulations of the original dataset. The statistical and observed data were binned into 10 groups and then subset into the total and free teriflunomide concentrations of the study participants, then the median and 5th and 95th percentiles were calculated without being influenced by regions of sparse data.

Simulations

The final population base model was used for simulation of total and free teriflunomide concentrations following various doses, and a local sensitivity analysis of single parameter perturbations on population total and free teriflunomide concentrations across a 1-year period was conducted (Supplementary Methods; adapted from Felmlee et al.25).

Covariate model

Missing values for continuous and categorical covariates were imputed. Continuous covariates were subjected to a last observation carried forward/previous observation carried backward approach. For patients missing all continuous data, the data were imputed with the median, while categorical variables were imputed with the mode. Covariates on relevant parameters were screened according to biological plausibility and prior knowledge of the factors considered to influence teriflunomide concentrations (Table 2).

Table 2

Biologically guided covariates screening plan

Parameter	Covariates screened
VBODY	Age, sex, and albumin
LGTFU	Age, sex, albumin, and creatinine clearance (CrCl)
CLINT	Age, sex, CrCl, bilirubin, ALT, AST, and methotrexate dose.
QBILE	Bilirubin, sulfasalazine dose
MFCLL	Smoking status, ALT, AST, CYP1A2*1F allele status (rs762551; AA versus CC and AC), CYP2C19 phenotype (rs4244285 and rs12248560; poor and intermediate metabolizers versus extensive metabolizer vs. ultra-rapid metabolizers)21

Continuous covariates were coded via a power model referenced to the median of the observed data, while the effects of categorical variants were assessed through a binary relationship. Concomitant drug effects were investigated via a linear function with their dose. The potential inclusion of covariates on the base model were selected based upon a significant decrease in OBJ (P < 0.01; SE < 50%). The final model was developed via forward inclusion/backward deletion, where each covariate was sequentially added, starting with the covariate that caused the largest drop in the OBJ, and the following covariate(s) were only retained if the selection criteria were met. The final covariate model was further subjected to a backward elimination procedure (P < 0.001 level).

RESULTS

Structural model development from literature data

A number of key physiological steps were identified for the semi-PBPK model including 1) the conversion of leflunomide to teriflunomide, 2) hepatic clearance of unbound teriflunomide, and 3) enterohepatic recycling of teriflunomide. The final compartmental structure is shown in Figure 1. Parameters such as cardiac output (QCO), hepatic blood flow (QH), lung volume (VLUNG), and hepatic volume (VLIV) were common to both leflunomide and teriflunomide (Figure 1), and were fixed to average physiological values.26,27 To accommodate drug absorption, hepatic-portal (75%) and hepatic-arterial (25%) blood flows were separated, with a hepatic-portal volume of 25 mL per 100 g of liver (i.e., QHp = 0.75*QH; volume fixed to 500 ml).26,28 The absorption rate of leflunomide (KAL) was initially based on the value used by Grabar.11 Conversion of leflunomide to its metabolites was controlled by CLM, and the initial estimate was based on in vitro data.10 CLL controls the hepatic conversion of leflunomide to teriflunomide, which was initially estimated to be 70% of CLM, and the remaining 30% representing the conversion of leflunomide to other metabolites.9 Teriflunomide was also formed from leflunomide by spontaneous blood-borne degradation (rate constant RATEBBM) based on in vitro data.10 In the absence of any suitable data, given similarities in structures and protein binding, the apparent volume of the body compartment (VBODYLEF) for leflunomide was assumed to be the same as teriflunomide's apparent volume of the body compartment (VBODY) (7.5 L; 11 L (Vd9,18) minus the volumes of the lung, liver, portal blood, and bile).9,18

Figure 1

Structural model of leflunomide and teriflunomide pharmacokinetics.

Representation of unbound drug is considered important for highly bound drugs, such as teriflunomide.3 This model controls drug entry into the liver via blood flow (QH), it is then separated into bound and unbound as controlled by the fraction unbound (FU), which was assumed to be the same as within plasma (0.00215).13 This dissociation equilibrium (KD; Eq. 1 was maintained by the dissociation (KOFF) and association (KON) rate constants. A “well-stirred” model was achieved by fixing KOFF to a nonlimiting value (2,880 h−1)3. KON was a function of FU and plasma protein concentration (CPROT29) (Eq. 2).

The enterohepatic recycling pathway replicated known physiology with 500 ml of bile produced daily, 95% of which was reabsorbed while the remaining 5% was excreted and may contain teriflunomide.30–32 Cholestyramine “washout” was assumed to completely stop enterohepatic recycling, as this made CLINT and QBILE identifiable via simulation.9 The initial parameter estimates of the model based on the literature are represented in Table 1. Respectively CLINT and QBILE represent a teriflunomide clearance parameter and a teriflunomide clearance into bile parameter, while MFCLL represents a fractional constant controlling the ratio of CLL and CLLM, as defined in Figure 1.

The first-order conditional estimation with interaction (FOCE-I) with ADVAN5 was initially used to fit the model against the mean (total) teriflunomide concentrations (272 observations) extracted from the literature (described in Supplementary Results).9,17,18 No between-study parameter random effects were estimated. Intrastudy variability was represented by a residual error model weighted by the inverse of the square root of the number of patients within each study (Eq. 3).

The estimation of the values of VBODY, CLINT, QBILE, and MFCLL significantly improved the fit over literature-based estimates (Table 3). Estimated concentrations from the best model (RUN8) showed a good fit after either leflunomide or teriflunomide dosing (Figure 2), indicating that the conversion process was described well. Similarly, the long half-life of teriflunomide significantly decreased when enterohepatic recycling was turned off (cholestyramine “washout”). The structural model therefore described the key physiological features of leflunomide and teriflunomide disposition.

Table 3

Parameter values for structural model when fitted to mean teriflunomide concentrations

Run	Obj	MFCLL	VBODY (L)	CLINT (L/h)	QBILE (L/h)
RUN1_LITERATURE_BASED_ESTIMATES	690.635	0.7	7.5	1.5	47.5
RUN8_NOFIX_MFCLL_VBODY_CLINT_QBILE	229.557	0.645	5.55	3.76	50.4

Figure 2

Mean teriflunomide concentrations extracted from the literature are represented in gray dots, while the blue lines indicate the population predicted curves from RUN8. Curves 1 to 17 represent mean teriflunomide concentrations following the oral/IV administration of teriflunomide in healthy volunteers and MS patients, three studies of which contained cholestyramine “washout.” Curves 18 to 20 represent mean teriflunomide concentrations following oral leflunomide (LEF) administration.

Population model development from literature and patient data

Total and free teriflunomide concentration data from 69 RA patients receiving leflunomide (156 samples; demographics in Supplementary Results) were best described by proportional residual error models, while the extracted mean (total) teriflunomide concentrations were represented as described above. The data were best described by the estimation of the population values of MFCLL, CLINT, VBODY, QBILE, and LGTFU with population parameter variability following a log-normal distribution for CLINT and VBODY and a normal distribution for the LGTFU parameter. LGTFU is a logit constraint parameter (Eq. 5) for the FU (Eq. 4),

Allometric scaling for weight, FFM,33 and NFM34 significantly reduced the OBJ value, and inclusion of FFM resulted in the greatest reduction in OBJ, and was thus accepted as the base model. The base model had an OBJ of −235.703 and precise parameter estimates (Table 1).

Diagnostic plots of total and free teriflunomide concentration predictions are represented in the Supplementary Data and were considered compatible with an unbiased model. The VPCs show the model provided a good description of the lower range and median total and free teriflunomide concentrations, as represented by a good overlay of the median and 5th percentiles (Figure 3). However, the model may tend to overpredict the upper 95% percentile free and total concentrations, although the observed data's 95th percentile remains within the wide empirical 95% CIs, suggesting between-subject variability is high and the amount of data informing the high concentrations is limited.

Figure 3

Visual predictive checks of the base model for total (above) and free (below) teriflunomide concentrations. The observed teriflunomide concentrations are represented by black circles, a black solid line (median), and the black dashed lines (5th and 95th percentiles). The simulated total and free teriflunomide concentrations are represented by the red (empirical 95% CI of median) and the blue shaded areas (empirical 95% CIs of 5th and 95th percentiles). The solid red line indicates the median of simulated concentrations, and the solid blue lines are the 5th and 95th percentiles of the simulated data.

Model simulations

Steady-state total and free teriflunomide concentrations were simulated for a population (n = 1,000) with a FFM of 44.55 kg receiving various leflunomide doses (5–20 mg/day) across a 1-year period. The 95% prediction interval of teriflunomide concentrations largely overlap throughout the 1-year period, indicating the extent of between subject variability (Figure 4).

Figure 4

Total (above) and free (below) teriflunomide concentrations for a simulated population from the base model. The shaded area indicates the 95% CI of the simulated data for each dose, with the central solid line indicating the median of simulated concentrations.

Sensitivity analysis

The population predicted and sensitivity indices of total and free teriflunomide concentrations following QBILE, CLINT, MFCLL, and FU parameter perturbation with a standard 20 mg daily dose are represented in the Supplementary Data. Increasing QBILE and CLINT decreased total and free teriflunomide concentrations, which is the converse to MFCLL perturbations. Decreased FU resulted in higher total teriflunomide concentrations. An increase in FU results in a transient increase in free concentration; however, after equilibration unbound steady-state concentrations returned to the prior level.35

Covariate analysis

Of the covariates screened, only the effect of alanine aminotransferase (ALT) on CLINT, aspartate aminotransferase (AST) on MFCLL, and sex on FU met the inclusion criteria, while there were no identifiable genetic influences. The influence of ALT on CLINT caused the greatest reduction in OBJ at 13.6 units; neither the forward inclusion of AST on MFCLL or sex on FU significantly reduced OBJ. The parameter estimates for the final covariate model are presented in Table 1, and the relationship between CLINT with FFM and ALT concentrations is represented by Eq. 6.

The typical value of CLINT for an individual with an FFM of 70 kg and ALT of 26 units/L is 176 L/day, with a between-subject variability CI of 108.6%. Therefore, as ALT concentrations increase the CLINT appears to decrease in a typical individual (e.g., if typical ALT doubles to 52 units/L, CLINT decreases to 110 L/day).

DISCUSSION

The population approach was used to describe free teriflunomide concentrations. This was accomplished by structuring the model around a physiological process in which the liver handles unbound drug.3 The model allows for both total and free teriflunomide concentrations to be described, which could lead to better studies investigating the total and/or free teriflunomide concentration relationship to response/toxicity.13,36 Previous small studies (23–71 participants) have indicated a relationship between total teriflunomide concentrations and effect in RA patients.11,15,37 However the reported response thresholds have been highly variable, with Van Roon et al.15 indicating improved response in those with total concentrations over 16 mg/L while Chan et al.37 recommended concentrations >50 mg/L. As such, therapeutic drug monitoring (TDM) of total concentrations is not currently used to guide leflunomide dosing. With the model described here, subsequent studies could use the population approach to consider the effects of total and free teriflunomide concentrations separately, and in combination with pharmacodynamic models account for demographic and pharmacogenetic factors in a physiologically appropriate fashion.

Simulations of normal clinical dosages were performed. It is acknowledged that the population variability of the CLINT and VBODY parameters is high, and the empirical 95% CIs of the 95th percentiles are wide in the VPC. As such, the CIs of the simulated concentrations should be considered with caution. However, it does appear that some individuals receiving 20 mg daily may have steady-state total and free teriflunomide concentrations lower than some individuals receiving 10 mg on alternate days. This apparent disparity highlights the currently restrictive nature of a 20 mg maximum dosing, and that TDM may be clinically useful. Conversely, to the wide empirical 95% CIs of the 95th percentiles, the empirical 95% CI of the median percentile is tighter and the median prediction was in good agreement with the observed data. This is not surprising, due to the large quantity of mean data informing the model, extracted from the literature. Therefore, mean prediction appears to be accurate, particularly within the dosage range of the study participants, hence the sensitivity analysis which aimed to inform future studies of the potential importance of covariate influences on certain parameters would appear appropriate.

In this study, a local sensitivity analysis (one at a time) of the QBILE, CLINT, MFCLL, and the FU parameters was performed, indicating the model to be least sensitive to perturbations of the QBILE parameter, and more sensitive to changes in CLINT and MFCLL (Supplementary Results). The primary aim of the analysis was to provide an indication of the most influential parameters on the model and the predictive power of the model for future studies. This indicated that when examining pharmacogenetic variables influencing the QBILE, CLINT, and MFCLL parameters, a more substantial effect (or a larger cohort of individuals) would be required to quantify the influence of the CLINT or MFCLL parameters compared to the QBILE. Confounding this finding is the lack of data informing the QBILE parameter within this study, where the estimation of the parameter was based on three studies with mean total teriflunomide concentrations after the administration of the cholestyramine washout.9,17 Future studies will ideally include individual participant data following the administration of cholestyramine to inform this physiological process and quantify the population variability. This lack of data would likewise have limited the ability to examine the potentially important effects of ABCG2 genotype on the QBILE parameter. FU perturbation exemplified the physiological nature of the model, where increases in the FU result in lower steady-state total concentrations, while there is a transient increase in free concentrations before returning to the original steady-state, in keeping with the concept that steady-state concentrations of the free drug do not depend on protein binding after oral administration of a hepatically cleared drug.35

In contrast to previous studies describing teriflunomide pharmacokinetics using one-compartment models with first-order inputs following leflunomide administration,11,37–40 we chose to represent leflunomide and teriflunomide pharmacokinetics with a semi-PBPK model. The complexity of the developed model required that we pool the mean teriflunomide concentrations from the literature to stabilize the parameter estimation, which included mean curves from healthy individuals as well as RA and MS patients. Using these potentially heterogeneous data may be a source of bias for the parameter estimates. However, previous studies show no appreciable differences in teriflunomide pharmacokinetics in MS compared to healthy subjects, and we assumed this to also apply to individuals with RA.9,17 The semi-PBPK approach was chosen based on the predicted use of the model within future studies, while the full PBPK approach appeared unrealistic, as emphasized by this study requiring the fixation of the enterohepatic process and the circulatory component of the model. Hence, the semiphysiological approach appears acceptable for the intended use of the model, which was to independently explore potential covariate effects on the conversion of leflunomide to teriflunomide, the teriflunomide enterohepatic recycling process, and protein binding effects. Ideally, however, the model would have sufficient data to describe (with acceptable precision) the population variability of the QBILE and MFCLL parameters. This is difficult to accomplish at present, as leflunomide has not been detected in plasma, while individual patient data describing teriflunomide concentrations after cholestyramine administration is not publicly available.9,17

The most comprehensive analysis to date of total teriflunomide concentrations used a one-compartment model and data from 435 patients receiving leflunomide.38 Age, sex, and body size appeared to influence teriflunomide concentrations, although the effects were not clinically significant. Studies with smaller cohorts of 71 and 23 patients reached similar conclusions.37,39 Recently, attention has turned to pharmacogenetic variants such as the 421C>A single nucleotide polymorphism in the ABCG2 transporter gene. Among 24 individuals, those who carried the A allele had 41% lower “clearance” and a 30% higher teriflunomide Cmax after a single dose.12 Grabar et al.11 used a one-compartment population model and data from 71 patients, showing the loss-of-function CYP2C19*2 allele was associated with a 71% increase in oral teriflunomide “clearance,” although this process likely represents altered first-pass metabolism. We used a similar sized cohort in our multivariate genetic analysis but did not identify an influence of ABCG2 or CYP2C19*2 genotype on total or free teriflunomide concentrations. Given that we used a relatively small dataset, it is difficult to conclude that these variants have no effect. But as the population variability of the CLINT and VBODY parameters is large and no influence was seen, it appears unlikely they will be of use in guiding clinical dosing, although population variability of QBILE was not accounted for in this study. Allometric scaling of the model with FFM (calculated using total body weight, height, and sex33) significantly improved the fit of the model, which is not surprising due to the apparent small volume of distribution of teriflunomide as evidenced by high plasma protein binding. Likewise, accounting for ALT on CLINT revealed a significantly improved model fit, although like previous models analyzing demographic and laboratory covariates, the clinical significance was minimal.

This study has conducted the most extensive analysis of genetic covariates on leflunomide kinetics to date; however, we concede that other factors not assessed may significantly affect between- and within-subject variability. Notably, we assessed the influence of methotrexate and sulfasalazine, the most commonly used disease-modifying antirheumatic drugs (DMARDs) within our RA cohort; however, many other medications were used, which may be important if the drug has enzyme-modulating effects or causes a rise in ALT concentrations. Furthermore, cytokines promote synovitis central to RA activity, and elevated levels have been associated with altered absorption and metabolism of drugs, which may occur through the downregulation of enzyme expression, particularly CYPs.41–43 Future studies should therefore consider the impact of cytokines and disease activity on teriflunomide concentrations, which may impact leflunomide conversion or teriflunomide reabsorption.

For the first time, the physiological processes that influence leflunomide and teriflunomide disposition have been represented in a semi-PBPK population model. This is considered important, as genetic variables have been shown to influence the conversion of leflunomide to teriflunomide, as well as the teriflunomide enterohepatic recycling process, while teriflunomide is highly bound within the plasma and small changes in its protein binding may have dramatic effects on it disposition and clearance. Furthermore, the model is capable of describing teriflunomide disposition following oral dosing of either leflunomide or the newly released teriflunomide, and it therefore appears to be an important step in quantifying the similarities and differences between these compounds. We anticipate this physiological fidelity may be important for future studies for optimizing the use of leflunomide and teriflunomide in RA, MS, and other diseases.