The sequence of disease-modifying anti-rheumatic drugs: pathways to and predictors of tocilizumab monotherapy

Background There are numerous non-biologic and biologic disease-modifying anti-rheumatic drugs (bDMARDs) for rheumatoid arthritis (RA). Typical sequences of bDMARDs are not clear. Future treatment policies and trials should be informed by quantitative estimates of current treatment practice. Methods We used data from Corrona, a large real-world RA registry, to develop a method for quantifying sequential patterns in treatment with bDMARDs. As a proof of concept, we study patients who eventually use tocilizumab monotherapy (TCZm), an IL-6 antagonist with similar benefits used as monotherapy or in combination. Patients starting a bDMARD were included and were followed using a discrete-state Markov model, observing changes in treatments every 6 months and determining whether they used TCZm. A supervised machine learning algorithm was then employed to determine longitudinal patient factors associated with TCZm use. Results 7300 patients starting a bDMARD were followed for up to 5 years. Their median age was 58 years, 78% were female, median disease duration was 5 years, and 57% were seropositive. During follow-up, 287 (3.9%) reported use of TCZm with median time until use of 25.6 (11.5, 56.0) months. Eighty-two percent of TCZm use began within 3 years of starting any bDMARD. Ninety-three percent of TCZm users switched from TCZ combination, a TNF inhibitor, or another bDMARD. Very few patients are given TCZm as their first DMARD (0.6%). Variables associated with the use of TCZm included prior use of TCZ combination therapy, older age, longer disease duration, seronegative, higher disease activity, and no prior use of a TNF inhibitor. Conclusions Improved understanding of treatment sequences in RA may help personalize care. These methods may help optimize treatment decisions using large-scale real-world data.


Background
Rheumatoid arthritis (RA) is a chronic systemic autoimmune inflammatory condition that involves the upregulation of various cytokines and auto-antibodies over time [1]. While the precise "cause" of RA is unknown, there are strong genetic alleles and several key environmental exposures that strongly associate with the risk of RA [2]. The immunobiology of RA has been well described, and many treatments developed that inhibit key aspects of the inflammatory cascade [1]. The majority of patients with RA find treatments that reduce disease activity and pain, but not all patients are so lucky and the path to successful treatment may be long and littered with disease-modifying anti-rheumatic drugs (DMARDs) that do not work [3]. Thus, it is surprising that relatively few studies have focused on sequences of RA treatments [4].
Recommendations regarding RA treatment suggest starting with conventional synthetic DMARDs (csDMARDs) [5,6], but then there is little evidence about which bDMARD should be tried next. The evidence grows even scanter as one considers second-or third-line bDMARDs [7]. Sequential decisions may be determined by previous treatment response, comorbidities, patient and provider preferences, and insurance. Optimizing sequential decisions through clinical trials is difficult since the number of possible sequences is large. A pragmatic approach is to first identify and quantify patterns in current treatment practice and then evaluate these in comparative trials.
In this work, we studied sequences of RA treatments using discrete-state Markov models for high-level statistics and logistic regression for determinants of treatment choice. Markov models describe the probabilities of moving from one state to another in a discrete dynamical system using a state transition matrix [8]. These have been widely used in econometric modeling, costeffectiveness modeling, and in some clinical areas [9][10][11], but they have played little role in personalizing treatment strategies in RA. We use logistic models to discriminate between patients on the same treatment who are moved to different drugs in the next stage.
In this set of analyses, we examined pathways to tocilizumab monotherapy (TCZm). Monotherapy with this agent has proven efficacy over other bDMARDs used without a csDMARD [12][13][14][15]. Since many patients find using monotherapy preferable, we focused on TCZm and examined the sequences to and predictors of its use in a large real-world dataset from the US, Corrona. While this set of analyses may have specific practical relevance for TCZm treatment, it should be considered as an exploration of methods that can be considered when evaluating sequences of treatment for a chronic disease like RA. To illustrate this point, for comparison, we apply the same methodology to briefly characterize common sequences leading to monotherapy with TNF inhibitors. We did not assess the optimization of treatment strategies but rather how to analyze sequential data in ways that can guide experiment design for such optimization.

