Rheumatoid arthritis (RA) is a chronic systemic autoimmune disease characterized by symmetric joint inflammation, progressive cartilage and bone destruction, and potential extra-articular manifestations [2]. The management of RA has undergone significant changes in both therapeutic targets and strategies in recent decades. In recent years, the prospects of achieving drug-free remission, as well as the use of genetically engineered biological drugs and Janus kinase inhibitors, have been actively discussed [3]. Despite this, methotrexate (MTX) remains the cornerstone disease-modifying antirheumatic drug (DMARD) for RA therapy [4]. Its widespread use is attributed not only to its clinical efficacy but also to its cost-effectiveness.

MTX is recommended by leading international rheumatology associations as a first-line agent for both monotherapy and in combination with other DMARDs [4, 5]. However, MTX efficacy is limited in 30% of patients due to an insufficient therapeutic response or side effects necessitating drug discontinuation [1]. The lack of tools to predict the response to therapy remains a significant challenge.

A key aspect of the "treat-to-target" strategy is the early initiation of DMARD therapy, particularly within the first six months of RA onset, which significantly increases the likelihood of long-term disease remission. Consequently, the early identification of risk factors for MTX inefficiency is crucial. An active search is underway for clinical, laboratory, and genetic markers capable of predicting the response to MTX. These include demographic characteristics such as sex and age, as well as smoking status, body weight, immunological parameters, and others [6–9]. Single nucleotide polymorphisms (SNPs) in genes encoding proteins involved in the pharmacokinetics and pharmacodynamics of MTX are being investigated (Fig. 1). The present study aims to develop a pharmacogenetic model for predicting the response to MTX in RA. While the creation and practical application of predictive models are used in oncology, psychiatry, and other areas of medical practice, this approach has not been widely adopted in clinical rheumatology, including for the treatment of the most common autoimmune disease, RA. A major limitation for implementing pharmacogenetic research is the lack of laboratory tests for polymorphism typing.

Fig. 1. Mechanism of action of methotrexate, adapted from [10]

Note: The proteins whose genes are studied in the current research are highlighted in red.

МТ (methotrexate); SLC19A1 (solute carrier family 19 member 1); ABCB1 (ATP Binding Cassette Subfamily B Member 1); FPGS (folylpolyglutamate synthase); GGH (gamma-glutamyl hydrolase); МТPG — polyglutamated МТ; DHFR (dihydrofolate reductase); THF (tetrahydrofolate); 5,10-СН2-THF (5,10-methylentetrahydrofolate); MTHFR (methylenetetrahydrofolate reductase); 5-СН3-THF (5-methyltetrahydrofolate); DHF (dihydrofolate); MS (methionine synthase); TS (thymidylate synthase); dUMP (deoxyuridine monophosphate); dTMP (deoxythymidine monophosphate); ATIC (5-aminoimidazole-4-carboxamide ribonucleotide transformylase); 10-CHO-THF (methylenetetrahydrofolate dehydrogenase); AICAR (5-aminoimidazole-4-carboxamide ribonucleotide); FAICAR (5-formylaminoimidazole-4-carboxamide ribonucleotide); ITPA (inosine triphosphatase); ITP (inosine triphosphate); IMP (inosine monophosphate); ATP (adenosine trihosphate); AMP (adenosine monophosphate); AMPD1 (adenosine monophosphate deaminase); ADA (adenosine deaminase); cAMP (cyclic adenosine monophosphate); IL-10 (interleukin-10).

The study was approved by the Local Ethics Committee of the Chelyabinsk State Medical Academy of the Ministry of Health of the Russian Federation (Protocol No. 10, November 25, 2012).

The study included 281 patients with a confirmed diagnosis of RA. It was designed as a prospective cohort study with participant enrollment over a ten-year period. After six months of MTX treatment, therapy efficacy was assessed using the Disease Activity Score-28 (DAS28), identifying groups of "responders" (170 patients) and "non-responders" (111 patients). Using a candidate gene approach, 12 SNPs from nine genes involved in MTX metabolism were selected. SNPs were chosen based on the following criteria: all are registered in the dbSNP database and have confirmed associations with MTX efficacy published in PubMed or PharmGKB. The minor allele frequency of each SNP in the population was at least 5%, ensuring statistical relevance.

