This article is mentioned in:

Introduction

Systemic-onset juvenile idiopathic arthritis (SOJIA) is a serious type of juvenile arthritis (1). It is driven by continuous activation of innate immune pathways producing pro-inflammatory cytokines (2).

Treatment with tacrolimus may suppress early activation of interleukin (IL)-2 gene transcription and inhibit the production of tumor necrosis factor-α (TNF-α), IL-1β and IL-6 during T-cell activation (3,4). Tacrolimus treatment in patients with SOJIA has been reported previously (5–8). However, as an immunosuppressive agent, the therapeutic range of tacrolimus is narrow, with considerable inter- and intra-individual variability (9,10).

Population pharmacokinetics (PPK) may be used to acquire PK information from sparse population data. Furthermore, PPK may differentiate between inter-individual and intra-individual variability and has considerable power to uncover the effects of confounding factors on PK behavior and to ensure that a treatment is suitable for personalized clinical therapy (11). In previous studies, PPK models of tacrolimus were set up among different populations (12–24). However, the PPK model of tacrolimus treatment in patients with SOJIA has remained to be established. The aim of the present study was to set up a tacrolimus PPK model in patients with SOJIA and to formulate initial dosage recommendations for personalized treatment.

Materials and methods

Patient data

Data from Chinese patients with SOJIA who attended the Children's Hospital of Fudan University (Shanghai, China) between January 2014 and December 2017 were retrospectively collected. The inclusion criteria were as follows: i) Pediatric patients with SOJIA (aged <16 years); ii) treatment with tacrolimus; iii) tacrolimus concentrations were routinely tested by therapeutic drug monitoring (TDM). Subjects with other concurrent serious clinical conditions were excluded (e.g. liver, kidney or bone marrow transplant). A schematic depicting the recruitment of the patients is provided in Fig. 1. Relevant clinical information and data on drug concentrations were gathered from medical records and TDM records, respectively. Demographic data of the patients and concomitant drugs were used as potential covariates to be analyzed in the current PPK model. The present study was a retrospective study and was approved by the ethics committee of the Children's Hospital of Fudan University (Shanghai, China) without the requirement for written informed consent.

Figure 1. Schematic of patient recruitment. SOJIA, systemic-onset juvenile idiopathic arthritis; TDM, therapeutic drug monitoring; NONMEM, nonlinear mixed-effects model.

Drug administration and analytical method

Oral tacrolimus dose adjustment was based on safety and effectiveness, along with the drug trough concentration from TDM. The whole-blood concentration of tacrolimus was analyzed using the Emit® 2000 Tacrolimus assay (Siemens Healthcare Diagnostics Inc.) according to the manufacturer's protocol.

PPK modeling

Patient data were analyzed using the nonlinear mixed-effects model (NONMEM, edition 7; ICON Development Solutions). PK parameters and their variability were estimated using the first-order conditional estimation method with interaction. The absorption phase was described using a first-order absorption and elimination one-compartment model. The PK parameters included apparent oral clearance (CL/F) and apparent volume of distribution (V/F). F was the bioavailability. The absorption rate constant (Ka) was set at 4.48/h (15).

Random-effects model

Inter-individual variability was analyzed using an exponential error model, as presented in equation i:

Pj is the value of the individual parameter. TV(P) is the parameter of the typical value. The individual deviation is represented by ηj, which is a symmetrically distributed, zero-mean random variable with variance terms.

The random residual variability was described using equation ii:

Y is the concentration observed and IPRED represents the individual predicted concentration. The variation is represented by ε, which is a symmetrically distributed, zero-mean random variable with variance terms.

Covariate model

Weight and PK parameters were modeled using equation iii:

Pj represents the PK parameter of the j-th individual, Weightj is the Weight of the j-th individual and Pnorm is the parameter of an individual with a normal Weight (Weightnorm) of 70 kg. The COE is the allometric coefficient: 0.75 for the CL/F and 1 for the V/F (25).

