Population Pharmacodynamic Modeling of Epoetin Alfa in End-Stage Renal Disease Patients Receiving Maintenance Treatment Using Bayesian Approach.

The ability to control dosage regimens of erythropoiesis-stimulating agents (ESAs) to maintain a desired hemoglobin (HGB) target is still elusive. We utilized a Bayesian approach and informative priors to characterize HGB profiles, using simulated drug concentrations, in patients with end-stage renal disease receiving maintenance doses of epoetin alfa. We also demonstrated an adaptive Bayesian method, applied to individual patients, to improve the accuracy of HGB predictions over time. The results showed that sparse HGB data from daily clinical practice were characterized successfully. The adaptive Bayesian method effectively improved the accuracy of HGB predictions by updating the individual model with new data accounting for within-subject changes over time. The Bayesian approach presented leverages existing knowledge of the model parameters and has a potential utility in clinical practice to individualize dosage regimens of epoetin alfa and ESAs to achieve target HGB. Further studies are warranted to develop an application for practical use.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

☑ The current pharmacokinetic/pharmacodynamic (PKPD) models are not able to describe sparse hemoglobin (HGB) concentration profiles in patients with end‐stage renal disease (ESRD) receiving maintenance doses of epoetin alfa in order to support dose optimization to achieve desired HGB concentrations.

WHAT QUESTION DID THIS STUDY ADDRESS?

☑ This study aimed to characterize HGB response to the maintenance treatment of epoetin alfa in the ESRD patient population, evaluate covariate effects on the PD parameters, and assess potential use of an adaptive Bayesian method to improve the accuracy of HGB predictions.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

☑ Population PD modeling using Bayesian approach and informative priors can be used to describe sparse HGB data without epoetin alfa concentrations. The PD parameters can be updated with new data to improve prediction accuracy. An adaptive Bayesian method applied to individual data addresses between‐subject and within‐subject variability.

HOW MIGHT THIS CHANGE DRUG DISCOVERY, DEVELOPMENT, AND/OR THERAPEUTICS?

☑ This approach has a potential utility to develop a dosing algorithm of epoetin alfa for dose recommendation in individual patients.

Anemia is a common complication of end‐stage renal disease (ESRD) associated with fatigue, weakness, dyspnea, increased morbidity, and mortality. Recombinant human erythropoietin (EPO; epoetin alfa) is the first erythropoiesis‐stimulating agent (ESA) introduced late 1980s for the treatment of anemia in patients with chronic kidney disease. Since then, it has revolutionized anemia management, helping to prevent or minimize the use of blood transfusions and improve quality of life. , Despite its long history, the ability to control dosage regimen of epoetin alfa as well as other ESAs to maintain a target hemoglobin (HGB) concentration of 10–12 g/dL is still elusive. HGB response is often observed with large oscillations and overshoots out of the desired range. , Furthermore, using high doses of ESAs to achieve the target HGB has been found associated with increased risks of death and cardiovascular events. Therefore, individualized dosage regimen is recommended but sufficient guidelines are not available. ,

Extensive studies have been conducted in healthy volunteers and patients to characterize pharmacokinetics (PKs) and pharmacodynamics (PDs) of epoetin alfa. , , , , However, the potential utilities of these studies have not been realized in clinical practice to individualize dosage regimen. There are challenges in bringing established PK/PD models into clinical use. These models are developed using rich data, including baseline information from individual subjects. However, in practice, only sparse HGB observations are collected, often once or twice a month, and PK data include only dose amounts without drug concentrations. For patients on maintenance doses, baseline HGB measurements are not available. In addition, clinical changes within patients over time may influence the response to drug treatment. , In order to utilize these PK/PD models in practice, an approach that can integrate available knowledge about the models and new sparse clinical data is needed.

The Bayesian approach offers an effective way to leverage existing PK/PD models of epoetin alfa combining with new data to update the model parameters. , Available information about model parameters can serve as priors in the analysis of new data to describe HGB response and ultimately to optimize doses in individual patients. This approach is computationally intensive, but with the advancements in computing power and Markov chain Monte Carlo algorithms it is becoming more accessible.