Genomic DNA was extracted from peripheral venous blood samples using the commercial kit "Protrans DNA Box 500" (Germany). Genotyping was performed by polymerase chain reaction (PCR). For 10 polymorphisms—SLC19A1 SNP rs1051266 (G80A); ABCB1 SNP rs1128503 (C1236T) and rs2032582 (A2677C); GGH SNP rs3758149 (C-401T); FPGS SNP rs84451422 (A>C) and rs1544105 (C105T); ATIC SNP rs2372536 (C347G); ADA SNP rs244076 (T>C); AMPD1 SNP rs17602729 (C34T); ITPA SNP rs1127354 (C94A)—primers were domestically developed and produced by TestGen LLC (Russia) for the first time. For the MTHFR SNPs rs1801131 (A1298C) and rs1801133 (C677T), ready-made primers were used. Detection of amplification products was performed using endpoint FLASH analysis on QuantStudio amplifiers (Applied Biosystems). Data interpretation was carried out using QuantStudio Design and Analysis Software (version 1.5.2).

Statistical analysis was performed using multifactor dimensionality reduction (MDR) and information analysis (Shannon entropy), resulting in several predictive models for MTX efficacy.

The development of predictive models requires careful selection of indicators to ensure high accuracy and a minimal number of variables. Model reliability implies its ability to accurately predict outcomes both in the studied sample and in other data from the general population. Optimality is achieved by using the smallest number of predictors that maintain high predictive power. The MDR method was applied to identify genes and their interactions affecting the risk of non-response to MTX in RA patients. The process of creating a predictive model involved several stages.

Automated Analysis — This was the initial stage in the search for optimal pharmacogenetic models, covering one to three genes (Table 1, Fig. 2, Fig. 3).

Fig. 3. Column charts of the number of patients in genotype combination cells for a number of two-gene models

Notes: Left-hand columns are the number of non-responders, right-hand columns are the number of responders to therapy; dark-gray cells are high-risk combinations, light-gray cells are low-risk combinations, white cells are no genotype combinations.

Table 1. Models for automatic prediction of non-response to therapy based on the results of MDR analysis (n_non-response=111, n_response=170)

ModelSensitivity, % [95% CI]Specificity, % [95% CI]Diagnostic AccuracyOR [95% CI]P-valueCross-validation ConsistencyATIC rs237253657.7 [48.4; 66.6]58.8 [51.3; 66.0]58.21.95 [1.20; 3.16]χ²(1)=7.31, p=0.00710/10ATIC rs2372536 +SLC19A1 rs105126665.8 [56.6; 74.1]55.3 [47.8; 62.6]60.52.38 [1.45; 3.90]χ²(1)=11.96, p<0.0016/10ATIC rs2372536 +SLC19A1 rs1051266 +MTHFR rs180113373.0 [64.2; 80.6]58.2 [50.7; 65.5]65.63.76 [2.24; 6.32]χ²(1)=26.34, p<0.0015/10Selected Two-Gene Models      ATIC rs2372536 +MTHFR rs180113359.5 [50.2; 68.3]59.4 [51.9; 66.6]59.42.15 [1.32; 3.49]χ²(1)=9.58, p=0.00210/10ADA rs244076 +MTHFR rs180113153.2 [43.9; 62.3]57.7 [48.4; 66.6]55.41.54 [0.95; 2.50]χ²(1)=3.15, p=0.07610/10SLC19A1 rs1051266 +FPGS rs445142245.0 [36.0; 54.3]72.9 [65.9; 79.2]59.02.21 [1.33; 3.66]χ²(1)=9.66, p=0.00210/10

The ATIC rs2372536 polymorphism showed the highest predictive significance as a single-gene model, demonstrating a sensitivity of 57.7% (95% CI: 48.4;66.6), specificity of 58.8% (95% CI: 51.3;66.0), diagnostic accuracy of 58.2%, and complete stability in tenfold cross-validation (10/10) (Table 1, Fig. 2A).

The next automatically generated model included the combination "ATIC rs2372536 + SLC19A1 rs1051266". This association demonstrated a moderate improvement in diagnostic accuracy with a sensitivity of 65.8%, but low cross-validation consistency (6/10), limiting its applicability to other samples (Table 1, Fig. 2B). The addition of MTHFR rs1801133 increased sensitivity to 73.0% and the odds ratio to 3.76. However, the cross-validation consistency of the "ATIC rs2372536 + SLC19A1 rs1051266 + MTHFR rs1801133" model was low (5/10), reducing its reliability for practical application (Table 1, Fig. 2C).