Continuous covariates and categorical covariates were used to describe the correlation between PK parameters using equations iv and v:

Pj and TV(P) are the individual parameter value and typical parameter value, respectively. Covj is the covariate of the j-th individual and θ is the parameter to be estimated. Covmed is the population median for the covariate.

This stepwise protocol was used to build the covariate model. The likelihood ratio test was used to compare hierarchical models. The covariate model was established in a stepwise manner, using the forward inclusion, backward elimination method (11,21,24,26,27). Changes in objective function values (OFV) were performed using covariate inclusions and a decrease of OFV >3.84 (P<0.05) was considered sufficient for inclusion in the base model (11,21,24,26,27). After establishing a full regression model, the model was further assessed by eliminating covariates from each PK parameter one-by-one to obtain the final model. An increase in OFV >6.64 (P<0.01) was considered sufficient for significance in the final model (11,21,24,26,27). This statistical method and its description have been published in numerous similar studies and may be considered as a fixed and applicable statistical method for PPK analysis (11,21,24,26,27).

Model validation

The stability and reliability of the final parameter estimates were evaluated using the internal validation method of bootstrap, which was produced using repeated random sampling, with replacement from the original data. The process was performed using the Wings package for NONMEM software and was repeated 1,000 times with different random draws. The median values and 2.5–97.5% percentile parameters from the bootstrap results were compared with the final PK parameters. Visual inspection of routine diagnostic plots along with prediction-corrected visual predictive check (VPC) plots were used to assess the final model.

Simulation of initial dosage recommendations

Monte Carlo simulation is an approach used to determine probability of target (28) and has been applied to determine the most suitable drug administration (22,26). In the present study, it was used to investigate the influence of covariates on the probability to achieve the target concentrations. A previous study reported that for safety reasons, the lower concentration for tacrolimus treatment for SOJIA was 1.7 ng/ml and the upper concentration was 5 ng/ml (7). Therefore, the probability to achieve 1.7 and 5 ng/ml concentration thresholds based on the established model without the combination with other drugs was estimated. Simulation was performed for each of the nine weight groups (5, 15, 20, 25, 30, 35, 40, 45 and 50 kg) and four dosing regimens [0.5 mg once every 24 h (q24h), 0.5 mg q12h, 1/0.5 mg q24h and 1 mg q12h] using 1,000 virtual patients with SOJIA.

Results

Data collection

Data from 17 Chinese patients with SOJIA (8 males and 9 females), aged 9.50 (3.20–14.60) years, were collected to build the population model. A total of 86 concentrations in the range of 1.3 to 9.2 ng/ml were used. Patient information and drug combinations are presented in Tables I and II.

Table I. Demographic data of the patients.

Table I.

Demographic data of the patients.

Characteristic	Mean ± SD	Median (range)
Age (years)	8.23±3.30	9.50 (3.20–14.60)
Weight (kg)	29.83±10.66	33.60 (13.50–46.00)
Duration of treatment with tacrolimus (days)	190.59±169.57	67.00 (5.00–535.00)
Daily dose of tacrolimus (mg)	1.66±0.71	1.50 (1.00–4.00)
Alanine transaminase (IU/l)	28.56±60.24	12.00 (2.00–458.00)
Aspartate transaminase (IU/l)	21.01±26.82	15.00 (7.00–235.00)
Creatinine (µmol/l)	35.36±10.32	35.50 (19.00–59.00)
Hematocrit (%)	37.35±3.43	37.45 (30.30–44.40)
Hemoglobin (g/l)	121.59±13.23	124.00 (90.20–152.00)
Mean corpuscular hemoglobin (pg)	27.13±2.01	27.00 (21.00–31.00)
Mean corpuscular hemoglobin concentration (g/l)	324.58±15.22	323.00 (285.00–358.00)

[i] SD, standard deviation.

Table II. Drug combinations with tacrolimus.

Table II.

Drug combinations with tacrolimus.