Adaptive Bayesian approach has been studied for dose optimization. , The principle of this approach is that the prior distributions of the model parameters are updated sequentially over time with new data and their corresponding posterior distributions are used to perform simulation to obtain individualized dosage regimens. This method can be applied to anemia treatment in patients with ESRD given long‐term ESA therapy. Because changes within patients can affect their response to the therapy, , an adaptive Bayesian model updated periodically with new HGB data is expected to improve HGB predictions and assist in dose optimization.

In this study, we applied Bayesian approach to characterize HGB profiles, without drug concentrations, in patients with ESRD given maintenance doses of epoetin alfa. We also assessed covariate effects on the PD parameters. Furthermore, we demonstrated the use of an adaptive Bayesian model for individual subjects to improve prediction accuracy of HGB over time for a potential application in individualizing dosage regimens of epoetin alfa and ESAs.

METHODS

Study population

Three groups of patients with ESRD were used in this analysis. The first group, including 22 subjects from a retrospective study (672972‐1), was used for model fitting and internal validation. The second group, including 25 subjects from a prospective study (STUDY00000261), was used for external validation. All patients in this prospective group underwent informed consent prior to research participation. The third group including four subjects (from the retrospective study 672972‐1) was used to demonstrate the adaptive Bayesian method. Both studies were approved by the University at Buffalo Institutional Review Board and in accordance with the 1964 Helsinki Declaration and its later amendments. All subjects were enrolled at the outpatient hemodialysis center, Erie County Medical Center in Buffalo, NY. The patients were given maintenance doses of epoetin alfa and did not receive blood transfusions. The epoetin alfa dosing algorithm was standard among all subjects over the study period using an internal protocol. This protocol was based on package insert recommendations. Demographic and clinical data were obtained up to 1 year for the first and third groups (January 2014 to September 2015), and up to 14 weeks for the second group (July 2017 to May 2018).

Data sets

The end point of interest was the HGB observations. Plasma concentrations of epoetin alfa were not available in any patient, only dose information was recorded. Three data sets were derived such that there were at least three HGB observations and multiple epoetin alfa doses per subject. Data A was extracted from the first patient group (study 672972‐1) within the first 22 weeks of the collected data. It was divided into two parts, part 1 (first 12 weeks) was used as a training set for modeling and part 2 (last 10 weeks) was used for internal model validation. Data B was extracted from the second patient group (study STUDY00000261) over 14 weeks and used for external validation. Data C was derived from the third patient group (study 672972‐1) over 22 weeks and used for the adaptive Bayesian method. In general, in all data sets, HGB concentrations were obtained biweekly. Demographics and clinical characteristics are summarized in Table .

Table 1

Summary of the data sets used in the study

Variable, unit	Median (min–max) or number
Data A a , part 1, 22 subjects, over 12 weeks:
Demographic
Sex, female/male b	8/14
Age, years	59.0 (35.0–87.0)
Weight, kg	85.0 (48.0–149)
Height, inches	65.0 (51.0–76.0)
Clinical
Serum creatinine, mg/dL	8.40 (4.20–13.4)
Glucose, mg/dL	164 (81.0–246)
Albumin, g/dL	4.00 (3.40–4.60)
Transferrin saturation, %	28.1 (21.0–43.5)
Ferritin, ng/mL	881 (573–1166)
Diabetes, yes/no b	13/9
Hyperlipidemia, yes/no	13/9
Coronary artery disease, yes/no	7/15
Heart failure, yes/no	8/14
ERI c IU/kg/week/g/dL	11.6 (0.606–96.0)
Pharmacokinetic/pharmacodynamic
HGB concentrations, g/dL, n = 131 d	10.9 (7.8–12.4)
Epoetin alfa subcutaneous doses, IU, n = 342 e	11,000 (1,000–35,000)
Data A a , part 2, 22 subjects, over 10 weeks:
HGB concentrations, g/dL, n = 95	10.9 (7.8–12.9)
Epoetin alfa subcutaneous doses, IU, n = 270	11,000 (1,000–35,000)
Data B a , 25 subjects, over 14 weeks:
HGB concentrations, g/dL, n = 127	10.7 (8.5–14)
Epoetin alfa subcutaneous doses, IU, n = 215	4,000 (1,000–48,000)
Data C a , 4 subjects, over 20 weeks:
HGB concentrations, g/dL, n = 98	10.8 (8.6–12.3)
Epoetin alfa intravenous doses, IU, n = 163	10,000 (4,000–22,000)