Study design and population
We examined transitions of RA treatments among a cohort of patients initiating a bDMARD. All patients were enrolled in the Corrona RA registry between 2002 and 2019, a large real-world data source in the USA [16]. Patients were followed longitudinally for up to 5 years to determine transitions between bDMARDs, specifically examining if they started TCZm, subcutaneous or intravenous. The visit prior to the first report of a bDMARD initiation was considered the baseline, and the visits during the ensuing 5 years were the follow-up period. The probabilities of transitioning from one type of bDMARD to another were calculated at 6-month intervals during follow-up. Predictors of TCZm and TCZ combination therapy (TCZc) were assessed using regression models and machine learning algorithms.

Outcome of interest
The primary outcome of interest was the initiation of TCZm and the sequence of RA treatments leading to TCZm. This was defined through the use of the case report forms that are collected every 4-6 months by the rheumatology clinical site in Corrona. Monotherapy was defined as no other DMARD treatment marked on the case report form at the time of TCZ use; we did not include sarilumab (another IL-6 antagonist) in these analyses.
Other RA treatments considered in the sequence analyses and the modeling included TCZc (TCZ with a csDMARD), TNF inhibitors without TCZ +/− csDMARD, non-TNF non-TCZ bDMARDs +/− csDMARDs, and csDMARDs without a bDMARD. csDMARDs were included as some patients transition back to a csDMARD after trying a bDMARD. JAK inhibitors (here Xeljanz) were counted as bDMARDs since they are used sequentially after synthetic DMARDs and target-specific immune pathways. We also included no DMARDs as a possible RA treatment state.

Covariates
Potential predictors of TCZm use were assessed at the visit prior to starting any bDMARD (baseline) and also during follow-up. We considered patient's age and sex at the baseline visit. We also considered a number of RA characteristics at baseline: disease duration, serologic status, erosion status, clinical disease activity index (CDAI) [17], health assessment questionnaire (HAQ) [18], and prior DMARDs. CDAI, HAQ, and prior DMARDs were updated during follow-up. In addition, we examined comorbidities at baseline. These were originally considered as individual comorbidities, and then, they were collapsed into a simple count variable, giving each condition a similar weight. These comorbidities are listed in Table 1.
At baseline, some patients had missing serologic status. We used all available information on rheumatoid factor and anti-CCP status up through the baseline visit.
However, serologic status remained missing on 8.9% and these values were categorized as missing in the analyses. The HAQ was missing at baseline for some patients; we used the non-missing value immediately prior to the baseline visit for such patients.

Statistical analyses
We described the baseline characteristics of patients included in the study cohort. Median values or numbers and percentages were assessed. The Markov transition Miscellaneous includes drug toxicity, fracture, and other less common comorbidities matrix for DMARD status was assembled using 6-month intervals over the 5 years of follow-up. We examined the transition probabilities for each bDMARD category. After starting a bDMARD, some patients had 6-month intervals with no DMARD or only csDMARDs noted; we also assessed the probability of no DMARD or only a csDMARD. Based on these probabilities, we assessed the most common sequences of medications used prior to TCZm. To further characterize patient characteristics associated with TCZm use, we constructed regression models for discriminating between patients who end up on different treatments. Since we noticed that most patients who used TCZm had used TCZc, we examined models that used both TCZ states as dependent variables. Variables were included based on forward selection: univariable logistic models considered one baseline variable at a time; variables with p < 0.2 were advanced to multivariable-adjusted regression. We built the baselineonly logistic regression model with selected baseline factors. To improve the prediction of TCZ use, we included variables characterized during follow-up using a linear logistic model with repeated measurements (PROC GLMMIX, SAS V9.0). We measured the within-person correlation with robust covariance estimates. Unlike a typical epidemiologic model focusing on baseline predictors, models with updated variables allow one to consider sequential changes in patient characteristics that may influence treatment decisions.