Drug/category	N
Ranitidine
0	16
1	1
Hydroxychloroquine
0	16
1	1
Ceftazidime
0	16
1	1
Cefmetazole
0	16
1	1
Ceftriaxone
0	16
1	1
Cefprozil
0	13
1	4
Cefixime
0	14
1	3
Cefdinir
0	11
1	6
Azithromycin
0	15
1	2
Methylprednisolone
0	12
1	5
Mycophenolate mofetil
0	16
1	1
Prednisone
0	2
1	15
Oxcarbazepine
0	16
1	1
Levetiracetam
0	16
1	1
Methotrexate
0	9
1	8
Omeprazole
0	11
1	6
Diltiazem
0	14
1	3
Felodipine
0	16
1	1
Montelukast
0	16
1	1
Aspirin
0	15
1	2
Loratadine
0	12
1	5

[i] Categories: 0, without drug; 1, with drug. N, number of patients.

Modeling

The first-order absorption and elimination one-compartment model was identified to fit the dataset. PK parameters from the final covariate models were as follows in equations vi and vii:

When patients were co-administered omeprazole, loratadine or diltiazem, the value of each was 1; otherwise, the value was 0. All weights in equations vi and vii were measured in kg.

Validation

The visual inspection of routine diagnostic plots is presented in Fig. 2 and the parameter estimates from the final model and bootstrap validation are presented in Table III. From 1,000 bootstrap runs, 988 runs were successfully minimized. Using Table III, the parameter estimate median values of bootstraps were found to be similar to the respective values determined with the final model, indicating that the final PPK model was accurate and reliable. The VPC plots for the final model (Fig. 3) demonstrate that most of the measured concentration data were included in the 95% prediction intervals of the simulation data, suggesting that the final PPK model is able to predict concentrations effectively.

Figure 2. Visual inspection of routine diagnostic plots. (A) Observations vs. population predictions of tacrolimus blood concentrations. (B) Observations vs. individual predictions of tacrolimus blood concentrations. For (A) and (B), black solid lines represent the line of unity (the y=x line, where predications match observed values), and the red smooth line represents the trend of the data. Hence, the closer the red smooth line is to the black solid line, the more predictive the model is. (C) |iWRES| vs. individual predictions. |iWRES| is the difference between the individual predictions and the observed tacrolimus blood concentrations. (D) weighted residuals vs. the time of tacrolimus blood concentration. black solid line represents the line of unity (the y=0 line). For (C) and (D), the red smooth line represents the trend of the data, therefore the closer the red smooth line is to the line of unity (the y=0 line), the more predictive the model is. |iWRES|, individual weighted residuals.

Figure 3. Prediction-corrected visual predictive check for the final model. The middle solid line represents the median of the prediction-corrected concentrations. The lower and upper dashed lines are the 2.5 and 97.5th percentiles of the prediction-corrected concentrations, respectively. The data-points indicate the measured concentrations. Theoretically, inclusion of the measured concentration in the 95% confidence interval of predicted values indicates good predictability of the model. Pink areas indicate the confidence interval of the middle solid line, purple areas indicate the confidence interval of the lower and upper dashed lines.

Table III. Parameter estimates of final model and bootstrap validation.

Table III.

Parameter estimates of final model and bootstrap validation.

			Bootstrap
Parameter	Estimate	SE (%)	Median	95% CI	Bias (%)
CL/F (l/h)	29.700	9.300	29.800	(24.300, 36.400)	0.340
V/F (l)	1120.000	27.900	1120.000	(604.000, 2188.000)	0
Ka (1/h)	4.480 (fixed)	–	–	–	–
θomeprazole	−0.362	16.800	−0.371	(−0.499, −0.192)	2.490
θloratadine	−0.322	23.800	−0.326	(−0.462, −0.081)	1.240
θdiltiazem	−0.307	34.200	−0.307	(−0.454, −0.006)	0
ωCL/F	0.265	18.400	0.243	(0.129, 0.352)	−8.300
σ1	1.229	5.100	1.200	(1.040, 1.326)	−2.360

