A Multi-Scale Computational Analysis of Bucillamine Neuroprotection Against Soman Toxicity

Abstract

This computational framework evaluates the neuroprotective mechanism of the dithiol antioxidant Bucillamine (and its active metabolite SA981) against Soman-induced oxidative damage. By integrating quantum chemical calculations, classical molecular dynamics, and systems-level ordinary differential equation (ODE) modeling, this work establishes a multi-scale predictive chain that identifies catalytic longevity as a key metric for antioxidant efficacy under acute nerve-agent stress. The framework is designed to generate falsifiable predictions testable by biological assay, not clinically predictive estimates. Key Methodological Layers: - Transport Energetics: Human LAT1 (PDB 7DSQ) and xCT (PDB 7P9U) channel binding affinities are computed using 18,300 ns of umbrella sampling PMF simulations (61 windows × 50 ns per window, 6 compounds). All tested compounds show barrierless LAT1 entry; the physically meaningful metric is well depth, reflecting substrate recognition. - QM/MM Disulfide-Exchange Energetics: The reaction coordinate for Somatostatin (SST) disulfide cleavage is modeled using gas-phase Density Functional Theory (DFT) relaxed scans at the PBE-D3/DZVP level of theory, comparing intact-disulfide and free-thiolate nucleophiles on the same 103-atom cluster. - Systems Pharmacology Modeling: A 28-variable ODE model of Soman neurotoxicity and antioxidant kinetics is constructed and calibrated to empirical literature targets using Bayesian Markov Chain Monte Carlo (MCMC) parameter estimation with the emcee ensemble sampler (7 parameters, 6 calibration targets, all within 1σ). Extensions model benzodiazepine pharmacoresistance via ATP-dependent seizure thresholds and GABA(A) surface-dependent diazepam efficacy, and NOX-derived ROS source decomposition. - Global Sensitivity & Counterfactual Analysis: Global sensitivity analysis via Sobol variance decomposition (N=22,528N = 22{,}528 N=22,528 model evaluations) isolates the primary drivers of ROS variance. Computational counterfactuals systematically disable individual pathway modules to map the causal hierarchy of neuroprotection. Core Scientific Findings: - The Catalytic Longevity Paradigm: Due to its dithiol structure, a single Bucillamine molecule can undergo an estimated 500 to 5,000 ROS-scavenging cycles under severe stress before irreversible overoxidation, compared to only ~11 cycles for single-thiol counterparts like cysteamine. Sobol analysis confirms that two parameters — scavenging rate and insult magnitude — account for ~85% of ROS variance, with all other parameters contributing negligibly. - The Reduced-Species Route: Direct exchange between intact cyclic disulfides is highly barrier-limited (barrier ≥\geq ≥ 40 kcal/mol). In contrast, the reduced Bucillamine thiolate presents a dramatically lower barrier (~2.5 kcal/mol), demonstrating that intracellular reduction is a strict mechanistic prerequisite for downstream somatostatin depletion. SST depletion itself is identified as a tolerable side effect — not a protective mechanism — that increases excitability via glutamate disinhibition and Ca²⁺-driven GABA(A) internalization. - GABA(A) Trafficking Is Ca²⁺-Driven: Endocytosis decomposition at calibrated MAP parameters attributes ~93% of GABA(A) receptor internalization to Ca²⁺-dependent signaling and ~0.2% to ROS, consistent with the calcineurin-mediated trafficking pathway established in the status epilepticus literature. SA981's primary therapeutic value operates through mitochondrial and protein sulfhydryl protection rather than GABA preservation. - Benzodiazepine Pharmacoresistance: The model predicts that SA981 prevents the transition from benzodiazepine-responsive to refractory status epilepticus by preserving ATP (1.00 vs 0.39 under standard care at 240 min), stabilizing the seizure threshold rather than preventing GABA(A) surface loss. This generates a falsifiable prediction: maintained diazepam responsiveness at late timepoints should correlate with ATP preservation, not GABA(A) receptor density. - The Therapeutic Window: Calibrated GSH-dependent dynamics predict that the catalytic cycle maintains high efficacy for 4 to 5 hours before GSH depletion terminates the cycle, with protein sulfhydryl levels remaining above 55% through 210 minutes of treatment delay, validating the clinical potential of delayed antioxidant intervention. NOX source decomposition (66% mitochondrial, 34% NOX-derived) suggests that combining SA981 with NOX inhibitors could extend the protective window beyond 4 hours by reducing the sustained enzymatic ROS production that eventually overwhelms the catalytic cycle.