Simultaneous Exposure-Response Modeling of ACR20, ACR50, and ACR70 Improvement Scores in Rheumatoid Arthritis Patients Treated With Certolizumab Pegol.

The Markovian approach has been proposed to model American College of Rheumatology's (ACR) response (ACR20, ACR50, or ACR70) reported in rheumatoid arthritis clinical trials to account for the dependency of the scores over time. However, dichotomizing the composite ACR assessment discards much information. Here, we propose a new approach for modeling together the three thresholds: a continuous-time Markov exposure-response model was developed, based on data from five placebo-controlled certolizumab pegol clinical trials. This approach allows adequate prediction of individual ACR20/50/70 time-response, even for non-periodic observations. An exposure-response was established over a large range of licensed and unlicensed doses including phase II dose-ranging data. Simulations from the model (50-400 mg every other week) illustrated the range and sustainability of response (ACR20: 56-68%, ACR50: 27-42%, ACR70: 11-22% at week 24) with maximum clinical effect achieved at the recommended maintenance dose of 200 mg every other week.

The American College of Rheumatology's (ACR) core dataset[1] has caused a major improvement in the standardization of clinical trials for rheumatoid arthritis (RA). Response criteria have been defined as 20%, 50%, and 70% improvement from baseline and named ACR20, ACR50, and ACR70. ACR20[2] has been accepted as the efficacy benchmark in RA clinical trials,[3] showing greater discriminant capacity over ACR50 and ACR70 to distinguish active treatment from placebo control.[4] However, 20% improvement does not represent optimal clinical progress and is not a meaningful clinical response for rheumatologists.[5] In clinical trials conducted since 1997, similar discriminant capacities were demonstrated for ACR20 and ACR50 to distinguish active versus control treatments, which seems to attest the superiority of newer therapeutic approaches and drugs, achieving higher level of ACR response.[6] ACR50 and ACR70 are more desirable targets for patients and provide useful information in addition to ACR20.[6,7] Dichotomizing the composite ACR assessment into binary variables such as ACR20, ACR50, and ACR70 discards information.[8,9] Classification as “ACR20 responder” ignores whether the subject actually reached 50 or 70% improvement over baseline. Thus, simultaneous modeling of the three available endpoints ACR20, ACR50, and ACR70 is a more effective use of the data.[10] ACR-N[11,12,13] or hybrid-ACR[14] analysis would represent a further step toward use of a continuous measure based on ACR assessment.

Successive ACR assessments collected repeatedly over time for each individual are not independent.[15] Not accounting for this property in the implementation of ordered logistic regression, the standard method to model ordered categorical data[16] is expected to lead to over-prediction of the number of transitions between the different grades.[17] The inclusion of Markov elements in the logistic regression[18] lead to better characterization of the transitions between response and non-response for the binary score ACR20.[15] Nevertheless, the complexity of the model increases dramatically with the number of categories and constant influence of preceding score(s) is implicitly assumed, whatever the observation schedule. The latent variable (LV) approach, that relates the observed scores to an underlying continuous measure of the disease activity,[19] was proposed to model the three ACR scores simultaneously.[20,21] This approach essentially focuses on describing the average probability of the population of subjects and was shown to over-predict the number of transitions between scores when applied to stand-alone ACR20 modeling in its simple implementation.[15] It was recently extended to accommodate extra-correlation between longitudinal dichotomous data[22] but remains complex to implement and is not integrated in standard software packages. The continuous-time Markov model was first proposed in the PKPD literature to model the tablet position in the gastro-intestinal tract.[23] It was applied to characterize various types of ordered categorical data, such as pain in animal models[24] and side effects in clinical trials: grades of proteinuria in cancer patients treated with anti-angiogenic drug[17] or extrapyramidal side effects in schizophrenic patients treated with antipsychotic drugs.[25] This approach addresses the dependence between successive observations with the influence of the previous score decreasing with increasing time between observations. It is more sparing than inclusion of Markov elements in proportional odds models when more than two categories are considered and allows for accurate simulations for various observation schedules.

In this article, we present the development of a continuous-time type mixed-effect Markov model characterizing the ACR20, ACR50, and ACR70 responses as a function of certolizumab pegol exposure from five phase II and phase III clinical trials including both licensed and unlicensed doses. Certolizumab pegol is PEGylated Fc free anti-TNF (Fc = crystallizable fragment) licensed for treatment of moderate to severely active RA in several countries, including United States and Europe, at the dose of 400 mg initially and at weeks 2 and 4, followed by 200 mg every 2 week (Q2W); for maintenance dosing, 400 mg every 4 weeks (Q4W) can also be considered.