ERI, erythropoietin resistance index; HGB, hemoglobin.

Data A and C from the retrospective study 672972‐1, data B from the prospective study STUDY00000261.

Binomial covariates coded as 1 for yes or female and 0 for no or male.

Calculated by weekly average epoetin alfa dose per body weight divided by average HGB concentration over the duration of data used for modeling, 12 weeks.

Number of HGB observations.

Number of epoetin alfa doses.

Structural PK/PD model

Because drug concentrations were not available, a validated population PK model of epoetin alfa in healthy subjects was utilized to simulate PK profiles of the patients. The simulated data were then used in an established PD model to describe HGB response to the epoetin alfa treatment. The resulting PK/PD model and parameters are presented in Figure . The PK simulation accounted for between‐subject variability (BSV) with the effect of sex, age, and body weight on k; body weight on F; sex on baseline endogenous EPO; age on fr; and age on V. The effect of baseline HGB on F was not applicable because data on this covariate were not available. The PD part included five transit compartments describing the red blood cell (RBC) lifespan. The amount of HGB in each compartment was the product of RBC count and mean corpuscular hemoglobin, which was assumed to be constant within a subject. The PD parameters to be estimated were T RBC, S max, SC 50, and initial HGB.

Figure 1

Schematic diagram of the pharmacokinetic/pharmacodynamic model of epoetin alfa given subcutaneously. A1, amount of epoetin alfa in the depot compartment; Ac, amount of epoetin alfa in the central compartment; Ap, amount of epoetin alfa in the peripheral compartment; C, concentration of epoetin alfa in the central compartment; CL, linear clearance; D 1, duration of zero‐order input into the depot compartment; D 2, duration of zero‐order input into the central compartment; F, absolute bioavailability; fr, bioavailability fraction; HGB1–HGB5 concentrations of hemoglobin in the five transit compartments, respectively; k a, first‐order absorption rate constant; K m, plasma concentration of epoetin alfa that gives 50% of the maximum elimination rate; Q, inter‐compartment clearance; SC 50, plasma concentration of epoetin alfa that gives 50% of the maximum production rate of HGB; Smax, maximum production rate of hemoglobin; T RBC, RBC lifespan; V c, volume of distribution of the central compartment; V max, maximum elimination rate of the saturable pathway; V p, volume of distribution of the peripheral compartment.

The model was described by differential Eqs. (1), (2), (3) for the PK and Eqs. 4 and 5 for the PD as follows:

Population PK parameters used for the PK simulation were obtained from Olsson‐Gisleskog

HGB1_IN to HGB5_IN were initial HGB concentrations in the five transit compartments, respectively. Because the patients were on maintenance doses, there were no baseline data for estimating baseline HGB (calculated secondarily using T RBC, S max, and SC 50 ) in the transit compartments. Instead, an additional parameter, HGBIN, was used to account for initial HGB in those compartments (i.e., HGB1_IN to HGB5_IN). For simplicity, it was assumed that:

The systemic concentration is the sum of HGB concentrations in the five compartments:

Population Bayesian PD model

HGB data were modeled using the Bayesian approach described previously for analysis of population PK data. The process included three stages:

Stage 1: Model for HGB concentrations

HGB concentrations were assumed to be log‐normally distributed with an additive error model: where y was the j th HGB observation of the i th subject at time t was the vector of individual PD parameters  = (S maxi, SC 50i, T RBCi, and HGB INi) for the i th subject, HGB( t) was the model‐predicted HGB for i th subject at time t, and ε was the residual error assumed to be independently and normally distributed with mean zero and variance σ 2.

Stage 2: Model for the BSV of the PD parameters

Individual subject PD parameters were modeled using a multivariate normal distribution on the logarithmic scale: where was the vector of the population PD parameters, Ω was the variance‐covariance matrix for the individual subject parameters, and N 4 was a 4‐dimensional multivariate normal distribution, ω was the vector of the SDs for the PD parameters representing BSV, Ρ was the correlation matrix of the individual parameters, and L (a 4‐by‐4 lower triangular matrix) was the Cholesky factor. The variance‐covariance matrix Ω was broken down so that ω and L were modeled instead of Ω to facilitate the convergence of the model parameters.