Fig. 2. Column charts of the number of patients in cells of genotype combinations for three (A, B and C) models of automatic construction

Notes: Columns on the left are the number of non-responders, columns on the right are the number of responders to therapy; dark gray cells are combinations of increased risk, light gray are of decreased risk, white are no combination of genotype combinations.

The automatically created two-gene models—"ATIC rs2372536 + MTHFR rs1801133", "ADA rs244076 + MTHFR rs1801131", and "SLC19A1 rs1051266 + FPGS rs4451422"—also with high consistency (10/10), did not demonstrate sufficient sensitivity and specificity, with diagnostic accuracies of 59.4, 55.4, and 59.0, respectively (Table 1, Fig. 3).

Thus, the automated search for optimal models failed to create a stable model with high diagnostic accuracy, as the approach relied solely on mathematical criteria without considering biochemical data.

Information Analysis. To improve the quality of predictive models, an entropy graph based on Shannon entropy was constructed at the next stage (Fig. 4), reflecting the contribution of genes to the lack of response to therapy. This measure of uncertainty allowed for the estimation of information gain as the difference between probability distributions of the system with and without certain elements. The graph consists of vertices (genes) and edges (their interactions), where information gain values (in percent) reflect the contribution of genes and their pairwise interactions to the total entropy. Edge thickness corresponds to the magnitude of the gain, and color indicates its nature: orange and red indicate synergistic, non-additive (epistatic) interaction that enhances the effect of genes; green and blue indicate additive interaction with information redundancy; brown lines indicate weak or independent influence. The combination of ATIC rs2372536 and SLC19A1 rs1051266 provided a gain of 3.11% (2.03% + 0.37% + 0.71%), surpassing the combination of SLC19A1 rs1051266 and FPGS rs4451422 (2.46%). The predominance of orange and red edges on the graph emphasizes the dominance of epistatic interactions, particularly between SLC19A1 rs1051266 and the FPGS polymorphisms rs4451422 and rs1544105.

Fig. 4. Entropy graph (top) and similarity dendrogram (bottom) for the contribution of genes and their interactions to the risk of non-response to therapy

The similarity dendrogram (Fig. 4) revealed close non-additive interaction between SLC19A1 rs1051266 and FPGS rs4451422, as well as a cluster of GGH rs3758149 and FPGS rs1544105. These connections formed the basis for a biochemically substantiated polyglutamylation model, including polymorphisms encoding the proteins for intracellular transport (SLC19A1), polyglutamylation (FPGS), and deglutamylation (GGH). The enzymes encoded by these genes are sequentially involved in the intracellular transport and metabolism of MTX, determining the concentration of its active form. The epistatic interactions of these genes are reflected in a model with high statistical significance (p<0.001) and stability (Table 2). However, the diagnostic accuracy was 65.1% with low specificity (60.0%), which is insufficient for precise prediction.

The final model combined two pairs of genes encoding the adenosine and folate pathways of MTX metabolism: ATIC rs2372536 with MTHFR rs1801133, and ADA rs244076 with MTHFR rs1801131. Informational redundancy was detected within the pairs, indicating additive interaction related to biochemical processes. MTHFR and ATIC are key enzymes determining the mechanisms of MTX action: regulating the synthesis of folates, adenosine, and purines. The SLC19A1 rs1051266 gene, included as a fifth component, confirmed the significant role of enzymes in MTX intracellular transport.

Table 2. Models for predicting non-response to therapy based on the results of MDR analysis, created on the basis of biological information (n_non-response=111, n_response=170)

ModelSensitivity, % [95% CI]Specificity, % [95% CI]Diagnostic AccuracyOR [95% CI]P-valueCross-validation ConsistencyPolyglutamylation ModelSLC19A1 rs1051266 +FPGS rs4451422 +FPGS rs1544105 +GGH rs375814970.3 [61.3; 78.2]60.0 [52.5; 67.1]65.13.55 [2.13; 5.90]χ²(1)=24.65, p<0.00110/10Final ModelATIC rs2372536 +ADA rs244076 +MTHFR rs1801133 +MTHFR rs1801131 +SLC19A1 rs105126680.2 [72.0; 86.8]69.4 [60.4; 77.4]74.89.18 [5.19; 16.22]χ²(1)=66.06, p<0.00110/10

