A new high-accuracy QSAR model based on ultra-curated data for predicting thyroid receptor-related endocrine activity

1 Introduction

In 2002, the World Health Organisation (WHO) defined an Endocrine Disruptor as “a substance or mixture that alters function(s) of the endocrine system and consequently causes adverse health effects in an intact organism, or its progeny, or (sub) population” (WHO and IPCS, 2002). Based on this definition, ECHA/EFSA defines in 2018 when a substance is considered as having an endocrine disrupting mode of action, i.e., there is a biologically plausible link between the adverse effect and the endocrine activity (ECHA/EFSA, 2018). The endocrine modulation can occur at multiple levels, including hormone synthesis, transport, receptor binding, and metabolism (La Merrill et al., 2020). A striking example is diethylstilbestrol, a synthetic non-steroidal oestrogen once prescribed for miscarriage prevention, which was later shown to induce foetal malformations through oestrogen receptor activation (Robotti, 2021). In addition to receptor-mediated effects, diethylstilbestrol also induces epigenetic modifications, leading to transgenerational outcomes such as carcinogenicity and teratogenicity (Klip et al., 2002). This illustrates the diversity of endocrine-disrupting mechanisms and has fuelled scientific and regulatory debate regarding their identification and classification (Autrup et al., 2020; Brescia, 2020). Regulatory agencies and authorities, however, have emphasised receptor-mediated mechanisms as a primary basis for classification. In 2023, new hazard classes for endocrine disruptors have been introduced in Classification, Labelling and Packaging (CLP) Regulation. In 2024 the European Chemicals Agency (ECHA) incorporated endocrine disruption potential into its CLP guidelines, further reinforcing the need for reliable methods to identify such chemicals (ECHA, 2024). In response, a variety of in vitro assays have been developed to detect agonist or antagonist activity, including hormonal pathways of high concern such as the thyroid axis (OECD, 2024a).

Currently, no new GHS hazard class for ED is available. In 2025, the OECD ad hoc group of various experts reviewed existing GHS hazard classes (e.g., carcinogenicity, reproductive toxicity, specific target organ toxicity, hazardous to the aquatic environment) that may include endocrine disrupters. The group agreed that there is a gap in GHS hazard classes regarding the possibility to identify EDs as such, according to the WHO/IPCS definition. In particular, inclusion of mechanistic data is not mentioned as part of the identification of most hazards (mutagenicity is an exception), which would be a gap in the current ability to identify a chemical as an ED under the GHS (OECD, 2025a). Moreover, using the IPCS/WHO definition of EDs, no stand-alone method can be used to identify a substance as an ED, and rather, methods must be used in combination. In this regard, in silico methods such as those described here, have their place in tier one alongside other studies in the battery to be implemented to identify EDs chemicals (ECHA/EFSA, 2018).

The thyroid gland produces the hormones triiodothyronine (T₃) and thyroxine (T₄), which regulate numerous physiological processes across multiple organ systems. Thyroid hormones control metabolic functions through effects on the liver, including regulation of thermogenesis, glucose production and uptake, cholesterol homeostasis, and fatty acid oxidation. They also modulate cardiovascular function by influencing heart rate and contractility, and contribute to musculoskeletal regulation via bone cell differentiation, muscle contraction, and creatinine metabolism. Most critically, thyroid hormones are indispensable for brain development, where they direct neuronal maturation, axonal myelination, and synaptic connectivity (Norman and Henry, 2015).

Consequently, disruption of thyroid signalling can have profound and lasting consequences, particularly during foetal and neonatal development (Feldt-Rasmussen and Mathiesen, 2011; UNEP/WHO, 2012). While hyper- and hypothyroidism are commonly managed as clinical conditions arising from autoimmune diseases (e.g., Hashimoto’s disease, Graves’ disease) or iodine imbalance, growing evidence implicates exposure to environmental chemicals as a cause of thyroid disruption (Fernandez et al., 2018; Taha et al., 2020; Marques et al., 2022). For example, certain benzophenones have been associated with thyroid-related effects, although their precise mechanisms of action remain poorly defined (Axelstad et al., 2013; Zhang et al., 2018; Mustieles et al., 2023) and may not be related to an interaction with the thyroid hormone receptors.

Historically, thyroid disruption was assessed through reproductive and repeated-dose studies in rodents (OECD, 2025b). However, in line with the 3R principles (Replacement, Reduction, Refinement), mechanistic in vitro approaches have gained importance. In 2024, the OECD reviewed eleven in vitro methods for thyroid disruption assessment, and large-scale screening programs such as the U.S. Tox21 initiative employed assays like the GH3 cell-based qHTS luciferase system (OECD, 2024a). While these assays improve mechanistic understanding and reduce animal use, they remain resource-intensive and insufficient to screen the vast number of chemicals in commerce or predict the activity of new ones. Moreover, they cannot fully replace in vivo testing.

To provide upstream strategic direction to experimental strategies, in silico approaches such as quantitative structure-activity relationship (QSAR) modelling is increasingly promoted within integrated approaches to testing and assessment (IATA). QSAR models offer the potential for high-throughput, cost-effective prediction of endocrine activity and when built and used appropriately, can indicate and limit future testing only to areas of potential concern. However, most existing models for thyroid receptor-related activity are limited in accuracy, rely on poorly curated datasets, and fall short of the standards required for regulatory use outlined in the OECD QSAR Assessment Framework (Sellami et al., 2022; OECD, 2024b; Vergauwen et al., 2024; Evangelista and Papa, 2025). These shortcomings result in frequent false positives and undermine regulatory acceptance.

In this context, the aim of this work was to develop a new high-accuracy QSAR model, further to a call for proposals in the NC3Rs “CRACK-IT Challenge” research program (NC3Rs, 2026), specifically designed to predict thyroid receptor-related endocrine activity. The model described here integrates curated datasets by toxicologists, offering improved predictive performance and enhanced regulatory applicability.

2 Materials and methods

2.1 Dataset construction and analysis

2.1.1 Dataset sources

