Non-animal testing to assess the efficacy and safety of chemical substances is one of the key concepts in balancing animal welfare and efficient development. Since the marketing ban in the EU in March 2013 [(EC) No. 1223/2009] (EU, 2009) of cosmetic products and ingredients tested on animal models, safety assessment methodologies independent of animal testing have attracted much attention. Simultaneously, the utilization of non-animal high-throughput technology for optimizing drug discovery processes is becoming highly important in pharmaceuticals (Loiodice et al., 2017; Rognan, 2017; Amano et al., 2020).

Read-across, a process that estimates substance toxicity based on the concept that substances with similar chemical structures have similar biological activities, represents a promising approach and has already been conceptually accepted as a reliable safety risk assessment by some regulatory authorities (ECHA, 2017; European Commission, 2018). Likewise, quantitative structure–activity relationship (QSAR) has been widely used, and impurity characterization received regulatory acceptance (ICH M7). However, since subtle structural differences may elicit different biological responses, supporting the read-across robustness by using biological similarities has been considered important (Ball et al., 2016, 2020; Zhu et al., 2016). Registration, Evaluation, Authorization, and Restriction of Chemicals (REACH) mentions that the read-across performed by registrants often fails to comply with the legal requirements due to defects in the hypothesis and justification of toxicological prediction (ECHA, 2020).

There are two approaches to enhance the reliability of read-across: 1) employment of in vitro data relevant to specific toxicity. Methodologies to incorporate in vitro data within read-across (Ball et al., 2016, 2020; ECHA, 2017; Guo et al., 2019) and some case studies (OECD, 2016a, 2016b, 2018; Nakagawa et al., 2020, 2021) have been reported. However, these approaches can be applied only to specific toxicity end points and substances with a known toxicity and mode of action. Such conditions were previously termed as “local validity” (Patlewicz et al., 2014). 2) The use of biologically similar substances based on their profiles obtained from a large number of bioassays. The United States Environmental Protection Agency’s (US EPA’s) research project, ToxCast and Tox21, provided hundreds of high-throughput screening assays, and several groups employed such biological activity data for toxicological evaluation (Sipes et al., 2013; Berggren et al., 2015; Richard et al., 2021). Although this concept could be applied to substances with little information to elucidate their entire toxicological profiles and find their key mode of action, it is time-consuming and expensive to conduct numerous bioassays for a new candidate substance. By contrast, transcriptome data containing approximately 30,000 gene expression values can be used to estimate perturbated mechanisms through enrichment analysis. Wang et al. (2016) tried to predict drug-induced adverse effects by employing LINCS L1000 data (Subramanian et al., 2017), whereas Iwata et al. (2019) developed a computational method to predict missing values from the LINCS L1000 transcriptomic profiles of various human cell lines and provided new drug therapeutic indications. Genomic data have been considered to be usable in read-across by Health Canada and a research group from the US FDA (Health Canada, 2019; Liu et al., 2019). However, several researchers have shown that in vitro gene expression values are not always highly correlated with in vivo data (Sutherland et al., 2016; Grinberg et al., 2018; Liu et al., 2018). Thus, interpreting toxicological meaning from the in vitro–in vivo relationship and in vitro to in vivo extrapolation (IVIVE) in omics data represents a big challenge for chemical risk assessment. IVIVE was originally researched in toxicokinetics, such as in hepatic clearance and metabolites using hepatocytes (Soars et al., 2007; Umehara and Camenisch, 2012); most recent studies on non-animal testing have focused on predicting plasma concentrations, which is relevant for identification of a margin of exposure in risk assessment (Thomas et al., 2013a; Bell et al., 2018; Li et al., 2021). However, IVIVE should be considered for both toxicokinetics and toxicodynamics. Understanding of the in vitro to in vivo relationship of bioactivity data is also essential for non-animal testing. As an IVIVE study in omics data, Liu et al. (2020) developed a useful in silico strategy to narrow the data gap between in vitro and in vivo conditions. They modified in vitro data using non-generative matrix factorization methods to improve the correlation with in vivo data, which overcame the shortcomings of previous large-scale genomic data predictions regarding the in vitro–in vivo data gap (Liu et al., 2020). Although non-generative matrix factorization enables macroscopic estimation based on a pattern recognition classifying chemical and biological responses, it does not focus on estimation of each gene. As an alternative solution, microscopic estimation of each gene expression was performed based on tensor-train weighted optimization using machine learning (Iwata et al., 2019); however, such comprehensive estimations have not been integrated within an IVIVE study. Therefore, predicting in vivo transcriptomic profiles from in vitro data for IVIVE might not only enhance the robustness of read-across but could also be utilized in other non-animal testing strategies as weight of evidence, such as in Integrated Approaches to Testing and Assessment (IATA) and new approach methods (NAMs) for safety and drug repositioning research.