Results

Data

The patient population and study design of the five double-blind placebo-controlled clinical trials[26,27,28,29] have been described previously.[15] Data from 2,380 subjects with RA were analyzed (24,519 observations); 1,747 subjects received certolizumab pegol at doses ranging from 50 to 800 mg, and 633 subjects received placebo. Treatment was administered subcutaneously Q2W or Q4W for 12 weeks up to 1 year, as monotherapy or in combination with methotrexate (patients were on stable dose at inclusion), with a loading dose for the two pivotal studies (). The predicted plasma concentrations at the time of ACR assessment (Cp) ranged from 0.4 to 174 μg/ml. Typical steady-state average concentration for the 200 mg Q2W and 400 mg Q4W label treatments was 23.2 μg/ml.

Structural ACR model

At each visit, subjects were classified as ACR20 non-responder (ACR-NR, score = 0), ACR20 but not ACR50 responder (ACR20-50, score = 1), ACR50 but not ACR70 responder (ACR50-70, score = 2), or ACR70 responder (score = 3). The probabilities of the scores, modeled as compartment amounts, were defined by four ordinary differential equations, the sum of the probabilities of the four responder states remaining equal to 1 at any time point (). At the time of observation, the actual score is known, i.e., the probability of observing this score is equal to 1. The probabilities of the other responder states and the integration time were reset to 0 (see example dataset in Supplementary Material). Between two observations, the probability “amount” distributes with time between the different states, with a rate determined by the transition parameters. At each time, the total probability “amount” in the different states/compartments is 1 as the system is closed ().

with K = TVK·exp(η) and TVK = f (Kp,0, gp(t), hp(qij))

where Pr(0) to Pr(3) are the probability of the scores from ACR20-NR to ACR70, the K are the transfer rate parameters between the compartments where the first digit of the rate parameter represents the current score and the second digit the preceding score (e.g., K21 can be read as the transfer rate to score 2 given current score of 1). K reflects the baseline value of the pth transfer rate parameter at the time of the first dose of study medication, g(t) and h(q) the time and drug effects on the transfer rate, t the continuous time and q the exposure measurement for the jth observation of ith subject. The predicted concentration at the time of observation, Cp, was selected as marker for drug exposure based on the objective function value (OFV).

The η, representing the random effect around parameter K for subject i, was assumed to have a normal probability distribution of mean zero and variance ω2; its empirical Bayes predictions allowed for subject-specific predictions of the corresponding transfer rate.

Final model parameters are summarized in . Whenever possible, the model was reduced by combining effect parameters on different transfer rates. All patients were considered as ACR20-NR at baseline (time = 0) and the score could improve naturally after this time point. Later in the study, scores tended to worsen with time in placebo treated patients: transfer rate to lower scores were greater than transfer rate to higher scores, i.e., an ACR20 responder subject receiving placebo was more likely to revert back to ACR20-NR than to become ACR50 responder.

A positive exposure–response relationship was identified. Emax-models with positive Emax value (K10) and linear functions with positive slope (K21 and K32, common slope) described the effect of the drug concentration on upward transfer rate parameters. Emax-models with negative estimates of Emax, described the drug effect on downward rate parameters (K01, K12, K23, common Emax). K were restricted to be positive by setting the lower boundaries of the initial parameter estimates.

where slt and sld represent the slopes of the linear relationship to continuous time and concentration, respectively, Emax is the maximum change in the pth transfer constant value as a function of exposure measure and EC50 the concentration to reach half of the maximum effect related to Cp, common to the relevant transfer rate parameters.

Incorporation of a full variance-covariance matrix characterizing the correlations between the transition rate parameters markedly improved the fit of the model to the data (decrease in OFV of 338 units) and the simulation performance of the model. Transfer rate parameters to higher ACR were highly correlated (K10, K21: 62%, K10, K32: 60%, K21, K32: 82%). Transfer rate parameters to lower ACR were less correlated, except 79% for K01, K12. However, inclusion of these 15 additional parameters dramatically increased the run time, and resulted in failure of the covariance step.

Upward transfer rates were decreased (lower ACR response) with increased baseline swollen joint count and concomitant corticosteroid treatment. Downward transfer rates increased with increased age, resulting in decreased ACR response for elderly subjects. Upward transfer rates were higher for subjects originated from Central and Southern America than for North Americans and Western, Eastern, and Northern Europeans.

