Background

Methotrexate (MTX) remains the anchor drug in rheumatoid arthritis (RA), yet inter-patient pharmacokinetic (PK) variability leads to unpredictable toxicity — hepatotoxicity, cytopenias, and pneumonitis — affecting 10–30% of patients. Traditional population PK models assume log-normal parameter distributions, poorly capturing the heavy-tailed, multimodal exposure distributions observed in polyglutamated MTX concentrations across patients with variable MTHFR, ABCB1, and SLC19A1 genotypes.

Hypothesis

Normalizing flow (NF) models — invertible neural networks that learn arbitrary probability distributions via successive bijective transformations — applied to serial intracellular MTX-polyglutamate (MTX-PGn) concentrations from dried blood spots will:

1. Capture non-Gaussian, multimodal PK variability across patient subpopulations defined by pharmacogenomic haplotypes (MTHFR C677T/A1298C, ABCB1 C3435T, SLC19A1 G80A)
2. Identify latent toxicity-predisposing exposure trajectories by detecting when a patient's PK trajectory migrates toward high-density regions of the learned toxicity distribution
3. Predict grade ≥2 MTX toxicity events 4–10 weeks before clinical manifestation with AUROC >0.82, outperforming standard population PK/PD models (expected AUROC ~0.65)

Theoretical Framework

Let x = (MTX-PG₁, MTX-PG₂, MTX-PG₃, MTX-PG₄, MTX-PG₅) represent the polyglutamation profile at each sampling timepoint. A normalizing flow f = f_K ∘ f_{K-1} ∘ ... ∘ f_1 maps x to a base Gaussian distribution z through K affine coupling layers. The exact log-likelihood is computed via the change-of-variables formula:

log p(x) = log p(z) + Σ log |det(∂f_k/∂x_k)|

Conditioning on pharmacogenomic covariates g (one-hot encoded haplotypes) via conditional normalizing flows p(x|g) allows the model to learn genotype-specific PK distributions. Toxicity prediction emerges by monitoring the Kullback-Leibler divergence between a patient's rolling PK trajectory distribution and the learned toxicity-associated density region.

Testable Predictions

1. NF-estimated PK distributions will exhibit ≥3 distinct modes corresponding to MTHFR/ABCB1 haplotype combinations (vs. unimodal assumption in standard models)
2. Patients whose serial MTX-PGn profiles show KL divergence >2.0 nats toward the toxicity-associated density region will develop grade ≥2 toxicity within 10 weeks (sensitivity >80%, specificity >75%)
3. The NF model will identify a previously unrecognized high-risk PK phenotype (estimated 8–12% of patients) characterized by paradoxically low total MTX-PGn but elevated PG₄/PG₅ ratio, predicting hepatotoxicity independent of cumulative dose
4. Conditional NFs stratified by SLC19A1 G80A will reveal that GG homozygotes require 15–25% lower MTX doses for equivalent therapeutic exposure, currently masked by population-average dosing

Proposed Validation

• Dataset: Prospective serial dried blood spot MTX-PGn sampling (biweekly × 24 weeks) from 400 RA patients initiating MTX, with paired pharmacogenomic genotyping
• Architecture: RealNVP with 8 affine coupling layers, conditional on genotype embeddings
• Benchmarks: Population PK (NONMEM), Bayesian hierarchical PK, and Gaussian mixture models
• Primary endpoint: Time-dependent AUROC for grade ≥2 toxicity prediction at 4, 6, 8, and 10 weeks prior to event
• Calibration: Expected calibration error (ECE) <0.05 across deciles

Limitations

• Requires serial MTX-PGn monitoring infrastructure (dried blood spots with LC-MS/MS), not yet standard of care
• NF models are computationally intensive and may require GPU inference for real-time clinical use
• Training requires sufficient toxicity events (~60–80) for stable density estimation in the toxicity region
• Pharmacogenomic conditioning assumes additive haplotype effects; epistatic interactions may require graph-based conditioning architectures
• External validation across diverse ancestral populations essential before clinical deployment

Clinical Significance

If validated, this approach transforms MTX monitoring from reactive (detect toxicity after it occurs) to predictive (identify toxicity-prone PK trajectories weeks in advance). The genotype-conditional density estimation enables precision dosing that accounts for the full distributional complexity of drug exposure — not just mean levels — potentially reducing MTX discontinuation rates (currently 30–40% at 2 years) while maintaining therapeutic efficacy.

RheumaAI Research • rheumai.xyz • DeSci Rheumatology