In this study, we developed a virtual DNA microarray that quantitatively predicts the in vivo gene expression profiles based on the chemical structure and/or in vitro transcriptome data. For each gene, elastic-net models, a regression analysis method that has been used in toxicity prediction with visualization of feature importance (e.g., Fujita et al., 2020), were constructed using chemical descriptors and in vitro transcriptome data. We named the set of prediction models “regression analysis–based inductive DNA microarray (RAID),” which inductively analyzes the mode of action and the key event in adverse effects with reference to the redundant arrays of inexpensive disks (also represented as RAID), a data storage virtualization technology that combines multiple physical disk drive components with the purpose of data redundancy. As RAID (storage technology) complements data based on information of multiple components, we hope that RAID (our microarray) will complement the relationships between multiple media (in vivo gene expression, in vitro gene expression, and chemical structure). Our RAID system achieved a quantitative in vitro to in vivo extrapolation (QIVIVE) by the integration of a structure-based approach (QSAR) with transcriptomic data. Whereas general “Q”IVIVE studies predict dose (or concentration) quantitatively in toxicological or toxicokinetic effects, our “Q”IVIVE predicts in vivo gene expression values quantitatively. Finally, the substance similarities were analyzed by principal component analysis (PCA), which proved useful in understanding the features of toxic substances based on their gene expression profile (Watanabe et al., 2012), using RAID (the virtual microarray) data, in vivo data, in vitro data, and chemical structure data to validate the usefulness of read-across.

The transcriptome data signatures derived from the RAID (the virtual microarray) system were in good agreement with those of in vivo data, and the technology provided an understanding of the features of hepatotoxic substances based on the toxicological mechanism interpretation. The mechanism of action of the two characteristic toxic substances separated using PCA analysis was shown to be achieved through Cyp1a induction via AHR and Cyp4a induction via PPARa (pathway and gene ontology analysis). The AHR-induced drugs raise safety concerns during developmental periods (Qin et al., 2019), and PPARa-induced drug toxicity requires species differentiation considerations (Ito et al., 2006). Therefore, predicting the involvement of these nuclear receptors and induction of metabolic enzymes is critical for understanding the molecular initiating events and the key events associated with adverse outcome pathway. RAID enables the prediction of gene expression levels, thus exhibiting properties required for the next-generation risk assessment methods.

The first substance group (TAA, MP, and HCB), representing toxic substances commonly differentiated from non-toxic substances using PCA on in vivo and RAID data, has been reported to have carcinogenicity with metabolic activation (Uehara et al., 2008; Hajovsky et al., 2012; US HSS, 2015). Furthermore, these substances have been shown to activate xenobiotic-related receptors, such as AHR inducing Cyp1a (Ushel et al., 2002; Yamashita et al., 2014; Clara et al., 2015). Moreover, in vivo transcriptome data in this study showed that TAA, MP, and HCB induce Cyp1a activation. AHR is known for mediating the toxicity and tumor promoting properties despite the mechanism through which AHR activates carcinogenesis needing to be elucidated (Safe et al., 2013; Murray et al., 2014).

The second substance group (WY, FFB, BBr, and GFZ) includes fibrates which are recognized as PPARa agonists (Schoonjans et al., 1996), implying that induction of Cyp4a via PPARa and perturbation of lipid-related genes are involved as a series of key events. Although another fibrate included in training data—clofibrate (CFB)—was classified as a non-toxic substance according to no serum chemistry findings from a previous study, CFB was shown to act as a PPARa agonist inducing peroxisomal proliferation on hepatocytes (Low et al., 2011) and was plotted around the second group in the PCA. Sustained activation of PPARa signaling and induction of enzymes, such as CYP4A, to increased fatty acid oxidation contributes to sustained oxidative stress in the liver. These changes lead to liver cell damage as hypertrophy and proliferation which contribute to the development of hepatocellular carcinoma (Parimal et al., 2013).

