Unique challenges and best practices for single cell transcriptomic analysis in toxicology.

The application and analysis of single-cell transcriptomics in toxicology presents unique challenges. These include identifying cell sub-populations sensitive to perturbation; interpreting dynamic shifts in cell type proportions in response to chemical exposures; and performing differential expression analysis in dose-response studies spanning multiple treatment conditions. This review examines these challenges while presenting best practices for critical single cell analysis tasks. This covers areas such as cell type identification; analysis of differential cell type abundance; differential gene expression; and cellular trajectories. Towards enhancing the use of single-cell transcriptomics in toxicology, this review aims to address key challenges in this field and offer practical analytical solutions. Overall, applying appropriate bioinformatic techniques to single-cell transcriptomic data can yield valuable insights into cellular responses to toxic exposures.

Interpretable predictive models of genome-wide aryl hydrocarbon receptor-DNA binding reveal tissue-specific binding determinants.

The aryl hydrocarbon receptor (AhR) is an inducible transcription factor whose ligands include the potent environmental contaminant 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD). Ligand-activated AhR binds to DNA at dioxin response elements (DREs) containing the core motif 5'-GCGTG-3'. However, AhR binding is highly tissue specific. Most DREs in accessible chromatin are not bound by TCDD-activated AhR, and DREs accessible in multiple tissues can be bound in some and unbound in others. As such, AhR functions similarly to many nuclear receptors. Given that AhR possesses a strong core motif, it is suited for a motif-centered analysis of its binding. We developed interpretable machine learning models predicting the AhR binding status of DREs in MCF-7, GM17212, and HepG2 cells, as well as primary human hepatocytes. Cross-tissue models predicting transcription factor (TF)-DNA binding generally perform poorly. However, reasons for the low performance remain unexplored. By interpreting the results of individual within-tissue models and by examining the features leading to low cross-tissue performance, we identified sequence and chromatin context patterns correlated with AhR binding. We conclude that AhR binding is driven by a complex interplay of tissue-agnostic DRE flanking DNA sequence and tissue-specific local chromatin context. Additionally, we demonstrate that interpretable machine learning models can provide novel and experimentally testable mechanistic insights into DNA binding by inducible TFs.

Generative modeling of single-cell gene expression for dose-dependent chemical perturbations.

Single-cell sequencing reveals the heterogeneity of cellular response to chemical perturbations. However, testing all relevant combinations of cell types, chemicals, and doses is a daunting task. A deep generative learning formalism called variational autoencoders (VAEs) has been effective in predicting single-cell gene expression perturbations for single doses. Here, we introduce single-cell variational inference of dose-response (scVIDR), a VAE-based model that predicts both single-dose and multiple-dose cellular responses better than existing models. We show that scVIDR can predict dose-dependent gene expression across mouse hepatocytes, human blood cells, and cancer cell lines. We biologically interpret the latent space of scVIDR using a regression model and use scVIDR to order individual cells based on their sensitivity to chemical perturbation by assigning each cell a "pseudo-dose" value. We envision that scVIDR can help reduce the need for repeated animal testing across tissues, chemicals, and doses.

Prediction of mammalian tissue-specific CLOCK-BMAL1 binding to E-box DNA motifs.

The Brain and Muscle ARNTL-Like 1 protein (BMAL1) forms a heterodimer with either Circadian Locomotor Output Cycles Kaput (CLOCK) or Neuronal PAS domain protein 2 (NPAS2) to act as a master regulator of the mammalian circadian clock gene network. The dimer binds to E-box gene regulatory elements on DNA, activating downstream transcription of clock genes. Identification of transcription factor binding sites and genomic features that correlate to DNA binding by BMAL1 is a challenging problem, given that CLOCK-BMAL1 or NPAS2-BMAL1 bind to several distinct binding motifs (CANNTG) on DNA. Using three different types of tissue-specific machine learning models with features based on (1) DNA sequence, (2) DNA sequence plus DNA shape, and (3) DNA sequence and shape plus histone modifications, we developed an interpretable predictive model of genome-wide BMAL1 binding to E-box motifs and dissected the mechanisms underlying BMAL1-DNA binding. Our results indicated that histone modifications, the local shape of the DNA, and the flanking sequence of the E-box motif are sufficient predictive features for BMAL1-DNA binding. Our models also provide mechanistic insights into tissue specificity of DNA binding by BMAL1.