The main source of data used to train and validate the model was the set of multi-concentration activity curves resulting from the GH3 Cell-Based qHTS luciferase assay to identify agonist or antagonist of the TR signalling pathway (Freitas et al., 2014). The dataset was primarily sourced from public invitroDB v3.3 (Williams et al., 2017; EPA’s National Center for Computational Toxicology, 2020).

(named “Tox21” in the following sections), and the following assays were used to provide data for this model.

• Tox21_TR_LUC_GH3_Agonist
• Tox21_TR_LUC_GH3_Antagonist
• Tox21_TR_LUC_GH3_Agonist
• Tox21_TR_LUC_GH3_Antagonist_viability

Since this dataset was found to contain only a small proportion of ‘actives’ (i.e., agonists or antagonists), it was complemented by the addition of positive in vitro data from the BindingDB database (Liu et al., 2025), where data come from manual curation of articles in selected scientific journals and US patents by the BindingDB team (named “BindingDB” in the following sections). Quantitative values of the activity and/or the affinity from in vitro experiments involving compounds and target proteins were extracted. The data in BindingDB do not necessarily come from the same assay type as those in the Tox21 dataset but data selected from this database are based on similar cell-based reporter assays providing EC50 values or IC50 values against a reference T3 activation of TR in all cases.

2.1.2 Curation method

Available SMILES were standardised by removal of salts, stereochemistry was not considered, and a selection of a single converging tautomeric state was carried out. Hence, using translation to InChiKey from the canonical SMILES of mono-constituents, alternative SMILES codes in the input data representing the same principal compounds were removed, resulting in 8,282 standardised SMILES allowing description of 7,155 mono-constituents and 1,816 multi-constituents (after salt removal). The subsequent input data analysis and model development used only records of mono-constituent substances.

The Tox21 actives were further validated by experts through the comparison of activity curves versus cytotoxicity curves. The rational of this comparison arose from detection of extreme cellular stress, ascertained by a decrease in intracellular protease level and cell membrane integrity (Paul-Friedman et al., 2019). In such cases the specific activity on THRs monitored by the ATP-dependent luciferase becomes irrelevant, since a dying cell can also be measured to have decreased ATP levels (Freitas et al., 2014). Hence, the decision to change the classification label from the previous ‘Active’ to curated ‘Inactive’ in such conditions is supported by the assumption that the monitored receptor transactivation signal depends on healthy cell metabolism.

We empirically calibrated the computational formula for the delay between activity and cytotoxicity signals to align with expert assessment (Supplementary Figure 6). A scalable threshold was established to monitor the overlap of dose-response curves of both endpoints, thereby hinting to contrast genuine antagonism from cytotoxicity artifacts. This distinction is critical, as cell death occurring within the antagonist concentration range can suppress transactivation signals, mimicking specific inhibition (Freitas et al., 2014).

The cytotoxic activity of tested compounds was assessed throughout the Tox21 program by monitoring cell protease levels and cell membrane integrity using CellTiter-Fluor technology in the same wells as the antagonist assay, prior to adding the reagent for luciferase activity measurement (National Center for Biotechnology Information, 2014; Paul-Friedman et al., 2019). However, this measurement is only a partial descriptor of cell viability. It was assumed that any specific activity of the test compounds would be relevant only at concentrations below the cytotoxicity limit observed in the satellite assay (Tox21_TR_LUC_GH3_Antagonist_viability), and it was also hypothesised that the cell could start to be non-viable at concentrations slightly lower than the measured level in the viability assay by monitoring protease levels, thus apparent endocrine activity starting at concentrations below observed cytotoxicity could be symptomatic of the loss of viability rather than a true specific interaction with THRs. An alternative criterion not used in this study was proposed by Freitas et al. (2014), where the authors applied a safety factor to limit the risk that observed antagonist activity is due to cytotoxicity: “[antagonist]IC50 values should be at least three times lower than the IC50 in the viability assay”.

The raw data quality assessment and filtering were applied by US EPA as defined in the database release documentation (EPA’s National Center for Computational Toxicology, 2020). Nevertheless, further curation of the input dataset was semi-automated to reproduce the process applied by human experts manually examining the quality of the Tox-21 assay data for a given tested substance. At this stage, various parameters of the models fitted to the raw signal in the original database were compared between activity assays and satellite viability assays (datapoints from the experimentally collected data were not individually considered). Substances remaining with the class label ‘Active’ (positive) after the automated curation were further validated by qualified experts examining the ensemble of experimental datapoints within an assay record. Subsequent manual curation of fitted models, preserving the positivity class label, considered additional external knowledge criteria for known biological or optical activity of substances interfering with the transactivation assay. For instance, specific luciferase inhibitors have been identified (Yonchev and Bajorath, 2020) that would lead to apparent antagonist signals in the assays considered in this article. A small sample of substances which were labelled as “inactive” after the automated curation were also reviewed by experts, confirming that all of them were also interpreted as “inactive” by experts.

In our curation protocol, a fitted model is considered representative of genuinely positive in vitro activity if it fulfils the following criteria:

• As a minimum it at least reaches the activity threshold for each endpoint defined by the US EPA in the design of the database. These thresholds were strictly accepted without modification.
• Its shape is: (i) a fully observed ascending sigmoidal curve (i.e., one can see the minimum activity plateau at ca. 0% and an ascending part) (ii) a partially observed sigmoidal curve (i.e., the minimum cannot be seen as activity appears to start below the lowest tested concentration); that can even appear as a plateau when the maximum activity was reached below the lowest tested concentration; (iii) a fully observed ‘gain-loss’ curve; (iv) a partially observed ‘gain-loss’ curve that appears to start with a plateau already reached before the lowest tested concentration and beginning to decrease after the third lowest experimentally tested concentration. “Gain-loss” models partially represented over the tested concentration range and exhibiting only the descending sigmoid section from the lowest tested concentration were not accepted as “Active”, because these fitted characteristics easily arise from optimal curve-fitting of noisy data.
• In the specific case of ‘gain-loss’ models the plateau must be represented by more than one tested concentration above the activity threshold.

