A total of 293 plasma concentrations from 48 patients with CPP were analyzed. The demographic and clinical characteristics of the patient population included in the pharmacokinetic analysis are summarized in Table 1.

A one-compartment disposition model with sequential first-order absorption processes (immediate and delayed) and first-order elimination best described the data. The schematic of the structural pharmacokinetic model is illustrated in Fig. 1. Implementation of a delayed absorption phase via lag time produced poorer fits compared with the transit compartment, based on a graphical assessment of their ability to describe the data and OFV (42 points higher). The inclusion of immediate zero-order absorption improved the goodness of fit; however, parameters such as the duration of input for the zero-order absorption process could not be estimated well. The final model was parameterized in terms of absorption rate constants for the immediate and delayed absorption processes (Ka1 and Ka2, respectively), mean transit time (MTT), number of transit compartments (N), apparent clearance (CL/F), and apparent volume of distribution (Vd/F), where F is the bioavailability. Mean transit time is the average time leuprolide spends in each transit compartment (Fig. 1). The number of transit compartments was fixed due to insufficient data to estimate the parameter with good precision. The number of compartments was fixed at three based on a successful model run that adequately described the absorption phase in the model diagnostics but failed to describe other phases. The final model included between-subject variability on Ka1, Ka2, clearance (CL), volume (V), and MTT, as well as a proportional residual error model. The structural model of the final model is shown in the following equations:

$$K_ = \theta_ \times e^ }}, \,\text \, \eta_ \sim N\left( ^ } \right)$$

$$K_ = \theta_ \times e^ }}, \, \,\text \, \eta_ \sim N\left( \omega_^ } \right)$$

Illustration of the structure of the final pharmacokinetic model. CL clearance, F1 fraction of dose absorbed via immediate first-order process, F2 fraction of dose absorbed via delayed first-order process, Ka1 absorption rate constants for the immediate absorption processes, Ka2 absorption rate constants for the delayed absorption processes, KTR transit rate between compartments

$$}_=_\times ^_}, \, \,\text \, _\sim N\left(0, _}^\right)$$

$$_=_\times ^_}, \, \,\text \, _\sim N\left(0, _}^\right)$$

$$}_=_\times ^_}, \, \,\text \, _\sim N\left(0, _}^\right)$$

Differential equation (1):

$$\frac}A_}1j}} }}}t_ }} = - K_ \times A_}1j}} $$

Differential equation (2):

$$\frac}A_}2j}} }}}t_ }} = } \times F_ \times K_}}} \times \frac}}} \times t_}}} } \right)^ }} \times e^ \times t_ }} }} \times N_^ + 0.5}} }} \times e^ }} }} - K_ \times A_}2j}}$$

Differential equation (3):

$$\frac}A_}j}} }}}t_ }} = K_ \times A_}1j}} + K_ \times A_}2j}} - K_ \times A_}j}}$$

where Ka1i, Ka2i, CLi, Vi, Ni, MTTi, and Ktri represent Ka1, Ka2, CL, V, N, MTT, and KTR in individual i, respectively. θ1, θ2, θ3, θ4, θ5, and θ6 are the population typical values for Ka1, Ka2, CL, V, N, and MTT, respectively. η1i,η2i, η3i, η4i, and η5i are the individual random effects in individual i on Ka1, Ka2, CL, V, and MTT, respectively. They follow a normal distribution with a mean of 0 and a variance of \(\omega_,}}^}\), \(\omega_,}}^}\), \(\omega_,}}^}\), \(\omega_,}}^}\), and \(\omega_,}}^}\), respectively. ADepot1, ADepot2, and ACent represent the amount of leuprolide in depot 1, depot 2, and central compartments, respectively. tad represents the time after dose. F2 represents the fraction of leuprolide absorbed via delayed absorption. The differential equations 1, 2, and 3 describe the change in the amount of leuprolide over time in depot 1, depot 2, and the central compartments, respectively. The first term of the differential equation (2), which describes the input from the transit compartments, was transformed to a logarithmic form to prevent numerical difficulties for a large N during the minimization of the model in NONMEM.

The residual unexplained variability was modeled using the following equation:

$$C_}ij}} = C_}ij}} \times \left( } \right)$$

$$\varepsilon_ \sim N\left( \sigma_^ } \right)$$

where Cobsij and Cpredij represent observed and predicted concentration for individual i at time j, respectively. ε1ij represents the proportional error for individual i at time j, respectively. ε1ij follows a normal distribution with means of 0 and variances of \(\sigma_,}}^}\).

In the covariates assessment, none of the covariates, including body size indices and age, was found to be significant. The relationship between body weight and age with pharmacokinetic parameters is shown in Fig. 2. The parameter estimates of the final model, along with the corresponding percentage relative standard errors, are summarized in Table 2. Basic goodness of fit to the observed concentration plots is shown in Fig. 3. Further assessment of the model using visual predictive checks showed close agreement of the 5th, 50th, and 95th quantiles of the prediction-corrected observations with the 95% confidence intervals of the corresponding quantiles of the prediction-corrected simulated data (see Fig. 4). This demonstrates the selected model’s ability to describe the central tendency and variability of the observed data.

Scatterplot of pharmacokinetic parameters (clearance [CL] and volume of distribution [V]) versus age and body weight. Circles are observations

Goodness-of-fit plots. Circles represent observations; the solid gray lines show the zero residual line (A and C) or the line of unity (B and D); the red dotted line is the locally weighted regression (LOESS) line. Conc. concentration

Prediction-corrected visual predictive check for the final model. Circles represent prediction-corrected observations, and shaded areas represent 95% confidence intervals of the 5th, 50th, and 95th percentiles of prediction-corrected simulated data. Black lines represent 5th (dashed), 50th (solid), and 95th (dashed) percentiles of the prediction-corrected observations. Some values may appear below the lower limit of quantitation/2 because of prediction correction