Stage 3: Model for the priors

Priors for σ, , and were assumed to be normally distributed on logarithmic scale, and L followed a Cholesky Lewandowski–Kurowicka–Joe (LKJ) correlation distribution : where, ln(μ) and b were the mean and SD of ln(σ), ln() and were the mean and variance‐covariance matrix of ln(), ln() and were the mean and variance‐covariance matrix of ln(), k was the parameter for the Cholesky LKJ correlation distribution denoted by lkj_corr_cholesky(k). Informative priors obtained from Wu et al. were used in this study. and were chosen as follows:

It should be noted that the values of 0.25 in and , used as the variances (or 0.5 as SDs) for the priors of the population PD parameters, created sample spaces much larger than using the relative standard errors from the reference study. As a result, the influence of the priors was reduced giving more weight to the HGB data in estimating the posterior distributions of the model parameters.

Parameter estimation

Markov chain Monte Carlo samples of the PD parameters were drawn from target distributions using No‐U‐Turn Sampler algorithm, an extension to Hamiltonian Monte Carlo. Model fitting was performed with four Markov chains run in parallel and 2,500 iterations per chain. For each chain, the first 500 iterations, warmup samples, were discarded and the last 2,000 iterations were kept for Bayesian inferences. Because the informative priors were used, 2,000 iterations per chain were sufficient to achieve stable posterior distributions. Convergence of the posterior distributions were assessed by Gelman‐Rubin statistic  < 1.05, trace plots showing random mixing horizontally of the four Markov chains. A successful run also showed no divergent transitions, no saturation of maximum tree depth, and Bayesian fraction of missing information > 0.2.

Covariate analysis

Information on covariate analysis is provided in Supporting Information file Doc .

Goodness of fit and evaluation of model performance

The goodness of fit of the model was evaluated using typical diagnostic plots, including observed vs. predicted plot, individual weighted residual plot, and individual and population visual predictive check (VPC) plots. The predictive performance was assessed using VPC plot, median prediction error (MPE), and median absolute prediction error (MAPE) for internal validation as well as external validation. These evaluations were conducted by simulations using 4,000 samples from the posterior distributions of the model parameters. Individual subject parameters were used for individual plots and internal validation. Population parameters were used for population plots and external validation. MPE and MAPE were calculated from prediction errors (PEs) of all HGB observations to be predicted in the corresponding situations using Eq. 8:where, and was the prediction error of the j th HGB concentration of the i th subject.

Adaptive Bayesian method for individual subjects

The posterior distributions of the population PD parameters obtained by fitting the base model were used as priors to fit the PD model to individual data (data C). At step 1, the individual model was fitted using HGB data in the first 12 weeks. Then, 4,000 samples of individual parameters drawn from the resulting posterior distributions were used to do simulation predicting HGB concentrations in the next 10 weeks. At step 2, HGB data in the first 5 weeks were excluded and the model was adaptively updated by refitting with the data from week 6 to week 17. Then, similar to step 1, simulation was performed to predict HGB concentrations of the last 5 weeks. The predictions of the last 5 weeks were compared between step 1 and step 2 using MPE, MAPE, and VPC plots generated as in the previous section.

Software