If activity did not start at a significantly lower concentration than cytotoxicity, it was reassigned as “not active up to the cytotoxicity limit” or more simply “inactive”. Confirmed agonists or antagonists below the cytotoxicity limit were labelled as “active”. Some examples of such situations that were interpreted as “not active up to the cytotoxicity limit” are shown in Figure 1. The calculation of the “delay” between apparent activity and cytotoxicity and the defined threshold to consider that the curves are sufficiently close to each other is described in the supplementary information.

Subsequently, substances for which the activity raised above noise level only at the highest concentration tested (between 50 and 100 µM depending on the samples) were also not accepted as positives, since cytotoxicity had a high chance of causing this apparent activity at such high concentrations (Supplementary Table S7: ∼30% of cytotoxicity labels assigned to substances at concentrations <15µM, 60% at concentrations <50µM, and ∼10% at concentrations >100 µM). The cytotoxicity assay occasionally can be extrapolated as positive at concentrations slightly higher than the maximum tested in the assay (e.g., see Figure 1C). In such cases, if the agonist/antagonist assay was positive within the assay’s concentration range and if higher concentrations were tested, then the observed activity would likely persist. However, as there are increasing signs of cytotoxicity, we nevertheless considered these cases as not active up to the cytotoxicity limit, as discussed with the curation rules that we established and applied to all assay results.

• Class labels of fitted models of activity were set to Negative if the activity class was assigned for the highest tested concentration and if in parallel at this concentration there is a cytotoxicity signal above the background noise (but below the cytotoxicity limit). This rule was used as a shorthand for extrapolation of the cytotoxicity signal to cytotoxicity limit and its fully defined comparison with the activity signal using our special points ratio procedure (see Supplementary Figure 6).

2.1.3 Data augmentation and dataset mergers

To be consistent with the curated subset of Tox21 in vitro data, data extracted from Binding DB was filtered down to records expressed as EC50 and IC50, because these data were obtained from cell-based in vitro assays, thus compatible with the Tox21 curated subset. The original data sources were not further investigated for veracity, consistency or data quality, because of the high number of data and because the peer-review of the BindingDB and primary sources were considered trustful. We included records described by EC50 or IC50 lower than 100 nM in our modelling dataset. According to descriptors collected from the curated Tox21 subset, we estimated a risk of 2% that a compound is cytotoxic when it is active in vitro at <100 nM. The threshold value was deduced from the Tox21 GH3 cell-based cytotoxicity assays where it was identified that cytotoxicity frequently started to occur at or above concentration of ca. 1 μM (see Supplementary Table S7, ∼10% of cytotoxicity labels given for concentrations <3 µM). The 100 nM threshold was obtained by further applying a safety factor of 10.

Finally, after data augmentation of the curated Tox21 dataset with the BindingDB database values, all compounds (with the same InChI) tested more than once and presenting divergent activity class labels (active and inactive) were filtered out from the curated dataset. Finally, after removal of salt counter-ions from all entries in the dataset, the following additional filters were applied: keep compounds represented by a single major constituent (mono-constituent like), and compounds containing atoms outside of the following list: Li, Be, C, N, O, F, Na, Mg, P, S, Cl, K, Ca, Se, Br, I.

As the BindingDB data was not systematically labelled as either agonist or antagonist, Tox21 data in the form of “agonist”, “antagonist” and “inactive” were relabelled as only two classes: “active” (grouping “agonists” and “antagonists”) and “inactive”. Several limitations also precluded the separation of the molecules into dedicate datasets for the agonist and the antagonist mechanisms of action. First, the Tox21 derived data inherently merges TRα and TRβ signals. Without the ability to distinguish between these biological targets, it is impossible to exclude mixed modes of action, wherein a molecule might exhibit activating influence on one receptor subtype and inhibiting influence on the other. Furthermore, partial or weak activators often present as competitive inhibitors when monitored in antagonist in vitro assays (such as TOX21_TR_Antagonist), making strict assignment to specific mechanistic training sets highly uncertain.

The final dataset is composed of 6,340 substances, of which 6,049 inactive substances and 291 active substances (containing 12 positives in agonist assay and 36 positives in antagonist assay from the Tox21 TR dataset, and 243 additional TR active substances from BindingDB). The dataset was highly imbalanced, with 95.4% of the data representing the inactive class and only 4.6% representing the active class. Given the small number of identified agonists and antagonists, it proved not feasible to train a model with these as separate classes, thus also supporting the choice of merging these two classes as “active”, and the choice of augmenting the dataset with BindingDB data despite the subsequent uncertainty in data quality.

2.2 Data-splitting procedure

All dataset partitions were made with in-house Python algorithms employing RDKit (Landrum et al., 2023) functionalities including the popular algorithms Butina clustering (Butina, 1999) and Kennard-Stone sampling (Kennard and Stone, 1969). Pairwise molecular similarities were computed using KREATiS′ proprietary hybrid 1D/2D fingerprint that is a modified version of PubChem fingerprint inspired by the implementation in PyBioMed (Dong et al., 2018) T c Similarity was quantified by the continuous Tanimoto coefficient (T c) (Tanimoto, 1958).

2.2.1 Split Level 1 - Hierarchical butina clustering

This methodology is schematised in Figure 2.

• Given the pairwise similarity T c, we assume working in the metric space X , T c where, a pairwise distance is defined as: Entries are merged in the same cluster when their pairwise distance is below an arbitrary defined threshold ε. The threshold at the first iteration of the Butina clustering was set toand the maximum distance at which Butina clustering would be performed was set to
• At initialisation step t = 0, all entries are part of their own cluster, e.g., all entries are singletons:
• The first Butina clustering step (t): We distinguish the clusters by their cardinality. Multimember clustersand singleton clusters
• If there were no singleton clusters in D t then go to Step 9, or else update the iteration counter
• In order to accomplish efficient iterative clustering, we needed to verify the starting point of each subsequent Butina clustering run. We started by selecting central entries (s e e d s) from each clusterand kept track of the remaining entries within each multimember cluster Naturally, all singletons are seeds in their own clusters. In addition, purely for algorithmic purposes, we divided any cluster into singletons if it was composed of pairwise equidistant sets because it is not pertinent to select a cluster seed in such conditions. We consider the family U α ∈ D t of finite equilateral subsets, where each U α satisfied d T x i , x j = δ α with ∀ i ≠ j. Hence, we defined the merged set: Ultimately, the total set of seeds for the next Butina clustering iteration became:
• We performed the next Butina clustering iteration:
• We propagated the new cluster identifiers to the previously masked entries (D m a s k t) through their association with s e e d i ∈ B t.
• If ε < ε l o w then we went to Step 4.
• Result: all entries were either assigned to a cluster of similar molecules or remained a singleton (not similar to any other molecule up to cutoff ε l o w)