The concomitant use of methotrexate was not found to affect the probability of ACR response. This is not surprising, as patients had to be on a stable methotrexate dosage upon enrollment in the relevant studies. The ACR assesses improvement from baseline, thus the methotrexate effect should already be included in the baseline assessment.

Dropout model

The subjects' dropout status (yes/no) was assessed at each visit. The dropout probability was modeled using a logistic model. No random effect around the dropout probability was estimated because subjects could only drop out once during the trial. The increase in dropout with time from the start of the trial was described by a Hill function on the logit scale where the time to reach half of the maximum dropout was estimated to approximately 3 weeks. In addition, ACR70 and ACR50 responders had the same, lower probability of dropping out than ACR20 responders, who had lower probability of dropping out than ACR20 non-responders:

where φ reflects the baseline (time = 0) logit score of the dropout, Emax T the maximum change in the logit related to time, T50 the time to reach half of this maximum change. The parameters θACR20, θACR50, and θACR70 reflects the effects of previous ACR outcome (ACR20 = 1 if the patient was an ACR20-50 and 0 otherwise, similarly for ACR50-70 and ACR70).

Model evaluation

The model simulated well the observed ACR levels () and the observed proportions of transitions between the ACR levels (, drug data in panel b and placebo data in panel c) for the whole dataset. For a few of the schedules, e.g., 200 mg Q2W, some small misfit is apparent when the ACR clinical scores by treatment were predicted (). Inter-study or inter-treatment variability was not included, that may impact the fit of individual studies.

Simulations

Simulations were performed to illustrate the exposure-response relationship. Typical PK exposure was derived from the population PK model (71.5 kg subject who did not develop anti-drug antibody) and the ACR scores simulated from the continuous Markov model. For example, at week 12, the rate of response ranged from 53.7 to 67.3% for ACR20, from 22.4 to 35.0% for ACR50 and from 6.5 to 14.4% for doses ranging from 50 mg Q2W to 400 mg Q2W (). This range of doses extending outside the label dosages was chosen for illustration purpose, as a plateau in the response is rapidly reached and the dose of 400 mg Q2W is not markedly better in term of efficacy than the 200 mg Q2W dose. Subjects that dropped out were considered as non-responders in the computation of these statistics, corresponding to what was done in the traditional statistical analysis of the performed clinical trials (worst case scenario). At week 24, the ACR20 outcome from the two pivotal trials RAPID 1 (ref. 28) and RAPID 2 (ref. 29) were somewhat lower than simulated from the model (for 200 mg Q2W: 58.8 and 57.3% observed ACR20 response vs. 65.2% simulated, 37.1 and 37.7% observed ACR50 response vs. 37.7% simulated, 21.4 and 15.9% ACR70 response vs. 16.7% simulated). An explanation could be that compared to the simulations that were based on a typical subject's PK, these studies included subjects with lower exposure due to higher body weight and/or formation of anti-drug antibodies.

To exemplify the potential use of the model for clinical purpose, some other outcomes were derived from the simulations (). At long term, the probability of responder subjects at week 12 still being responder at week 48 increased with the dose: 72.1–77.3% for ACR20, 71.4–78.6% for ACR50 and 56.9–68.8% for ACR70 for treatment from 50 to 400 mg Q2W. At week 24, the likelihood of remaining responder was high, between 77 and 94%, for the three ACR scores whatever the dose, characterizing a sustained response. On the other hand, non-responders at W4 had between 41 and 48% chance to achieve ACR20 at week 12, 12–17% to achieve ACR50 and 2–5% to achieve ACR70 depending on the dose.

Discussion

This pooled analysis established the relationships between exposure to the anti-TNFα drug certolizumab pegol and the ACR20, ACR50, and ACR70 responses in five phase II and III clinical trials, including both licensed and unlicensed doses. The model development was based on the database that was previously used to develop stand-alone models for ACR20 (ref. 15) and ACR50 (ref. 10) responses; however, this new approach has the advantage of gaining more information from the same amount of data, as various levels of response are characterized simultaneously. The model adequately predicts both the proportion of response for the three ACR thresholds and the transition between the four statuses of response.

The exposure–response relationship was established over a large range of doses, including doses from the dose-ranging phase II study that are outside of the range of recommended doses for maintenance certolizumab pegol treatment in the RA indication, 200 mg Q2W and 400 mg Q4W. A plateau in the probability of becoming a responder and achieving a higher ACR response is reached, with the 200 mg and 400 mg Q2W doses providing similar level of response. The dropout is also affected, as subjects with higher level of ACR response are less likely to drop out. The model characterized a sustained response to treatment, a factor that is recognized as key in evaluating the treatment response and disease status in patients with RA,[7] the aim being to prevent joint damage and achieve sustained low disease activity or remission.[30] As such, the new EULAR/ACR recommendations request to include time to onset of response and sustainability of response when evaluating treatment response in clinical trials.[31,32]