From the perspective of capturing individual gene responses, RAID could detect gene expressions related to major drug metabolism responses in in vivo more broadly (more common principal component–related gene number; Figure 4) and quantitatively (lower RMSE values; Figure 6) than in vitro. The 36 genes that were commonly related to the principal components of in vivo and RAID data contained genes that were known to be involved in drug metabolism and hepatotoxicity. In addition to the genes described above (Cyp1a and Cyp4a), Acot1 acts as an auxiliary enzyme in the oxidation process of various lipids in peroxisomes (Hunt et al., 2012). Furthermore, Vnn1 is expressed by the centrilobular hepatocytes and is involved in lipid and xenobiotic metabolism (Bartucci et al., 2019), whereas Pex11a (peroxisomal biogenesis factor 11 alpha) is involved in peroxisome maintenance and proliferation associated with dyslipidemia (Chen et al., 2018). All of these genes are known as PPARa target genes (Rakhshandehroo et al., 2010; Lake et al., 2016). Thus, these features indicate that RAID can predict possible toxicity by taking into account a broader range of mechanisms than the range of in vitro data. Indeed, the in vivo changes detected using the in vitro data were limited (Figure 4), and the PCA showed that most of the differentially expressed genes were associated with irrelevant nonphysiological conditions. Thus, the IVIVE effect combining the QSAR technique and in vitro data would allow for more precise predictions through de-noising these types of in vitro specific biological responses.

In vitro data contribute to accurate gene expression predictions that could not be achieved with QSAR alone (Figure 2D). In vitro data contributed to the prediction of the mechanism shown in Figure 5. The biological mechanisms related to metabolic processes were consistent with the key mechanisms of characteristic hepatotoxic substances described above, which indicates that in vitro data contributes to the precise predictions obtained using RAID. In addition, whether in vitro responses were observed in the suggested mode of action predicted by the RAID system or not is an important point in terms of weight of evidence. This study provides valuable evidence supporting that transcriptome data should be considered in light of previous reports indicating that in vitro data does not necessarily reflect in vivo conditions (Tamura et al., 2006; Sutherland et al., 2016). Simultaneously, in vitro studies focusing on a specific mechanism should consider the external validity of their findings and whether the findings reflect in vivo situations.

Evaluating the read-across performance using external substances, such as 3,4,5,3′,4′-pentachlorobiphenyl, 2,2′,4,4′-tetrachlorobiphenyl (a type of polychlorinated biphenyl) and pentachlorodibenzofuran (dioxin-like compounds) (Figure 7A), which are known as IARC group 1 carcinogens and Cyp1a1 inducers (EPA,U S, 1996; Walker et al., 2005; National; Toxicology Program, 2006); these were separated as toxic substances. Additionally, benzo(a)pyrene, 3-methylcholanthrene, and 9,10-dimethyl-1,2-benzanthracene plotted apart from the origin of coordinates (PC1 = 0 and PC2 = 0), and are polycyclic aromatic hydrocarbons inducing Cyp1a1 (Moorthy et al., 2007; Pushparajah et al., 2008). Non-carcinogenic chemical substances, such as food components or preservatives, were positioned near the origin, second quadrant or third quadrant, indicating low risk. Furthermore, substances interacting with Cyp4a (Figure 7B), such as ciprofibrate, nafenopin, clofenapate, clofibric acid, and di(2-ethylhexyl) phthalate, which plotted in the area of the 2nd substance group (PC1 > 0), are also known as PPARa agonists (Bocos et al., 1995; Roberts et al., 2002; Yadetie et al., 2003; Currie et al., 2005; Pyper et al., 2010). Chemicals that were not characterized by the PC1 component (PC1 < 0) are not hyperlipidemia drugs. These results suggest that the RAID system effectively classifies substances based on their mode of action as well as the strength of their toxicity, and ultimately contributes to precise read-across. Thus, the RAID system provides a new method for read-across in line with IATA that should be called “a virtual functional read-across”. Here, we showed that substances without high structural similarities might have similar toxicological properties, and our new approach interpreted the shared mechanism of action. This means that RAID considers the qualitative and quantitative similarities of biological responses, which was one of the major issues of QSAR-based read-across. The structural similarities of TAA, MP, and HCB observed using correlation coefficients of the chemical descriptor used for the predictive model, and the maximum common substructure (MCS) similarities with the Tanimoto coefficient, were less than 0.5; however, the homology of RAID and in vivo data was as high as a 0.8 Pearson’s correlation coefficient. Furthermore, achieving such an accurate read-across without using in vitro data will provide a new perspective on the structural information-based predictions.