Contrarily to the classic Butina algorithm applied with a single ε l o w, the iterative scheme splits this toy-dataset into two clusters and a singleton. While, in the classic method there would be one large cluster whose seed would be the element highlighted in blue in Cluster 1, and includes the four closest elements, two from the depicted Cluster one and two from Cluster 2. Then, a singleton and two smaller clusters would be composed of two and three elements.

2.2.2 Split Level 2 - Stratification by experimental activity

According to reviewed experimental THRBα/β activity, each cluster C i ∈ D t (where i enumerate cluster identifiers) was split by positive activity (including both agonists and antagonists) class label P, and negative activity class label N:

P i ← C i \ N i

2.2.3 Split Level 3 - Isolation of the external validation set

The rules listed in Table 1 were applied to subsets P i and N i independently:

TABLE 1. Data splitting strategy for external validation set depending on the size of clusters.

After the preprocessing in Splitting Level 1 and Splitting Level 2, the desired exact splitting of the curated dataset in 20%:80% subsets (classical Train/Test splitting scheme) was not straight forward. An exact 20%:80% ratio is only possible if all P i and N i clusters fulfil simultaneously the conditions:

c a r d P i ≡ 0 mod 5

c a r d N i ≡ 0 mod 5

Since these conditions were not exactly met, the global outcome of the Train/Test splitting workflow was.

• Prospective external validation set ≈20% (immediately locked)
• Training + validation pool ≈80%

2.2.4 Split Level 4 - Creation of final training and internal validation sets

After isolation of the Train subset (Split Level 3), the procedure described in Split Level 3 was applied on the Train subset to isolate a validation pool, aiming at the splitting 80%:20%. Eventually the exact desired splitting ratio was not possible to achieve due to the fact that some P i and N i clusters have less than ten members.

Final global partition (typical values):

• Final training set 64% of original dataset
• Internal validation set 16% of original dataset
• Prospective external validation set 20% (fully independent, never used in Levels 1–4 operations on the other subsets)

2.3 Calculation of molecular fingerprints

Several molecular fingerprints were computed for each molecule using a combination of cheminformatics libraries. Substructure-based fingerprints (e.g., PubChem, MACCS keys) were generated using PyBioMed, which internally relies on RDKit for molecular handling. Circular fingerprints such as ECFP6 and FCFP6 were computed directly with RDKit, using the Morgan algorithm and appropriate atom invariant definitions (Landrum et al., 2025). Additional fingerprints commonly used in QSAR modelling, such as complexity-related or topological descriptors, were also included depending on the model requirements.

A customised version of the PubChem fingerprint was developed. Starting from the implementation in PyBioMed, the features-counting scheme was adjusted. This fingerprint was used in the dataset splitting framework.

Each fingerprint type captures different structural or functional aspects of the molecule, providing complementary representations for machine learning.

2.4 Model building methods

2.4.1 SVM model

Support Vector Machines (SVMs) are supervised learning techniques widely used for classification tasks. First, Introduced by Vapnik and colleagues (Cortes and Vapnik, 1995), SVMs aim to identify an optimal decision boundary, a hyperplane, which separates the data into distinct classes (Burges, 1998). However, real-world datasets are often not linearly separable in their original feature space. To address this, SVMs use a kernel, a mathematical function that measures similarity between points and implicitly maps the data into a higher-dimensional feature space in which a linear separation becomes possible. In the present model, the following hyperparameters were fine-tuned:

• C: This parameter controls how much the model allows for errors during training. By setting C to a lower value, the model focuses on creating a wider margin between the classes, even if it means misclassifying some points. This helps avoid overfitting and makes the model more general.
• Gamma: Gamma is particularly critical for nonlinear kernels, as it governs the flexibility and smoothness of the decision surface.
• Kernel: Specifies the type of transformation applied to the data to enable the model to capture linear or nonlinear relationships depending on the problem.

2.4.2 Optimisation of model hyperparameters

The SVM hyperparameters, namely, the penalty parameter C and the kernel parameter γ, were optimised using grid search combined with five-fold cross-validation on the training set, where the five folds were generated using the splitting method provided by KREATiS, in order to identify the configuration offering the best predictive performance. Once the optimal hyperparameter configuration was identified, the model was evaluated on the internal validation set, and the best model was selected based on its performance on this validation set.

The optimisation was performed using the GridSearchCV function from Scikit-learn (Pedregosa et al., 2011), a widely used Python library for machine learning. This approach allows for the systematic exploration of a predefined parameter grid while evaluating each configuration via cross-validation, thus ensuring improved model performance while preserving its generalisation ability.

The optimisation process for a binary classification (active or inactive) also evaluated different types of molecular fingerprints, including substructure-based fingerprints (e.g., PubChem and MACCS keys) and circular fingerprints (e.g., Feature-Class Fingerprints FCFP6 and Extended-Connectivity Fingerprints ECFP6). The best predictive performance was obtained using circular fingerprints (FCFP6) as input features with C = 0.1 and kernel = “linear”. During the grid search, the F1-score was used as the evaluation metric to select the optimal hyperparameter configuration, ensuring a balanced assessment of both precision and recall.

2.5 Evaluation of model performance