Results
We identified 7300 patients with RA in the Corrona registry who initiated a bDMARD after entry. Characteristics of these patients are described in Table 1. The median age at study baseline was 57.3 years, 77.6% were female, the median disease duration was 5 years, and 56.8% were seropositive. At baseline, 82.1% reported the use of another csDMARD (e.g., methotrexate) and 33.1% reported the use of glucocorticoids.
Of the patients starting bDMARDs during Corrona follow-up, 676 (9.3%) initiated TCZ at some point and 287 (3.9%) reported use of TCZm. The median time until any use of TCZ was 17.5 months and 25.6 (IQR 11.5, 56.0) months until TCZm. The median number of bDMARDs between initiation and TCZ was 1 (range 0-6); for TCZm the median was also 1 (range 1-5). Eighty-two percent of TCZm use began within 3 years of starting any bDMARD; 93% of TCZm users switched from TCZ combination, a TNF inhibitor, or another bDMARD. Common sequences of DMARDs until TCZm use, going back at most 4 treatments, are illustrated in Fig. 1. Each line represents a different trajectory of treatments (color) and a different number of instances (thickness). Only trajectories ending with TCZm are included. Common sequences from the initiation of a bDMARD through TCZm use include TNFi to TCZ combination to TCZm (5.9%), TNFi to another bDMARD to TCZm (9.8%), and TNFi to TCZm (13.9%). Very few patients are given TCZm as their first DMARD (0.6%) or their first bDMARD (0.9%). The full 5 years of transition probabilities after initiation of a bDMARD are described in Table 2. For comparison with Fig. 1, we applied the same methodology to illustrate common DMARD sequences leading to TNF inhibitor monotherapy (see Fig. 2). We see that such sequences in current practice are vastly dominated by transitions from TNFi combination therapy with nonbiologic DMARDs to TNFi monotherapy (72%). Sequences in which bDMARDs precede TNFi monotherapy are comparatively rare in this sample (9.6%).
To better understand the factors associated with the use of TCZ and specifically TCZm, we conducted two different sets of regression models. First, we examined baseline variables only as predictors of any TCZ use and of TCZm. Similar variables were found to be associated with both outcomes (see Tables 3 and 4). These include age, serologic status, and use of csDMARDs at baseline.
In a separate set of regression models, we allowed post-baseline variables (accumulation of comorbidities, disease activity, and prior DMARD use) to enter the regression, updating these variables over time. Variables associated with the use of any TCZ included older age, longer disease duration, seronegative, higher disease activity, no prior use of a TNF inhibitor, and more comorbidities (see Table 4).
Finally, we focused on predictors of TCZm in longitudinal analyses (see Table 3). As with any use of TCZ, predictors of TCZm included older age, longer disease duration, higher disease activity, and no prior use of a TNF inhibitor. The use of TCZc was a very strong predictor of TCZm use. We stratified the population into subsets with no prior TCZc use and ever TCZ use (see Table 3, right columns). After removing prior TCZc users, more comorbidities and non-TNF biologic use demonstrated a significant association with future TCZm use.