[i] 95% CI was displayed as the 2.5th, 97.5th percentile of bootstrap estimates. SE, standard error; CL/F, apparent oral clearance (l/h); V/F, apparent volume of distribution (l); Ka, absorption rate constant (1/h); θ_omeprazole, coefficient of omeprazole; θ_loratadine, coefficient of loratadine; θ_diltiazem, coefficient of diltiazem; ω_CL/F, inter-individual variability of CL/F; σ₁, residual variability, additive error. Bias, prediction error, calculated as Bias=(Median-Estimate)/Estimate ×100%.

Simulation

In the present study, the initial tacrolimus dose without drug combination was predicted. In clinical practice, combination therapy is not common at the time of initial administration. Therefore, the probability to achieve the target concentrations based on the established model without any combined drugs was estimated. The predicted median, along with the 2.5–97.5% percentile parameters and the probability of achieving the target concentration were presented in Table IV. According to the simulation dataset, the 0.5 mg q24h regimen appeared to be most suitable for pediatric patients with 5 kg body weight, the 0.5 mg q12h regimen was appropriate for patients with 15–25 kg body weight, the 1/0.5 mg q24 h regimen was most suitable for subjects with a body weight of 26–35 kg and the 1 mg q12h regimen was fit for the group with a body weight of 36–50 kg.

Table IV. Predicted median tacrolimus concentration (ng/ml), 95% CI and probability (%) of achieving the target concentration with respect to body weight for different dosing regimens.

Table IV.

Predicted median tacrolimus concentration (ng/ml), 95% CI and probability (%) of achieving the target concentration with respect to body weight for different dosing regimens.

	Body weight (kg)
Regimen	5	15	20	25	30	35	40	45	50
0.5 mg q24h	2.57 (1.00–5.82)	1.33 (0.59–2.74)	1.11 (0.51–2.23)	0.96 (0.45–1.90)	0.86 (0.41–1.66)	0.77 (0.37–1.48)	0.71 (0.34–1.34)	0.65 (0.32–1.23)	0.61 (0.30–1.13)
	77.2%	22.5%	12.0%	5.3%	1.9%	0.9%	0%	0%	0%
0.5 mg q12h	7.59 (3.54–14.64)	3.57 (1.77–6.52)	2.92 (1.47–5.24)	2.49 (1.27–4.42)	2.19 (1.12–3.85)	1.96 (1.01–3.41)	1.78 (0.93–3.08)	1.63 (0.85–2.81)	1.51 (0.80–2.58)
	12.2%	82.7%	89.7%	88.0%	79.4%	67.7%	55.9%	45.4%	35.3%
1/0.5 mg q24h	10.25 (4.45–20.66)	4.94 (2.32–9.33)	4.06 (1.94–7.54)	3.48 (1.69–6.37)	3.07 (1.50–5.55)	2.75 (1.36–4.93)	2.51 (1.25–4.45)	2.31 (1.16–4.07)	2.14 (1.08–3.75)
	4.9%	50.9%	71.6%	83.1%	87.9%	89.1%	87.5%	83.3%	77.1%
1 mg q12h	15.18 (7.08–29.27)	7.13 (3.54–13.03)	5.83 (2.94–10.49)	4.98 (2.53–8.85)	4.37 (2.25–7.69)	3.92 (2.03–6.83)	3.56 (1.85–6.15)	3.27 (1.71–5.61)	3.03 (1.59–5.17)
	0.1%	13.7%	31.7%	50.5%	65.3%	77.3%	89.5%	89.8%	92.1%

[i] 95% CI was displayed as the 2.5–97.5th percentile. The frames indicate the optimal dose for the respective body weight. q24h, once every 24 h.

Discussion

To control SOJIA disease, a number of patients require long-term corticosteroid treatment (5). However, prolonged and repeated steroid treatment increases the risk of adverse reactions, including obesity, cushingoid appearance, hypertension, growth retardation, osteoporosis, infections and psychological problems (29). Thus, a safe and effective therapeutic method to treat patients with SOJIA remains to be explored (30).