To assess the robustness and unbiased predictive performance of the models, five-fold cross-validation and external validation were conducted. The predictive performance and robustness of the model were evaluated using a range of metrics, including classification accuracy (CA), sensitivity (SE), specificity (SP), precision, F1 score, and balanced accuracy (BA), as shown in Equations 1–6. In this study, an active compound misclassified as non-active is defined as a false negative (FN), whereas a non-active compound misclassified as active is defined as a false positive (FP). Correctly classified active and non-active compounds are referred to as true positives (TP) and true negatives (TN), respectively. Sensitivity and specificity correspond, respectively, to the proportion of correctly identified active compounds among all active compounds and the proportion of correctly identified non-active compounds among all non-active compounds. Classification accuracy represents the overall proportion of correctly classified compounds.

The formulas used to calculate these metrics are given below:

S E = T P / T P + F N (1)

S P = T N / T N + F P (2)

C A = T P + T N / T P + T N + F P + F N (3)

P r e c i s i o n = T P / T P + F P (4)

B A = S E + S P / 2 (5)

F 1 = 2 * T P / 2 * T P + F P + F N (6)

Additionally, F1 score and precision were calculated separately for the 1 (active) and 0 (inactive) to provide a more detailed evaluation of model performance in imbalanced datasets.

In addition, the receiver operating characteristic (ROC) curve and the precision-recall (PR) curve were used to assess the performance of the binary classification model. An area under the ROC curve (AUC-ROC) or PR curve (AUC-PR) of one indicates a perfect classifier, whereas a value of 0.5 for AUC-ROC (or a low AUC-PR) indicates poor discriminative ability.

2.6 The combined applicability domain and reliability assessment model

In this present work, “the combined applicability domain and reliability assessment model” refers to a combination of the applicability domain (AD) and the reliability assessment of the SVM model’s predictions.

2.6.1 Assessment method for applicability domain

The applicability domain defines the region of chemical space and response space in which a model can make reliable predictions. Two complementary approaches were used to assess the domain: where:

• Assessment of structural domain using Atom-Centred Fragmentation: For each molecule in the training set, atom-centred molecular fragments were generated and stored using RDKit. Each fragment included all atoms within a topological radius of 1, capturing the local chemical environment. To assess whether a new compound falls within the structural domain of the training set, its fragments are generated in the same way and compared against the stored training fragments. If all fragments are present in the training set, the compound is considered within the structural domain of the model.
• Assessment of the descriptors domain using the k-nearest neighbours: To define the applicability domain based on molecular descriptors, the k-nearest neighbours (kNN) method first was calculated, for each compound of the training set, the average distance to its k nearest neighbours. From these average distances, a reference value was then determined using the interquartile range of the distances, as proposed by Sahigara et al. (2013):

d ∼ k = Q 3 d ¯ k + 1.5 × Q 3 d ¯ k − Q 1 d ¯ k

• d ¯ k is the average distance to the k nearest neighbours.
• Q1 and Q3 are the 1st and 3rd quartiles, respectively, of the distances d ¯ k.
• d ∼ k is a reference value.