Depending on the audience, each of the three ACR levels may provide the most relevant information. Regulatory authorities may be interested in the ACR20 and ACR50 responses to discriminate between drugs, rheumatologists may consider higher thresholds as a more relevant target to be achieved in clinical trials and might rather consider the time to reach ACR50, while patients will appreciate sustained ACR70 response to treatment. The model can be utilized to provide clinically interesting information through population simulations, with a large population, or through trial simulations, with a limited number of subjects and a specific design, according to the question of interest. For example, this could be the prediction of the long term response or the evaluation of the chance for an ACR20 non-responder at an early time point to achieve an ACR20, ACR50, or ACR70 response at a later time.

Accounting for the dropout and simultaneous modeling of the three ACR thresholds enables accurate trial simulations with specific study design. For example, in two of the studies used for the present analysis, a study-design dropout was defined based on the ACR20 response at weeks 12 and 14, which was an issue for developing a stand-alone ACR50 model.[10] It also allows for applying different post-processing rules to compute the response rate, e.g., to consider dropped out subjects as non-responders.

The dependency between the successive ACR assessments was accounted for by mimicking a compartment model where the probability flows between the compartments are estimated. Any transitions between the four levels of ACR response are allowed, but the model is parameterized such that only transitions between neighboring states are estimated, making it more parsimonious than a proportional odds model with Markov elements. For nonconsecutive observed scores (e.g., ACR20 to ACR70), it is assumed that the intermediate unobserved score (ACR50 in this example) must have been reached at an intermediate time. The influence of previous observations decreases with increasing time between observations, i.e., if observations are made with large time interval, the influence of the previous observation will be small compared if two observations are made close in time. The half-life for transfer between the compartments were found to be dependents on drug concentration and time in the study. In addition, the inter-individual variability allowed for different patients to have different transfer rates.

The latent variable (LV) approach has also been proposed to model the ACR20, ACR50 and ACR70 scores simultaneously.[20] The concept was, however, different, aiming at describing the average probability of ACR20/50/70 response in the population rather than an accurate description of an individual's evolution of the response with time. It had the appeal to be parsimonious while flexibility can be added by fitting separate intercepts for the different thresholds. The continuous Markov model enables larger flexibility by allowing different covariate and drug effects for each transition rate, but this increases further the complexity of the model. We tried, however, to keep the model as simple as possible by using the same effects on the different transfer rates whenever possible.

Similar to models for ordered categorical data and the ACR20 Markov model, the parameters of the continuous Markov model may be difficult to interpret directly, and their clinical relevance is most easily interpreted through simulations from the model. The latent variable model may be considered as a step further for a more mechanistic interpretation, with parameters such as Emax that can be interpreted directly. Run time may be another limitation of the continuous Markov model—coded as differential equations—particularly with large datasets and increased number of categories.

In conclusion, an exposure–response model that simultaneously characterized the temporal course of the ACR20, ACR50, and ACR70 assessments in patients suffering from RA enrolled in phase II and phase III clinical trials, treated with certolizumab pegol and placebo, was successfully developed. The model predicted that increased certolizumab pegol exposure resulted in an increased probability of attaining higher level of ACR response, with the majority of the clinical effect being attained at the exposure achieved with the 200 mg Q2W label maintenance dose. This approach is applicable to model any ordered categorical score where consecutive observations are dependent.

Methods

The present analysis was performed retrospectively on pooled data from five clinical trials, not following the protocol analysis and hence the results cannot be directly compared with the published results of the respective studies.

Clinical score: ACR20, ACR50, and ACR70 responses are defined as a decrease of 20, 50, or 70%, respectively, in both tender and swollen joint counts, as well as the same magnitude of improvement in at least three of the five other measures of the ACR core dataset: the patient's and physician's global assessments of disease activity, the patient's assessment of pain and physical function, and an acute phase reactant (C-reactive protein level or erythrocyte sedimentation rate). The C-reactive protein level was used to compute the ACR scores in the present analysis.

Study population and design: Subjects with active RA diagnosed at adult age since at least 6 months and inadequate response to methotrexate monotherapy were enrolled in the five clinical trials (). Study protocols were approved by the medical ethics committees or institutional review boards of the participating centers, and all subjects provided written informed consent prior to initiation of study procedures.[33]