Recent studies revealed the beneficial impact of suppressing IL-6 and other pathogenic pro-inflammatory cytokines for controlling SOJIA (30–32). Tacrolimus potently suppresses the production of TNF-α, IL-1β and IL-6 through T-cell activation (4,33), and it has therefore been administered to patients with SOJIA (5–7).

However, the therapeutic range of tacrolimus is narrow, with considerable inter-individual and intra-individual variability (9,10). Thus, it is necessary to build a tacrolimus PPK model for patients with SOJIA and to formulate initial dosage recommendations for personalized treatment.

To the best of our knowledge, the present study was the first to provide a PPK model of tacrolimus for patients with SOJIA. The PPK model was established for SOJIA patients by using a population modeling method. The approach was necessary, as logistic and ethical restrictions prohibit excessive blood sampling when studying pediatric patients (34). The tacrolimus PPK model is able to predict the PK process in patients with SOJIA and it therefore has important clinical value.

In the present study, the first-order absorption and elimination one-compartment model fitted the dataset, as all of the tacrolimus concentrations were trough concentrations and the Ka was fixed at a value from the literature of 4.48/h (15). It was not possible to estimate the area under the curve, minimum concentration and maximum concentration of tacrolimus, as the drug was orally administered and tacrolimus concentration data were insufficient. The typical CL/F and V/F values of the final tacrolimus PPK model were 29.7 l/h and 1,120 l. In the present PPK model, drug combinations were used as categorical variables. The present study also tested the influence of the following various covariates on different parameters: Weight, omeprazole, loratadine and diltiazem on CL/F, as well as Weight on V/F. Numerous studies have determined a non-linear association between drug clearance and body weight in pediatric patients, and it may be well described with allometric scaling using a coefficient of 0.75 for clearance and 1 for volume (25,26,35,36). Body weight is the most important predictor of clearance and volume in pediatric patients with maturation of elimination processes (35), and is also considered to be the primary factor determining clearance and volume based on the theory explaining the link between mass, function and structure; this theory is valid across numerous orders of magnitude of body weight (37). Important factors that also impacted tacrolimus clearance were omeprazole, loratadine and diltiazem, possibly due to tacrolimus being a substrate of the cytochrome P450 3A (CYP3A) enzyme (38), and omeprazole and diltiazem inhibit CYP3A activity. In addition, loratadine is a CYP3A substrate that is able to compete with tacrolimus for the binding site on the enzyme and lead to a decrease in tacrolimus clearance. Thus, concomitant medication with omeprazole, loratadine or diltiazem may reduce tacrolimus clearance in patients with SOJIA.

In terms of model application, Monte Carlo simulations based on the established model were used to investigate the influence of covariates on the probability to achieve the target concentration. The probability to achieve 1.7 and 5 ng/ml concentration thresholds based on the established model without drug combinations was estimated. In addition, it appears that IL-6 is a marker of pharmacodynamics of tacrolimus (8). In future studies by our group, a population pharmacodynamics model will be built to analyze the association between drug exposure and IL-6.

In conclusion, to the best of our knowledge, the present study provided the first PPK model of tacrolimus in patients with SOJIA, and may be used for precision therapy in pediatric patients. A large external evaluation of this model will be performed in future studies.

Acknowledgements

Not applicable.

Funding

The present study was supported by the Clinical Pharmacy Key Specialty Construction Project of Shanghai (grant no. YZ2017/5), Important Weak Subject Construction Project of Shanghai (grant no. 2016ZB0305), the Scientific Research Project of Science and Technology Commission of Shanghai Municipality (grant no. 18DZ1910604) and the Fudan University Hospital Management Construction Project (Fudan medical administration 2018).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

ZL and HX conceived and designed the study. DW and XC collected and analyzed the data. DW and XC wrote the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The study was approved by the Research Ethics Committee of the Children's Hospital of Fudan University (Shanghai, China).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References