Discussion
While there are many good treatment options for RA, the multitude of treatments has led to greater confusion. Which treatments are best for which patients in which order? These questions are complicated to approach but must begin with an understanding of how different treatments are sequenced in RA. We used a large realworld data set from the Corrona RA registry to examine methods for describing the path to a given treatment. We embarked on this sequence analysis focusing on TCZ, and specifically TCZm, as a "proof of concept" study. We found several dominant sequences of treatments that led to TCZm and some patient characteristics associated with TCZm use over time, including prior TCZc use, older age, longer disease duration, seronegative status, higher disease activity, and no immediate past prior use of a TNF inhibitor. Further work will analyze sequences of competing treatments. The goal of these analyses was not to derive a better treatment sequence (we did not focus on clinical outcomes), but rather to develop ways of assessing sequences of DMARDs. The typical method for examining comparative effectiveness in RA has focused on comparing one drug with another, sometimes in randomized controlled trials but often using observational epidemiologic methods [19,20]. While comparing active drugs with one another (versus placebo) is critical, the therapeutic armamentarium for RA has at least five dominant mechanisms each with several agents. As the current analyses demonstrate (Fig. 1 and Table 2), patients switch drugs commonly. The switching to an option like TCZm has strong correlates, including the use of specific prior DMARDs. Thus, to evaluate current practice, sequential trials must be considered. For example, early RA is typically responsive to multiple different therapies, but the effectiveness of treatments after initial therapy varies. Identifying patient characteristics (e.g., serologic status, gender, disease duration, HLA status) that might identify more or less successful treatment sequences would assist clinicians and patients determine the next treatment of choice. However, structuring these decision nodes in a trial format requires an adaptive trial design.
Other implications follow from these findings. First, the differences in predictors between baseline and follow-up were subtle and somewhat to be expected. The variables that became important during follow-up were prior treatments and disease activity, emphasizing the importance of DMARD sequence and DMARD response. Second, the variables associated with TCZ and TCZm use were not substantially different. This likely reflects the fact that almost one-third of patients who try TCZm have used TCZc, possibly because TCZm is similarly effective to TCZc [12][13][14][15]. Third, the path to TCZm typically takes 25.6 months and a median of 1 other bDMARD and up to 5 others. This illustrates the tremendous amount of trial and error that is typical of the treatment course for RA. Several factors may contribute to the relative infrequent use of TCZm as firstline bDMARD. First, in the USA, drug insurance often dictates which drugs are used first-line as a bDMARD. The structure and rules of patients' drug insurance is not a variable contained in the Corrona database. Second, it may be that rheumatologists are less comfortable using TCZm since the drugs are newer to the market than TNFi's or abatacept. Finally, it is also possible that patients are reluctant to use a drug that may worsen their lipid profile.
While this study did not aim to define which treatments are best for a given set of patients, we focused on describing the complex sequence of RA treatments and how patients transition between treatments toward TCZm. TCZm is just one example of a bDMARD that helps illustrate the haphazard treatment sequence that patients with RA experience. We believe that there needs to be a better framework for explicitly testing treatment sequences in RA. The concept of the dynamic treatment regime describes treatment decisions based on patient states that are recognized to evolve over time [21]. For example, early RA is typically responsive to multiple different therapies, but the effectiveness of treatments after initial therapy varies. Identifying patient characteristics (e.g., serologic status, gender, disease duration, HLA status) that might identify more or less successful treatment sequences would assist clinicians and patients determine the next treatment of choice. However, structuring these decision nodes in a trial format The dynamic treatment regimens can be tested in the setting of sequential multiple assignment randomized (SMART) trials [22]. As the name implies, this type of a randomized trial allows for sequences of treatments that may be different across patient subgroups and are determined based on treatment response [23]. SMART trials have grown in popularity but we are not aware of any such trials occurring in RA.
The current set of analyses has important limitations. Similar to most registries, Corrona has longitudinal data but it is "left-censored," people often enter as prevalent cases of RA, and thus, we may not have full information about their treatment history. Corrona is a very large registry but represents patients with RA in rheumatology practices. Thus, patients whose RA is cared for by primary care clinicians are not well represented. As noted, we did not focus on clinical outcomes, making it impossible for us to comment on which treatment is better or worse for a given patient. Like all good science, this type of analysis should be repeated in other datasets.
Strengths of the analysis include the large size of Corrona, and the fact that it reflects real-world evidence. The dataset contains many longitudinal variables, allowing us to consider sequential predictors of DMARD treatment. Including variables that change over time, such as prior treatments and disease activity, had a major impact on the regression coefficients. This suggests that responses to treatment are critical for determining sequences of treatment. Notes: Odds ratios (95% confidence intervals). Baseline models used logistic regression. Models with follow-up data used mixed generalized linear regression TCZ tocilizumab; DMARD disease-modifying anti-rheumatic drug; csDMARDs conventional synthetic DMARDs; nonTNFi nonTCZ bDMARD includes all JAK inhibitors, abatacept, and rituximab ---Not considered at baseline or not significant in univariable screen so not considered in multivariable model or withheld from the model because of problems with convergence In conclusion, we examined sequences of RA treatments in a large USA-based real-world dataset, focusing on TCZm. We characterized time until TCZm, dominant paths to TCZm, and longitudinal predictors of TCZm. This work highlights the variable treatment sequences experienced by patients with RA starting a bDMARD. This variability is likely because of an evidence deficit regarding the comparative effectiveness of different treatment sequences in RA. While observational datasets may provide useful information regarding dynamic treatment regimes, it is more likely that SMAR T trials could substantially impact future care in RA.