In studies 4 and 5, subjects assessed ACR20 non-responders at both week 12 and week 14 withdrew at week 16 with the option of entering an open-label extension study (data not included in the analysis). This design feature was reproduced in the model evaluation.

Covariate data: Covariates tested during model development are summarized in .

Certolizumab pegol exposure measures: Certolizumab pegol plasma concentrations were determined using a validated enzyme-linked immunosorbent assay (ELISA) method. Individual exposure measurements (individual certolizumab pegol concentration at the time of clinical observation Cp and average certolizumab pegol concentration over the last dosing interval Cavg) were derived from a previously developed one-compartment population model using the empirical Bayes post hoc estimates from NONMEM. These exposure measures were set to 0 for placebo data.

Software and data analysis

The model was developed in NONMEM 7.1.2 (ref. 34) (ICON Development Solution, Hanover, MD) and finalized in NONMEM 7.3.0. launched with PsN (Uppsala University, Uppala, Sweden).[35] The Laplacian method was used to approximate the marginal likelihood; it uses a second-order expansion around the empirical Bayes predictions of the inter-individual random effects. R version 2.11 and 2.15 (ref. 36) (R Foundation for Statistical Computing, Vienna, Austria) were used for data management and graphical outputs.

Model development

Model structure: See Results section.

Model building and selection: During the model building, efforts were made to keep the model as simple and parsimonious as possible. Modeling was initiated based on placebo data. Various functions of increasing complexity, such as linear, Emax, Hill, second order and quadratic polynomial were tested to describe the influence of time since start of the study on the transfer rates between the four scores and on the dropout probability. After addition of data from certolizumab pegol-treated subjects, drug-effect models were evaluated. The effect of drug on the transition rates were modeled as functions of exposure, using a treatment flag, the dose, Cp or Cavg as markers for certolizumab pegol exposure. The exposure relationships for each of the transfer rate were explored using linear, Emax and Hill equations. Effects were tested separately on each transition rate and the parameters were shared between rate constants if found to be similar and this did not lead to a significant increase in OFV. Equations retained in the final model are presented in the results section.

Inter-individual variability was assessed for the six transfer rate parameters.

Covariate model: Covariates were included in the model based on their statistical significance. The potential covariate effects were added on all transition rate constants (see example in the NONMEM control file). Continuous covariates were introduced as

where θcov is the estimated covariate effect, CONT the continuous covariate value and CONTmed the median of the covariate in the dataset.

Co-medications were evaluated as binary covariates (use/no-use). Categorical covariates were introduced as

where the θcov,CATx are the estimated covariate effects, CATx = 1 for the corresponding covariate, CATx = 0 otherwise.

A manual stepwise approach was used, beginning with a univariate testing (α = 0.05). Significant covariates were included in both forward inclusion (α = 0.05) and backward elimination (α = 0.001) steps, which continued until only significant covariates remained in the model.

Model evaluation: Model selection was based on changes in the NONMEM OFV evaluated using a likelihood ratio test (nested models, P = 0.05) and simulations from the models (visual predictive checks) to evaluate the model's capacity to capture the observed frequency of scores and transitions between the scores over time. For that purpose, 250 datasets were simulated from the model of interest. Additional rows at planned visits were added in case of early dropout in the original dataset in order to reproduce the theoretical observation scheme and allow subjects to drop out according to model simulations. The frequency of ACR outcomes and transitions between the four ACR scores were plotted as a function of time and compared with the observed data. Residuals were not available for evaluating the goodness of fit because the model estimated the probability of scores and not a prediction of the scores themselves.

A bootstrap with 100 samples, stratified on study and treatment, was performed to evaluate the precision of estimation of the parameter estimates. For each model parameter, the median and the standard deviation of the parameters' distributions from the bootstrap runs were computed and reported.

Simulations: For performing simulations, a probability amount of 1 was allocated in the compartment corresponding to the simulated score (subjects are set to non-responder at baseline, see model file and an example simulation dataset in Supplementary Material).

Author contributions

B.D.L., M.O.K., and L.E.F. wrote the manuscript, designed the research, performed the research, and analyzed the data.

Conflict of interest/Disclosure

The clinical trials on which the manuscript is based were sponsored by UCB Pharma. Brigitte D. Lacroix is full-time employee of UCB Pharma. Mats O. Karlsson has received research funding from UCB Pharma, but not related to this work. L.E.F. declared no conflict of interest. Deputy Editor in Chief L.E.F. was not involved in the review or decision process for this paper.

Study Highlights

Table 1

Study design and available data

Table 2

Final model parameter estimates

Table 3

Covariates