Next, the ordered distances of each i-th compound of the training set from all other n-1 compound of the training set were compared with the reference value. If the distance of the i-th sample from its j-th training neighbour (where1 ≤ j ≤ n-1 was less than or equal to the reference value, it was retained otherwise, it was discarded. Finally, each i-th compound of the training set was associated with a threshold t_i which defined the width of its neighbourhood as:

K i : D i j ≤ ~ d k , ∀ j = 1 , … , n ‐ 1 t i = 1 K i ∑ j = 1 K i D i j

For a test compound, its distance to all compounds of the training set was calculated and compared with the thresholds t_i. If the distance to at least one compound of the training set was less than or equal to its local thresholds, the compound was considered inside the domain otherwise it is outside the domain.

More formally, given the training set (TS), for each test compound j, the decision rule was:

n j = ∑ i ∈ T S b 1 D i j ≤ t i , j ∈ A D ⇔ n j ≥ 1

• When n_j = 0 (number of training set neighbours that are sufficiently close according to the threshold t_i) the compound was considered outside the descriptor domain.

In conclusion, a compound was considered outside the applicability domain if it is identified as out-of-domain by either the descriptor-based domain or the atom-centred molecular fragment approach.

2.6.2 Assessment method for the reliability of prediction

The process was based on a binary classification model that categorises compounds prediction into two classes: “Reliable” and “Unreliable”. Analogues of the test compound were identified using a “k-Nearest Neighbours” approach based on the Tanimoto similarity index using a molecular fingerprint (Figure 3). Only compounds exceeding a predefined similarity threshold were retained. If no analogue was found, that is, if k = 0, the compound prediction was immediately classified as “Unreliable” (STEP 1).

When analogues were available, i.e., when k ≥ 1, two criteria were evaluated. The first criterion, accuracy, measured the agreement between the predicted and experimental values of the analogues. Accuracy was equal to one only if the predicted value of each analogue exactly matched its experimental value; otherwise, accuracy ≠ 1. The second criterion, concordance, reflected the consistency between the predicted value of the test compound and the experimental values of its analogues. Concordance was equal to one only if the predicted value of the test compound fully agreed with all experimental values of its analogues; otherwise, concordance ≠ 1. A compound prediction is classified as “Reliable only if both accuracy and concordance are equal to one; otherwise, it was classified as “unreliable” (STEP 2).

2.6.3 Statistical evaluation of the combined applicability domain and reliability assessment model

To quantify the performance of this model, the exclusion rate and a confusion matrix were analysed. The ideal goal is to maintain a low exclusion rate while excluding all misclassified examples, while at the same time retaining correctly predicted compounds. This corresponds to maximising True Unreliable (TU), minimising False Unreliable (FU), and reducing the number of compounds considered out-of-domain by the applicability domain.

• exclusion rate: retains the majority of compounds to maximise coverage and allows the SVM model to make predictions across a wide range of molecules.

E x c l u s i o n r a t e = n b o f c o m p o u n d O u t ‐ o f ‐ d o m a i n + T U + F U T o t a l c o m p o u n d s × 100

The optimisation of the exclusion rate was performed using an internal validation set, which depends on both the descriptor-domain assessment and the reliability prediction assessment.

• Method of descriptor domain, the number of neighbours used to calculate the mean distance d ¯ k was defined based on the FCFP6 molecular descriptors and was set to k = 3. In addition, K_i≤3 was set as the value for calculating the local threshold.
• For the assessment method of reliability prediction, several molecular fingerprint types were evaluated, among which MACCS fingerprints showed the best performance. In this study, between one and three structural analogues (k = 1–3) with a Tanimoto similarity greater than 0.7 were considered to ensure reliable predictions.

The performance of the combined applicability domain and reliability assessment model was validated on a test set using the confusion matrix metrics (exclusion rate and proportion of mis-predicted compounds correctly flagged) (Table 2).

TABLE 2. Confusion matrix of the combined applicability domain and reliability assessment model.

FR (False Reliable) represents the most critical cases where molecules are mis-predicted by the SVM but not detected so by the workflow, posing a potential risk for model reliability. FU (False Unreliable) are molecules correctly predicted by SVM but rejected by the workflow, which reduces model coverage. TU (True Unreliable) and TRe (True Reliable) correspond to correct decisions, indicating reliable detection.

3 Results

3.1 Curation method

The standardisation of input structures discarded 434 entries (associated with activity measures) due to inability to process a chemical structure into an in silico object. Moreover; 1816 entries were filtered out because remaining multi-constituents after salt removal and/or contained atoms not permitted by the defined list (2.1.3). The dataset obtained after the automatic curation workflow had significantly lower number of entries labelled Agonist and/or Antagonist, compared to the original dataset (Table 3). The projection of the dataset before and after curation clearly demonstrates the distribution of active and cytotoxic labelled samples in both variants of the dataset (Figure 4). This massive shift from the “active” class to the “inactive” class stems from the fact that the original class labels in the Tox21 database are purely based on individual curve-fitting of agonist and antagonist assays, without comparison with corresponding cytotoxicity curves.

TABLE 3. Results of the automated curation method.

After data augmentation with positive entries from BindingDB, the final number of substances in each class was as follows: 291 substances active on thyroid receptors in vitro and 6049 substances inactive on thyroid receptors in vitro.

3.2 Splitting method

The optimisation of the splitting procedure was made in respect to the classic Butina clustering method (detailed comparison not shown). For different clustering thresholds in the range 0.5 ≤ ε ≤ 0.35, the classical method issued fewer singletons than the modified scheme; the classical method formed larger clusters, while the number of clusters issued by the iterative method appeared to converge. When examined by human experts, often the large clusters from the iterative method appeared much more chemically coherent in respect to equivalent clusters (both methods grouped substances around the same seeds) resulting from the classical algorithm, where substances seen as chemically not compatible with the rest of the cluster were present.

The distribution of the reference dataset among the Train and External validation subsets is detailed in Table 4. The counts demonstrated the severe disbalance between active and inactive class, 5.84% vs. 94.16% within the 6,351 entries available.

TABLE 4. Counts by subset during the dataset splitting process.

(*) as fraction of entries in subsets at Level 2, (^) as fraction of entries of the total set at Level 1.

3.3 SVM hyperparameter optimisation for optimal model performance

A set of different molecular fingerprints, including MACCS keys, PubChem, FCFP6, and CFFP6, was employed in the present SVM model to classify substances as active (class 0) or inactive (class 1). The results obtained from the internal validation set are summarised in Supplementary Table 1. In this study, all fingerprints demonstrated good performance. However, the FCFP6 fingerprint, which achieved superior results, particularly in terms of sensitivity, was chosen for the optimal model. FCFP6 fingerprints are also considered advantageous, as they provide additional information about molecular functions that ECFP6 fingerprints may not capture (Supplementary Table 1).

3.4 Dataset representativity and consistency analysis

A t-SNE projection was performed on the training set, the test set, and the combined training + test sets (Figure 5). The visualisation clearly showed that active compounds were well separated from inactive ones in the model’s descriptor space, which explained the model’s strong ability to discriminate between the two classes. Moreover, training set and internal validation set points consistently overlapped for both classes, confirming that the test set was structurally representative of the training data in terms of model descriptors and chemical space coverage.

3.5 Statistics of the training set

A confusion matrix was generated for the training set, comprising 3,344 substances, to evaluate the model’s performance (Supplementary Table 2). The corresponding set of performance metrics is summarised in Table 5. The model correctly predicted all examples in the majority class (3,193 substances) while effectively detecting the minority class, with only 12 false negatives out of 151 instances. This resulted in a sensitivity of 92.05%, a specificity of 100%, a precision (positive predictive value) of 100%, and a negative predictive value (NPV) of 99.63%. These results demonstrate a good balance between sensitivity and specificity with high reliability in predictions for both classes.

TABLE 6. Confusion matrix for the combined applicability domain and reliability assessment model.

In contrast, a naïve baseline model that always predicts the majority (inactive) class would achieve an accuracy of 95.48% (3,193/3,344), which simply reflects the class proportion. The proposed model, with an accuracy of 99.64%, substantially outperforms this baseline and more importantly successfully identifies 92.05% of the rare active compounds while producing zero false positives. This clearly demonstrates the model’s ability to effectively address the severe class imbalance.

The balanced accuracy of 96.03% further highlights the model’s robust performance under class imbalance, as it assigns equal weight to both classes unlike standard accuracy.

The AUC-ROC score (99.91%) (Supplementary Figure 1) highlights the model’s ability to distinguish between classes across all decision thresholds, with a high value indicating excellent discriminative power. Meanwhile, the AUC-PR score (99.21%) (Supplementary Figure 2), which is particularly informative in imbalanced datasets, shows the precision-recall trade-off, confirming that the model excels in identifying the minority class without compromising precision.

3.6 Statistics of the external test set

The model was evaluated on 1,812 previously unseen compounds. It correctly classified 1,729 inactive compounds (2 false positives only) and 71 of 81 active compounds (sensitivity 87.65%), with 10 active compounds classified as false negatives (Supplementary Table 3). The high F1 scores for both classes demonstrate a good precision-recall balance, providing reliable identification of both active and inactive compounds while keeping the false positive rate low (Table 2).

3.7 Comprehensive summary of predictive performance and model generalisation

The statistics from the training folds set, validation folds set, the full training set, and the external validation set are presented in Figure 6. Of the training folds, all metrics (sensitivity, specificity, precision, and negative predictive value) exhibited very high mean values (∼92–100%) with narrow 95% confidence intervals (CI). Mean performance on the validation folds remained high and comparable to that observed on the training folds, indicating no substantial overfitting. The 95% confidence intervals for the metrics of the active class were wider, as expected in a highly imbalanced dataset combined with the limited number of folds, whereas those for the inactive class remained narrow. On the full training set, the same pattern was observed: excellent performance on the inactive class and slightly lower, yet still highly satisfactory, performance on the active class; all metrics fell within the confidence intervals derived from cross-validation. Finally, on the external validation set, performance metrics were very similar to those observed during cross-validation, with results for both the active and inactive classes falling entirely within the corresponding 95% confidence intervals.

3.8 Validation of the combined applicability domain and reliability assessment model on the external validation set

The applicability domain (AD) was computed on the external validation set containing 1,812 substances. When considering the descriptor domain alone, 111 substances were classified as outside of the descriptors domain, whereas all substances were considered in-domain according to the atom-centred fragment approach.

For the remaining substances, the method for assessing the reliability of predictions was applied. Of the 12 substances mispredicted by the SVM model, the prediction of only 7 substances was considered unreliable, while 5 remained reliable. Similarly, of the 1,649 substances correctly predicted by the SVM model, the prediction of 177 substances was classified as unreliable and 1,512 remained reliable (Table 6).

TABLE 5. Summary of model performance statistics for training set, cross-validation set with CIs and external validation before and after the combined applicability domain and reliability assessment model approach.

Results are rounded to two decimal places. *Values of 100.00% indicate performances ≥ 99.95% (or equivalent after rounding).

The study shows that, among these 184 molecules for which predictions were either unreliable, or no confirmation could be obtained, 125 are singletons, meaning that they have no close analogues in the training set with a similarity threshold greater than 70%. They therefore do not share enough characteristics with the other molecules. Next, 59 of the 184 molecules were excluded during the second stage of assessment for the reliability of prediction.

The rate of molecules excluded by the model is 16.2%, while the exclusion rate for misclassified compounds is 58.33%.

The compounds related to a prediction identified as unreliable were removed from the test set prior to recalculating the statistics. The majority class (inactive) maintains an exceptional specificity of 100%, with 1452 true negatives without any false positive (Supplementary Table 4). The minority class (active) also shows a notable improvement with a sensitivity of 92.31% (59 true positives) and a reduction of false negatives to just 5. These results demonstrate that the model improved and is now better able to correctly identify the minority class while maintaining very high precision. Balanced accuracy, reflecting the balance between the two classes, also improves, rising from 92.46% to 96.15%, confirming an overall performance boost (Table 6).

4 Discussion

4.1 Curation and dataset splitting

The ultra-curation of the Tox21 TR dataset started with the standardisation of the description of chemicals. At this point of the workflow filters were applied to exclude metal-organic and multi-constituent entries in the dataset. Although numerous procedures are proposed in the literature, it appears that a “gold-standard” approach is not yet agreed upon, hence we chose to develop our model using the methods and data that are understood by the community. We are willing to address the limitations of these methodologies in the future, especially because many of the compounds discarded from the reference dataset are of interest to industry.

While our expert-driven curation related to cytotoxicity has the advantage of generating a dataset more representative of the “real” activity of chemicals to Thyroid Receptors, it results in difficulties comparing statistics between our model and other models for Thyroid Receptor-mediated activity. Most publicly available models are based on poorly curated datasets, often not fully accounting for the misleading cytotoxicity effect (e.g., see Sapounidou et al., 2023). Thus, these models are trained to detect such cytotoxic substances as “active” and may achieve good statistics on a dataset with a similar level of curation. But such models would perform poorly on our curated dataset where many cytotoxic substances have been assigned as “inactive”, and conversely our model would obtain lower statistics on a less well curated dataset. The importance of creating a truly valid external dataset is thus highlighted.

Similarly, the choice that was taken to merge agonists and antagonists as a single “active” class (because of the low number of positives in each category and data augmentation with BindingDB) makes comparisons with other models more difficult, since other models often predict either agonist or antagonist. A kNN (k-Nearest Neighbours) method is envisaged as a further step to be applied to our model to distinguish between agonists and antagonists in the case where all nearest neighbours of a given test substance predicted as “active” would be identified as either agonists or antagonists.

On the other hand, the integration of Tox21 screening results with BindingDB records inherently introduces experimental heterogeneity, as the active molecules were evaluated across diverse assay types and laboratory conditions. While this variability presents a potential limitation regarding data noise, it also serves as a methodological strength: observing molecular activity across different experimental perspectives prevents the model from overfitting to the specific biological or technical biases of a single assay, thereby enhancing its generalizability. To mitigate the risks associated with this data merging, the dataset was rigorously curated to eliminate any compounds presenting conflicting activity records. Furthermore, it is important to note that the impact of this heterogeneity is confined strictly to the active minority class. The inactive class, which comprises 95.4% of the modelling dataset, was sourced exclusively from standardized Tox21 data and remains unaffected by cross-assay variance. Ultimately, we accepted the residual risk of potential experimental false positives within the active dataset, as this aligns with the conservative nature of our predictive modelling approach. Integrating the BindingDB data also lead to a wider applicability domain. Indeed, the Tox21 dataset is lacking pharmaceutically active substances that are on the other hand well represented in BindingDB.

Our novel iterative Butina clustering scheme proved to be a highly effective approach for partitioning structurally diverse chemical datasets. We did not include comparisons with random or kernel-based splitting strategies, as such methods do not enforce structural coherence and are therefore not appropriate baselines for evaluating clustering performance in highly heterogeneous chemical spaces. Hence, our objective was specifically to focus on improvements over the classical Butina algorithm within a chemically informed clustering framework. Instead, we benchmarked our approach against the classical Butina clustering method as implemented in RDKit.

The iterative scheme resulted in an early depletion of compounds surrounding singletons, highlighting a fundamental methodological difference between the two approaches. This behaviour accounts for the improved chemical coherence of the clusters produced by the iterative method, which, upon expert inspection, were substantially more consistent than those generated by the classical algorithm. In addition, the iterative approach yielded approximately 11% fewer singletons (678 substances) compared with the classical method (761 substances), thereby reducing the number of compounds that had to be assigned to training and test subsets in a quasi-random manner. Because singleton compounds are at high risk of being classified as outside the model’s applicability domain by the kNN-based approach, this reduction represents an additional and practically relevant advantage of the proposed scheme. While this behaviour warrants further investigation through application of the workflow to other datasets, the present results clearly indicate a methodological improvement over the classical Butina clustering strategy.

4.2 New version of Tox21 dataset

During the course of this work, the US EPA released invitroDB v4.1 (Feshuk et al., 2023, Center for Computational Toxicology US EPA ORD, 2023), which is based on the same raw data but with an updated data processing pipeline, compared to invitroDB v3.3 (EPA’s National Center for Computational Toxicology, 2020) that was used in the present work. Only a limited number of substances exhibited changes in activity class labels relative to the previous release (Supplementary Table S6), while the primary updates to the Tox21 TR dataset concerned the inclusion of additional fitted models and modifications to data representation. Despite these enhancements, we retained invitroDB v3.3, as our computational workflows were specifically designed for its data structures.

A key feature of our ultra-curationmethodology is the systematic reassignment of activity labels to Inactive when activity signals observed at concentrations comparable to those associated with significant cytotoxicity. Manual quality assessment of positive substances persisting after the automated curation, further incorporated external criteria including auto-fluorescence and well-documented assay biases. To our knowledge, this integrated and large-scale curation strategy has not been previously applied in the development of predictive models based on the Tox21 TR dataset and thus represents a substantial contribution to QSAR modelling of substances’ activity on the Thyroid Receptors.

4.3 Model and the combined applicability domain and reliability assessment model

Within the regulatory requirements framework, it is essential to have an unambiguous, transparent model supported by a clearly defined scope of applicability. Each step in the modelling process must therefore be explicitly documented, starting with the selection and description of molecular descriptors. Taking these requirements into account, the SVM model for classifying compounds, except organometallic and salt derivatives, as active or inactive not only demonstrated high predictive performance, but also excellent robustness when using FCFP6 fingerprints. In addition, the application of the combined applicability domain and reliability assessment model significantly improved the model’s performance, particularly for the identification of active compounds. Moreover, the overall exclusion rate of 16% is relatively low, demonstrating the model’s relevance and effectiveness, as it retains the majority of compounds while successfully identifying and filtering out unreliable predictions.

However, FCFP6 fingerprints may be more complex for regulatory evaluators to interpret, as they are based on an abstract representation of atomic motifs encoded as bits, with no direct correspondence to explicit chemical descriptors or easily identifiable functional groups. In this context, the use of predefined structural fingerprints (PubChem, MACCS, etc.) is often preferred to facilitate understanding by non-specialists and to enhance transparency during regulatory evaluation. The results presented in Supplementary Table 1 indicate that these fingerprints offer generally satisfactory performance but remain inferior to those obtained with FCFP6. However, they can be optimised, particularly through approaches based on genetic algorithms, to enhance their robustness and predictive efficiency.

In a screening context, where sensitivity and rapid identification of potentially active compounds are priorities, the model based on FCFP6 remains particularly relevant and effective.

The approach to defining the applicability domain (AD) is based on a combination of two complementary methods: one based on the structural domain, the other on the descriptors domain. The descriptor-based approach assesses whether a new compound is sufficiently represented in the local neighbourhood of the training set, thereby increasing confidence in the model’s predictions. Then, the reliability of predictions is assessed using three criteria: (i) the similarity of analogues, based on a predefined threshold and a fixed number of neighbours (k); (ii) accuracy, ensuring thatthe predicted value for each analogue exactly matches its experimental value; and (iii) concordance, by comparing the experimental values of the analogues with the predicted value for the query compound. This reliability assessment method is conservative, considering a prediction reliable only when all identified analogues strictly meet both the accuracy and concordance criteria (Figure 3). It should be noted that this methodology strictly follows the recommendations of the QSAR Applicability Framework (QAF) (OECD, 2024b), a regulatory framework designed to ensure that QSAR models provide reliable and relevant predictions.

Thus, the best performance for evaluating the combined applicability domain and reliability assessment model was obtained using MACCS molecular fingerprints, which are therefore recommended as the default option, particularly for non-expert users. However, for expert users, the methodology remains flexible: other, more specific types of fingerprints can be selected, and the number of neighbours (k) can be adjusted to recompute prediction reliability. Since MACCS key fingerprints provide a relatively general description of molecular similarity, the use of more detailed fingerprints and alternative k values may allow a more tailored and refined assessment of prediction reliability in specific applications. Although this strategy guarantees a high level of confidence, it leads to a very limited scope of applicability, restricting the number of compounds for which a prediction can be considered reliable. It would therefore be particularly relevant to introduce an overall reliability index, incorporating the three key parameters, similarity, accuracy, and concordance in a weighted form. Such an approach would allow the influence of each criterion to be adjusted based on the relevance of the analogues. A similar approach is implemented in the VEGA software to assess the reliability of predictions (Benfenati et al., 2013; Danieli et al., 2023).

5 Conclusion

In conclusion, despite the modelling challenges expected due to severe dataset disbalance, the ultra-curation and advanced dataset splitting provided high quality training opportunities for the SVM model to effectively combines predictive performance, robustness. Moreover, the rigorous definition of the domain of applicability completes the presented model with a conservative assessment of the reliability of predictions. Put together, the iVEMPS-THR complies with newest regulatory requirements for in silico predictions, being characterised with outstanding accuracy, and concordance it ensures a high level of confidence while offering methodological flexibility tailored to the needs of both expert and non-expert users.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ftox.2026.1852502/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Rex FitzGerald, University of Basel, Switzerland

Reviewed by: Fernando Calahorro, University of Southampton, United Kingdom

Sagir Yusuf Ismail, Kano State University of Technology, Nigeria