Using MDR, the contribution of 12 SNPs from 9 genes to the risk of non-response to MTX was analyzed. The final model "ATIC rs2372536 + MTHFR rs1801133 + ADA rs244076 + MTHFR rs1801131 + SLC19A1 rs1051266" showed high statistical significance (p<0.001) and cross-validation consistency (10/10) (Table 2). Its "if-then" decision rule provides a risk prediction with a sensitivity of 80.2% (95% CI: 72.0;86.8) and specificity of 69.4% (95% CI: 60.4;77.4). The practical application of the final model is presented in a fragment of the "if-then" decision rule table (Table 3). The model is patented and ready for practical application [11].

Table 3. The decision rule "if, then" for the final model of predicting the risk of non-response to methotrexate by genetic indicators: 1 — non-response, 0 — response (fragment)

IfAndAndAndAndThen (Risk =)SLC19A1rs1051266ADArs244076ATICrs2372536MTHFRrs1801131MTHFRrs1801133 GGTTCCCCCC1GGTTCCAACC0GGTTCCAACT0GGTTCCAATT0GGTTCCACCC1GGTTCCACCT0GGTTGGCCCC1GGTTGGAACC0

Over the last decade, medical science has been striving for new ways to improve the efficacy of treatment for various diseases. In the Russian Federation, the direction "Transition to personalized medicine, high-tech healthcare, and health-saving technologies, including through the rational use of medications" has been recognized as a priority [12]. In rheumatology, such tasks concern the therapy of RA as the most common and socially significant disease. A personalized approach to treatment can be realized by studying the genes regulating the metabolism of MTX, as the first-line pharmaceutical for RA. Our team's research activity in the pharmacogenetics of MTX efficacy and tolerability in RA has been ongoing since 2013. During this time, we have investigated the complex mechanism of MTX action and obtained data on certain genetic polymorphisms influencing the therapeutic response [13–15]. Understanding the complexity of MTX's effects served as an incentive for further work and the implementation of pharmacogenetic testing in practice. For this purpose, a domestically developed validated test system was preliminarily created, and key biomarkers of therapeutic response to MTX were predetermined using a candidate gene approach. The next stage was the search for a predictive pharmacogenetic model. The first attempt to create such a model was based on an automated (machine) analysis without considering the principles of MTX metabolism. This method did not yield the expected results. The second attempt at model creation was based on the non-additive interaction of genes, taking into account MTX metabolism. The resulting final model showed good predictive value and stability upon cross-validation. The result of this stepwise approach was the creation of a computer program for deciding on the advisability of initiating baseline treatment with MTX. Since the data were obtained from a European population, the reliability and reproducibility of the results across most territories of the Russian Federation are likely to be high.

Thus, the decade-long study on predicting MTX efficacy in RA has led to: the algorithmization of the process for building a drug efficacy prediction model; the development of a domestic test system for genetic typing of "ATIC rs2372536 + MTHFR rs1801133 + ADA rs244076 + MTHFR rs1801131 + SLC19A1 rs1051266"; and the practical implementation of the research results via a computer program.

A limitation of the study is the moderate specificity of the final model (69.4%), which may reduce its predictive value in some clinical scenarios. The limited sample size and its single-center nature may not fully reflect the genetic heterogeneity of the population. The lack of accounting for epigenetic factors and gene expression data also limits the comprehensiveness of the analysis.

The final model, combining five polymorphisms—"ATIC rs2372536 + MTHFR rs1801133 + ADA rs244076 + MTHFR rs1801131 + SLC19A1 rs1051266"—demonstrates high sensitivity and stability, opening prospects for a personalized treatment approach. Despite limitations related to moderate specificity and the single-center design, the results of the pharmacogenetic approach create a foundation for further research aimed at optimizing RA therapy.

Accounting for biochemical processes in the development of pharmacogenetic models allows for a deeper understanding of gene interactions, which is crucial for improving the diagnostic accuracy of the final model. Ultimately, this increases the precision of predicting the therapeutic response to MTX. A domestically produced test system was developed to implement the study, which will allow for the introduction of pharmacogenetic typing into clinical practice.