Pharmacokinetic and toxicodynamic concepts in idiosyncratic, drug-induced liver injury.

INTRODUCTION: Idiosyncratic drug-induced liver injury (IDILI) causes morbidity and mortality in patients and leads to curtailed use of efficacious pharmaceuticals. Unlike intrinsically toxic reactions, which depend on dose, IDILI occurs in a minority of patients at therapeutic doses. Much remains unknown about causal links among drug exposure, a mode of action, and liver injury. Consequently, numerous hypotheses about IDILI pathogenesis have arisen. AREAS COVERED: Pharmacokinetic and toxicodynamic characteristics underlying current hypotheses of IDILI etiology are discussed and illustrated graphically. EXPERT OPINION: Hypotheses to explain IDILI etiology all involve alterations in pharmacokinetics, which lead to plasma drug concentrations that rise above a threshold for toxicity, or in toxicodynamics, which result in a lowering of the toxicity threshold. Altered pharmacokinetics arise, for example, from changes in drug metabolism or from transporter polymorphisms. A lowered toxicity threshold can arise from drug-induced mitochondrial injury, accumulation of toxic endogenous factors or harmful immune responses. Newly developed, interactive freeware (DemoTox-PK; https://bit.ly/DemoTox-PK) allows the user to visualize how such alterations might lead to a toxic reaction. The illustrations presented provide a framework for conceptualizing idiosyncratic reactions and could serve as a stimulus for future discussion, education, and research into modes of action of IDILI.

A Negative Feedback Loop and Transcription Factor Cooperation Regulate Zonal Gene Induction by 2, 3, 7, 8-Tetrachlorodibenzo-p-Dioxin in the Mouse Liver.

The cytochrome P450 (Cyp) proteins Cyp1A1 and Cyp1A2 are strongly induced in the mouse liver by the potent environmental toxicant 2, 3, 7, 8-tetrachlorodibenzo-p-dioxin (TCDD), acting through the aryl hydrocarbon receptor (AHR). The induction of Cyp1A1 is localized within the centrilobular regions of the mouse liver at low doses of TCDD, progressing to pan-lobular induction at higher doses. Even without chemical perturbation, metabolic functions and associated genes are basally zonated in the liver lobule along the central-to-portal axis. To investigate the mechanistic basis of spatially restricted gene induction by TCDD, we have developed a multiscale computational model of the mouse liver lobule with single-cell resolution. The spatial location of individual hepatocytes in the model was calibrated from previously published high-resolution images. A systems biology model of the network of biochemical signaling pathways underlying Cyp1A1 and Cyp1A2 induction was then incorporated into each hepatocyte in the model. Model simulations showed that a negative feedback loop formed by binding of the induced Cyp1A2 protein to TCDD, together with cooperative gene induction by the ß-catenin/AHR/TCDD transcription factor complex and ß-catenin, help produce the spatially localized induction pattern of Cyp1A1. Although endogenous WNT regulates the metabolic zonation of many genes, it was not a driver of zonal Cyp1A1 induction in our model. Conclusion: In this work, we used data-driven computational modeling to identify the mechanistic basis of zonally restricted gene expression induced by the potent and persistent environmental pollutant TCDD. The multiscale model and derived results clarify the mechanisms of dose-dependent hepatic gene induction responses to TCDD. Additionally, this work contributes to our broader understanding of spatial gene regulation along the liver lobule.

Pregnancy-specific physiologically-based toxicokinetic models for bisphenol A and bisphenol S.

Predictions from physiologically based toxicokinetic (PBTK) models can help inform human health risk assessment for potentially toxic chemicals in the environment. Bisphenol S (BPS) is the second most abundant bisphenol detected in humans in the United States, after bisphenol A (BPA). We have recently demonstrated that BPS, much like BPA, can cross the placental barrier and disrupt placental function. Differences in physicochemical properties, toxicokinetics, and exposure outcomes between BPA and other bisphenols prevent direct extrapolation of existing BPA PBTK models to BPS. The current study aimed to develop pregnancy-specific PBTK (p-PBTK) models for BPA and BPS, using a common p-PBTK model structure. Novel paired maternal and fetal pregnancy data sets for total, unconjugated, and conjugated BPA and BPS plasma concentrations from three independent studies in pregnant sheep were used for model calibration. The nine-compartment (maternal blood, liver, kidney, fat, placenta and rest of body, and fetal liver, blood and rest of body) models simulated maternal and fetal experimental data for both BPA and BPS within one standard deviation for the majority of the experimental data points, highlighting the robustness of both models. Simulations were run to examine fetal exposure following daily maternal exposure to BPA or BPS at their tolerable daily intake dose over a two-week period. These predictive simulations show fetal accumulation of both bisphenols over time. Interestingly, the steady-state approximation following this dosing strategy achieved a fetal concentration of unconjugated BPA to levels observed in cord blood from human biomonitoring studies. These models advance our understanding of bisphenolic compound toxicokinetics during pregnancy and may be used as a quantitative comparison tool in future p-PBTK models for related chemicals.

Quadriceps mechano- and electromyographic time-frequency responses during muscular contractions to volitional exhaustion.

INTRODUCTION: Surface electromyography (SEMG) and mechanomyography (SMMG) responses of the quadriceps during muscular contractions to exhaustion were computed and analyzed by analysis of variance and polynomial regression analyses. METHODS: Participants performed maximum flexion-extension movements at 180°/s until volitional exhaustion, rested for 2 minutes, and then completed a second bout of movements until exhaustion. Torque and SEMG/SMMG median frequencies and amplitudes were examined at 9 points across repetitions completed. RESULTS: (1) Torque decreased precipitously; (2) SEMG amplitude displayed an initial increase, then a steady decrease, and SMMG amplitude showed a continuous decrease; and (3) SEMG and SMMG median frequencies displayed a continual decrease over repetitions completed. Fractional polynomial and quadratic models explained the fatigue process with the highest precision. CONCLUSIONS: Changes in electrical and mechanical properties of the quadriceps during fatigue reflect alterations in neuromuscular activation strategies and/or muscle wisdom. SEMG frequency modeled muscle fatigue more effectively than amplitude, whereas SMMG frequency and amplitude were equally effective.

Wavelet analysis of quadriceps power spectra and amplitude under varying levels of contraction intensity and velocity.

INTRODUCTION: We investigated the effect of contraction intensity [100%, 75%, 50%, and 25% maximum voluntary contraction (MVC)] and movement velocity (50°, 100°, 200°, and 400°/s) on surface electromyography root mean square amplitude (SEMGRMS ) and median frequency (SEMGMDF ) of rectus femoris (RF), vastus lateralis (VL), and vastus medialis (VM). METHODS: SEMGs during knee extension were resolved into their respective frequencies using wavelet transformations. RESULTS: RF, VL, and VM muscles displayed increased SEMGMDF as contraction intensity increased from 25% to 50% MVC and from 75% to 100% MVC, and each muscle displayed its own unique frequency shifting patterns. The SEMGMDF was not influenced by movement velocity. SEMGRMS increased in all 3 muscles as contraction intensity increased and was influenced by movement velocity, with the highest values observed at 400° and 200°/s. CONCLUSIONS: We infer that increasing contraction intensity facilitates greater recruitment of fast-twitch muscle fibers, but there are differing responses in RF, VL, and VM muscles.