Statistical and graphical analyses were conducted using RStudio version 1.2.1335 (https://www.rstudio.com/) together with R version 3.5.2. The Bayesian PD model was implemented using Stan with PK/PD library Torsten 0.85 (Metrum Research Group LLC, Tariffville, CT; https://github.com/metrumresearchgroup/Torsten) via rstan version 2.18.2 package in RStudio.

RESULTS

The Base Bayesian PD model

The modeling process started by fitting the base PD model to data A, part 1. The model included BSV for all four parameters (i.e., HGB IN, T RBC, S max, and SC 50). The results showed that the model was able to capture HGB profiles of all the subjects (Figure ). Convergence of the parameters was achieved with all and good diagnostic trace plots (Supplementary Information Figure ). There were no divergent transitions or saturations of maximum tree depth, and the Bayesian fraction of missing information values for four chains were > 0.2, indicating the parameter space was explored sufficiently. The observed vs. predicted plot showed a good agreement between observed and model‐predicted HGB concentrations (Figure ). The posterior distributions of the parameters for the base model are summarized in Table . The model code and data structure are provided in Supplemental Material Code and Supplemental Material Data (Supporting Information files).

Figure 2

Prediction of new hemoglobin (HGB) concentrations 10 weeks in advance in subjects used for modeling. Empty circles are HGB concentrations used for modeling (data A, part 1), large, solid triangles are HGB concentrations (data A, part 2) to be predicted using the final pharmacodynamic model, lines are median predictions, shaded bands are 95% prediction intervals, and small, solid triangles are individual doses of epoetin alfa.

Figure 3

Diagnostic and model evaluation plots. Observed vs. median individual prediction (a) and median population prediction (b); solid lines are identity lines, circles are observed hemoglobin (HGB) concentrations, dashed lines are loess curves. Population visual predictive check plots for subjects used for modeling (c) and subjects used for external validation (d), two representative individual subjects for each plot; circles are observed HGB concentrations, triangles are individual epoetin alfa doses, lines are median predictions, bands are 95% prediction intervals.

Table 2

Posterior distributions of the pharmacodynamic parameters from the base model

Parameter, unit	Description	Median (95% CrI b )	Prior c (typical value)
θ HGBIN, g/dL	Typical initial HGB concentration	10.8 (10.3–11.4)	10.8
θ TRBC, days	Typical RBC lifespan	65 (41–114)	73.6
θ Smax, g/dL/hour	Typical maximum stimulatory effect of epoetin alfa	0.00762 (0.00443–0.0118)	0.00858
θ SC50, mIU/mL	Typical potency of epoetin alfa	3.71 (1.30–8.11)	5.35
ω HGBIN	Between subject variability of log(HGBIN)	0.105 (0.0759–0.151)	0.10
ω TRBC	Between subject variability of log(TRBC)	0.284 (0.133–0.573)	0.38
ω Smax	Between subject variability of log(Smax)	0.242 (0.109–0.515)	0.38
ω SC50	Between subject variability of log(SC50)	1.09 (0.409–2.14)	1.05
σ	Residual variability (log scale)	0.0416 (0.0358–0.0492)	0.05

CrI, credible interval; HGB, hemoglobin; RBC, red blood cell.

The base model is also the final model.

The 95% credible interval.

Typical values of the prior distributions of the model parameters.

Covariate analysis

Covariate effect on three parameters, T RBC, S max, and SC 50, was analyzed for demographic and clinical variables (Table ) except erythropoietin resistance index (ERI). However, no meaningful covariate effect was found. The inclusion of the covariates in the model did not meet the predetermined criteria (BSV reduced by 20%, 95% credible interval excludes zero, lower leave‐one‐out cross‐validation criteria value). A post hoc analysis using simple linear regression showed a significant association between ERI and median individual SC50 values with the slope = 0.117, P < 0.001, R 2 = 0.56 (Figure ).

Figure 4

Relationship between drug potency (SC 50) and erythropoietin resistance index (ERI). Slope = 0.117, P < 0.001, R 2 = 0.5595. Median SC50 values are obtained from the base pharmacodynamic model (also the final model).

Evaluation of the model performance

The performance of the final PD model, which was also the base model, was adequately validated. The individual prediction and population prediction VPC plots demonstrated that the model well‐described HGB profiles of the individual subjects (Figures , ). The prediction bands capturing the range of HGB concentrations indicated that BSVs and residual variability were well‐estimated. Internal validation results showed a good agreement between predictions and observed HGB 10 weeks in advance in the same subjects used for modeling (Figure ). The accuracy of the prediction was represented by MPE = −1.3% (interquartile range (IQR): −5.7 to 3.5%) and MAPE = 4.9% (IQR: 1.9–8.2%). External validation also shows a reasonable predictive performance of the model (Figure ). The HGB concentrations up to 14 weeks in a new group of subjects who were not used for the modeling process were predicted with MPE = −6.7% (IQR: −11.5% to 0.68%) and MAPE = 8.3% (IQR: 4.0–12.4%).

The PD model was fitted to individual data separately for four subjects (Figure ). The results showed that, at step 1 (Figure , A panels), the predictions of HGB concentrations in the last 5 weeks gave MPE = 11.5% (IQR: −7.6% to 16.6%) and MAPE = 13.5% (IQR: 10.8–17.8%); at step 2 (Figure , B panels), the predictions of the same HGB concentrations gave MPE = 3.4% (IQR: 1.7–6.3%) and MAPE = 4.6% (IQR: 3.3–6.7%). By fitting the model adaptively with new data, the HGB predictions were improved with better accuracy.

Figure 5

Adaptive Bayesian modeling of individual subjects. Data include four patients with end‐stage renal disease given intravenous doses of epoetin alfa (data C). At step 1 (A panels), the pharmacodynamics (PDs) model was fitted to the hemoglobin (HGB) data over 12 weeks (black circles), then HGB concentrations in the next 10 weeks were predicted (empty circles). At step 2 (B panels), the data in the first 5 weeks were excluded (grey squares and triangles) and the PD model was fitted to the data from week 6 to week 17 (black circles), then HGB concentrations in the last 5 weeks were predicted (empty circles). Black and grey triangles represent individual doses of epoetin alfa.

DISCUSSION

The Bayesian approach and previously published PK and PD models of epoetin alfa were utilized to characterize the HGB profiles in patients with ESRD given maintenance doses of epoetin alfa. The use of the PK model in healthy subjects is supported by previous reports , showing that PK of epoetin alfa in healthy subjects and patients with ESRD are similar. We chose the Bayesian approach because we were not successful in using traditional methods (maximum likelihood or extended least squares) to model these data. The reason was that HGB data (data A, part 1) collected from clinical practice were sparse with median of 5 observations per patient over 12 weeks (range: 3–12 observations) and baseline measurements were not available. The Bayesian approach allowed for incorporating informative prior distributions of the PD parameters from Wu et al. with similar patient population. In addition, the uncertainty of the parameters was accounted for by specifying the SDs of the priors. As described in the model for the priors, we specified larger SDs for the priors than the results from Wu et al. (6.7% vs. 2.3%, 10.5% vs. 3.6%, and 29.8% vs. 6.8% for T RBC, S max, and SC 50, respectively; and 21.1% for HBGIN). These assignments create larger parameter space but less weight for the priors and more weight for the data to inform the parameters. Increased SDs for the priors cause loss of information from previous studies and wider posterior distributions of the parameters, whereas narrower priors cause the reversed effect. This is in contrast with traditional methods where point estimates are used and uncertainties of the parameters are not conveniently taken into consideration in new analyses.

The PD model has an additional parameter, HGBIN, compared with the model by Wu et al. Because the data do not contain baseline HGB, this parameter is needed to represent the initial HGB at time zero. In the study by Wu et al., with baseline data, the concentration at time zero is the baseline that can be calculated secondarily using the other parameters (i.e., T RBC, S max, and SC 50). The use of HGB IN allows for the PD model to be fitted to the HGB data at any time during the treatment. This is particularly useful because we often do not know the baseline HGB and patients with ESRD are dependent on long‐term anemia treatments with within‐subject changes over time. ,

The population modeling results show that the PD model well captures all HGB profiles in individual patients. The individual median predictions well describe the curvature of the observed data and the individual prediction bands are narrow. These indicate that the sparse HGB data are informative for the PD parameters. Median RBC lifespan estimate of 65 days is comparable with reported results in similar patient populations, such as 64 days by Uehlinger et al., 62.8 days by Korell et al., and 58.1 days by Vos et al. Because other published PD models have different structures and parameterizations, the comparisons for S max, SC 50, and HGB IN are not reasonable. Compared with the results from Wu et al., the parameter estimates are in agreement. It should be noted that the study by Wu et al. did not use the Bayesian approach, instead it used the maximum likelihood estimation methods in NONMEM. The point estimates of the parameters from their study are well‐contained around the middle of the posterior distributions from our study. This is expected because the patient populations in both studies are similar. The study by Wu et al. was conducted using extensive HGB data from 118 subjects and was externally validated with data from 201 subjects. Therefore, the use of the PD parameters from this study as informative priors in our study is an effective way to leverage existing knowledge to analyze sparse HGB data from clinical practice.

The covariate analysis was conducted but no meaningful effect on the PD parameters was found. It is well‐known that iron deficiency (either absolute or functional) causes suboptimal response to EPO treatment in patients with ESRD. , Absolute iron deficiency can be addressed by iron supplementation. In fact, 8 of 22 subjects in data A were administered iron therapy. Levels of serum ferritin and transferrin saturation (Table ) are much greater than the references of 200 ng/mL and 20%, respectively, indicating that the patients have enough iron supply for erythropoiesis. This may be the reason that the covariate analysis did not identify ferritin or transferrin saturation as a meaningful covariate.

A post hoc analysis showed a strong positive association between SC 50 and ERI (Figure ). ERI is a complex metric representing the responsiveness to the treatment of anemia. Patients with higher ERIs are less responsive to the drug effect. In the PD model, SC 50 reflects the sensitivity of HGB response to the drug effect. Patients with lower SC 50 are more responsive to the therapy, whereas patients with increased SC 50 indicate weaker response or resistance. Therefore, the relationship between ERI and SC 50 can be considered as an additional test to see whether the model is appropriate to describe HGB data. Despite the strong association, ERI was not evaluated for covariate effect on SC 50 because it is calculated by the average weekly dose of epoetin alfa per body weight divided by the average HGB concentration over the study period. Therefore, using it as a covariate is not desired because dose information and HGB data are already used in the modeling process.

A practical utilization of the results from this study is the use of the parameter posterior distributions as priors to fit the PD model to individual data to improve HGB predictions. We demonstrated this application with four patients with ESRD given maintenance intravenous doses of epoetin alfa using an adaptive Bayesian method. This is especially relevant in patients with ESRD because they are on long‐term ESA treatment with within‐subject clinical changes over time, including blood loss and hemolysis, invasive procedures with blood loss, gastrointestinal bleeding, acute infection and inflammatory condition, iron therapy, etc. , The Bayesian model applied to individual data allows the informative priors of the population parameters to inform the estimation of the individual parameters. In addition, by adaptively updating the model parameters over time individually, clinical changes within subjects that may affect the PD parameters are accounted for and BSV is irrelevant. The duration of 5 weeks with new data to update the model corresponds to 2–4 HGB observations collected in clinical practice that might be informative for dose adjustment in patients. The results show that this approach is promising and feasible, which warrants further investigation with a larger number of subjects. When the ability to predict HGB concentrations in individual subjects is improved, simulation can be performed to select optimal dosage regimens to achieve desired target HGB.

To achieve 95% of the HGB observations to be predicted within the target range of 10–12 g/dL, the maximum PE should be around 4.6% (~ 0.5 g/dL/11 g/dL). The MAPE results from the internal validation using individual parameters are 4.9% (IQR: 1.9–8.2%) for predictions 10 weeks in advance with 95 HGB observations, and 4.5% (IQR: 1.9–7.8%) for predictions 5 weeks in advance with 45 observations. The MAPE from the adaptive Bayesian method for individual subjects is 4.6% (IQR: 3.3–6.7%) for predictions 5 weeks in advance with 16 observations. These results show that the approach used in this study is quite close to achieve the target and gives improved precision compared to the results from Wu et al. (parameters for the proportional plus additive error model are 6.1% and 0.5, respectively). However, we acknowledge that the number of subjects in this study is limited and further studies with larger sample size are needed to validate the approach.

We demonstrated that Bayesian approach with informative priors can be used effectively to characterize sparse HGB profiles in patients with ESRD receiving maintenance treatment of epoetin alfa. The model can be fitted to individual data adaptively over time to improve prediction accuracy of new HGB concentrations. Adaptive Bayesian modeling applied to individual data addresses within‐subject variability by updating the model with new data periodically. A dose adjustment algorithm, which is a critical component of the approach, needs to be developed in future studies for dose recommendation in individual patients.

Funding

This study was supported by the Bristol Myers Squibb fellowship to L.M.N. and a grant from the American Association of Colleges of Pharmacy to C.J.M. The study was performed independently of all funders.

Conflict of Interest

All authors declared no competing interests for this work.

Author Contributions

W.K., L.M.N., C.J.M., G.G.R., and M.P. wrote the manuscript. W.K., L.M.N., and C.J.M. designed the research. W.K., L.M.N., and C.J.M. performed the research. W.K. and L.M.N. analyzed the data.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.