PCA analysis was used to understand the features of substances to predict the modes of action and identify biologically similar substances for read-across in this study. The examples of applications of RAID for read-across described above were compared to other methods (Table 5). The RAID system could enhance read-across reliability by estimating toxicity including modes of action, while this was difficult by other methods (e.g., QSAR or read-across using chemical structure data) (Supplementary Figure S1). On the other hand, focusing on certain specific toxicities, discriminant analysis, classifier model, or biomarker analysis might improve the separation of toxic substances. Indeed, as shown in Supplementary Table S1, when the RAID system was applied to discriminating hepatotoxicity, PLS-DA using RAID data showed a good predictive performance, indicating usefulness for the specified toxicity prediction. Thus, the use of RAID data instead of experimental transcriptome data would achieve previously reported biomarker-based classification without using animals. For example, Liu et al. (2017) indicated that certain genes were associated with hepatocellular hypertrophy and hepatocarcinogenesis, as well as markers such as Cyp1a1, Acot1, Stac3 (SH3 and cysteine rich domain 3), and Hdc (histidine decarboxylase), which were correctly evaluated in the present study to characterize hepatotoxic substances. Similarly, the constructed RAID system could be applied to previous studies to predict carcinogenicity or estimate transcriptional benchmark dose by toxicogenomics analysis of short term in vivo studies (Ellinger-ziegelbauer et al., 2008; Thomas et al., 2013b; Matsumoto et al., 2014; Kawamoto et al., 2017).

One important issue that should be considered in toxicological evaluation using the RAID system is consideration of species differences. The RAID system provides mechanistic insights on repeated-dose toxicity in animal models; however, since some species differences have been observed, the suggested mode of action and the corresponding molecules need to be confirmed by toxicologists. The interspecies extrapolations, such as rat-to-human extrapolations, could be achieved by further experiments to construct new RAID systems with these different species’ transcriptome data. In addition, RAID data of substances that were separated as toxic substances in PCA (e.g., TAA, MP, HCB, WY, FFB, BBr, and GFZ) showed high similarity to in vivo data (Supplementary Table S2). Since regression analysis requires certain levels of standard deviation of training data, the RAID accuracy for substances may be related to the number of substances with similar modes of action. Thus, database expansion for several substance groups with minor modes of action would contribute to further improving the accuracy and applicability domain. In addition, evaluation of RAID usefulness for various toxicities is required.

The present approach integrates QSAR and IVIVE and will contribute to other areas of research, such as drug repositioning, which has recently attracted attention toward pharmaceuticals that are available on the market and might be repurposed for new diseases (Jourdan et al., 2020). However, the previously proposed methodologies (Iwata et al., 2018; Lippmann et al., 2018; Zhu et al., 2020; He et al., 2021) have room for improving the IVIVE aspect of in vivo predictions. Thus, our system provides an alternative to screening candidate drugs and exploring new biologically similar drugs at a low cost.

In conclusion, we developed a virtual DNA microarray system that quantitatively predicts in vivo gene expression profiles based on the chemical structure and/or in vitro transcriptome data. Estimated transcriptomes are considered scientifically relevant from PCA data interpretation as well as pathway and GO analysis. Based on its external validation, our system works as an alternative test for repeated dose toxicity tests with toxicogenomics analysis enabling IVIVE and mechanism estimation. Although our technology might have limited applicability domain due to the small data size of chemical substances and their characteristics (using hepatotoxic substances), the concept of the virtual microarray analysis contributes to the 3Rs (reduction, refinement, and replacement) and might benefit much future animal testing.