The reproductive toxicity dataset assembled in this study includes 2,154 small-molecule compounds from the ECHA-C&L Inventory, OECD-eChemPortal and previous literature (Ren et al., 2024; Feng et al., 2021; Zhang et al., 2020, Jiang et cal. 2019). The lists of chemicals are shown in Supplementary Table S1. This dataset includes a varied mix of chemical types—industrial substances, environmental pollutants, and pharmaceuticals. In the European Chemicals Agency(ECHA) database, substances receive Category 1A (known human toxicant), 1B (presumed), or 2 (suspected) based on a weight-of-evidence approach incorporating epidemiological, in vivo, in vitro, and structural-activity data—excluding effects merely secondary to general toxicity (European Chemicals Agency, 2017). Likewise, The Organisation for Economic Co-operation and Development (OECD)- eChemPortal aggregates reproductive hazard classifications (1A/1B/2 and lactation effects) from national and international regulatory sources (OECD, 2008; 2023). Most of the underlying data were derived from rodent studies measuring endpoints like fertility rates, implantation success, offspring development, and congenital malformations, with human evidence primarily influencing Category 1A. Therefore, this dataset provides a good basis for distinguishing reproductively toxic from non-reproductively-toxic compounds to develop computational methods.

This study was conducted in a Python 3.12.2 environment installed with RDKit 2025.3.2, Chemprop 2.1.0, PyTorch 2.7.1, scikit-learn 1.6.1, Pandas 2.2.2, and NumPy 1.26.4. SMILES strings of the 2,154 chemicals were converted with the standard SMILES parser into RDKit molecular objects. Morgan fingerprints (ECFP4) were generated with a radius of 2 and a total length of 2,048 bits in binary (presence/absence) mode, producing ExplicitBitVect-type representations that provided a reliable data foundation for subsequent classification modeling.

This study used eleven classical machine learning methods to predict chemical compounds’ reproductive toxicity. Decision Tree constructs an interpretable tree by recursively splitting on feature thresholds, k-Nearest Neighbors is a “lazy” learner that assigns class based on the majority vote of nearest training samples in Euclidean space, Linear SVM finds a maximum-margin hyperplane for linear separation, Naive Bayes applies Bayes’ theorem under a conditional feature-independence assumption, Logistic Regression fits a sigmoid function to estimate class probabilities, Random Forest ensembles multiple decision trees built on random features and sample subsets to reduce overfitting, Adaptive Boosting (AdaBoost) trains weighted weak learners iteratively and combines them into a strong classifier, Gradient Boosted Decision Tree (GBDT) builds learners sequentially by fitting residual errors, Extra Trees further randomizes split thresholds and features to increase model diversity, Light Gradient Boosting Machine (LightGBM) employs a leaf-wise tree growth strategy for faster training, and XGBoost uses second-order derivative information and regularization to optimize both speed and accuracy. To ensure fair comparison and reproducibility, all models were tuned and evaluated within a unified pipeline. Algorithm implementation and key hyperparameters are shown in Supplementary Table S2.

In this study, we selected the CMPNN framework (Figure 1), given its architecture of enhanced message exchange and boosting, to develop our reproductive toxicity prediction model (ReproTox-CMPNN) and compared it with machine learning methods.

To provide a more formal and precise description, we decomposed CMPNN into five steps, each with its governing equations.

Using Python and RDKit, we generated classical molecular fingerprints and physicochemical descriptors as input features for various machine learning models like Random Forest (RF) and Support Vector Machines (SVM). To overcome the limitations of fingerprints, such as difficulty with activity cliffs and neglecting 3D conformational details, we constructed molecular graphs G and applied ReproTox-CMPNN to learn end-to-end embeddings based on automatic extraction of rich topological and chemical context from the molecular structure (Figure 2).

Model training and evaluation consisted of a rigorous repeated nested cross-validation scheme. In the outer loop, the full dataset was randomly partitioned into five distinct folds, with one fold each time serving as the test set for performance evaluation. In the inner loop, a similar procedure was also repeated five times, with 12.5% each time (10% of the total data) serving for validation/hyperparameter tuning. Nested cross-validation is more robust than 1-layer five-fold cross-validation as it uses inner loop for tuning and a separate set in outer loop for unbiased testing. The 80/20 split in the outer loop is a widely adopted guideline that balances the need for ample training data to learn patterns while retaining a sufficiently large test set to assess generalization; Gholamy et al. (2018) concluded that “p ≈ 80% is empirically the best division into the training and the testing sets.” This system ensured our results’ robustness and reproducibility, laying a solid foundation for future virtual screening and toxicological risk assessment.

Binary cross-entropy loss was minimized using the Adam optimizer, batch size = 50, maximum epochs = 60, with early stopping patience = 10 on validation AUC. Hyperparameters (hidden dimension = 256) were tuned via grid search within each training fold. All experiments were implemented in PyTorch and run on a Linux server equipped with 6 NVIDIA 4090 GPUs.

The comparison of machine learning and ReproTox-CMPNN models was based on accuracy score (ACC), balanced accuracy (BA), Cohen’s Kappa, Matthews correlation coefficient (MCC), F1 score, and the area under the curve (AUC) of the receiver operating characteristic curve (ROC). The designations TP, TN, FP, and FN in calculations refer to the number of true positives, true negatives, false positives, and false negatives, respectively. Compared to accuracy, Cohen’s kappa accounts for the possibility of agreements due to randomness, with p o defined as the observed agreement between two classifiers and p e as the expected agreement by chance. S e n s i t i v i t y = T r u e   p o s i t i v e   r a t e   T P R = T P T P + F N S p e c i f i c i t y = 1 − F a l s e   p o s i t i v e   r a t e   F P R = 1 − F P T N + F P = T N T N + F P A c c u r a c y = T P + T N T P + T N + F P + F N B a l a n c e d   A c c u r a c y = 0.5 ∗   s e n s i t i v i t y + s p e c i f i c i t y F 1   s c o r e = 2 ∗ T P T P + F P ∗ T P T P + F N T P T P + F P + T P T P + F N C o h e n ′ s   K a p p a = p o − p e 1 − p e   w h e r e p o = T P + T N T P + T N + F P + F N   a n d     p e = T P + F P T P + F N + F N + T N F P + T N T P + T N + F P + F N 2 M C C = T P ∗ T N − F P ∗ F N T P + F P T P + F N T N + F P T N + F N

AUC is calculated as the following with FPR as X-axis and TPR as Y-axis. A higher AUC value indicates a greater ability of the model to distinguish between positive and negative cases, i.e., better performance: A U C = ∫ 0 1 T P R t ∗ d F P R t or using trapezoidal rule, A U C ≈ ∑ i = 1 n − 1 F P R i + 1 − F P R i ∗ T P R i + 1 +   T P R i 2

Similarly, for other metrics, a value closer to 1 indicates superior model performance.

The dataset used in this study contains Simplified Molecular Input Line Entry Specifications (SMILES) and a binary classification of reproductive toxicity (1 = Yes, 0 = No). Of the 2,154 chemicals, mainly organic compounds, 1,091 (51%) and 1,063 (49%) are classified as reproductively toxic and non-toxic, respectively. This balanced distribution provides a solid foundation for subsequent QSAR model training and validation.

To better understand the physicochemical features that distinguish reproductive toxicants from non-toxic compounds, we first assessed six common molecular descriptors across our dataset. These six descriptors reflect basic physicochemical characteristics of molecules as well as widely used in QSAR models and drug discovery. As shown in Figure 3, we compared them between toxic and non-toxic compounds. Overall, toxicants exhibit higher median values of and larger dispersion in molecular weight (Weight) and topological polar surface area (TPSA) compared to non-toxic molecules, suggesting that larger and more polar structures are more likely reproductive toxicants. Specifically, the median molecular weight in the non-toxic group is 200 Da compared to 314 Da in the toxic group; similarly, the median non-toxic TPSA is 37 Å² compared to 58 Å² for toxicants with a more right-skewed distribution with more extreme outliers exceeding 200 Å². In terms of lipophilicity, the median of logarithm of the partition coefficient (Log P) is 2.43 in the non-toxic class, compared to 2.76 in the toxic class, although the non-toxic class displays more outliers beyond a very high value (>10), reflecting that enhanced lipophilicity may facilitate membrane permeation and toxic bioaccumulation.

For hydrogen bond donors and acceptors, a slightly more right-skewed distribution in the toxic class indicates potential involvement of hydrogen-bonding interactions in mediating reproductive toxicity. However, comparable distributions between the two classes suggest that these descriptors alone are insufficient for clear-cut classification. Similarly, distributions of rotatable bonds imply limited impact of molecular flexibility on toxicity risk. Taken together, while several molecular descriptors show discernible trends between toxic and non-toxic compounds, these distributions’ overall comparability indicates a need for multivariate modeling. Integrating these descriptors within a comprehensive machine learning framework would be key to robust and generalizable reproductive toxicity predictions.

ReproTox-CMPNN uses a hyperparameter configuration carefully designed to improve predictive performance by balancing convergence and efficiency (Table 1). We used 60 training epochs with a batch size of 50, random enough to avoid local optima while effectively taking advantage of GPU memory to handle molecules of different sizes. We additionally employed an adaptive learning rate through a two-epoch warm-up phase followed by a decay phase.

Learning rate increases linearly from 1e-4 to 1e-3 in the warm-up stage and then decreases exponentially to 1e-4. This method reduces instability during early training, accelerates convergence in the middle stage, and allows for fine-tuning in later stages.

In ReproTox-CMPNN’s architecture, we set the hidden layer dimension to 300 to sufficiently capture complex structural information in molecular graphs while avoiding over-parameterization; in molecular graph neural networks, this hidden layer dimension determines the richness of node (atom) and edge (bond) representations, an important part of the model. ReproTox-CMPNN’s linear layers default to not using bias terms, reducing the number of model parameters and lowering overfitting risk. For molecular graph representations, relative relationships are typically more important than absolute offsets, and this bias-free design enables the model to focus more on learning the relative feature importance in molecular structures. The message passing depth is set to 3, allowing each atom to perceive neighboring atoms up to three hops away. For most drug molecules, this depth adequately captures key structural information and chemical environments while avoiding over-smoothing and overfitting problems that deeper networks might introduce.

We obtained results presented in Figure 4A which show the ROC curves of the ReproTox-CMPNN model under five-fold cross-validation. The AUC values for folds 0–4 are 0.929, 0.942, 0.962, 0.956, and 0.939, respectively, resulting in a mean AUC of 0.946 with a standard deviation of 0.013. These results demonstrate excellent discriminative power and low variability across different folds, indicating robust stability and generalization performance. Figure 4B, showing the evolution of AUC during training epochs, illustrates the model’s rapid improvement within the first 10 epochs, with AUC rising from approximately 0.85 to over 0.92, after which the curve gradually flattens and stabilizes around 0.93. This behavior suggests fast convergence and no significant overfitting during later stages of training coupled with a consistent predictive accuracy throughout.

Overall, the CMPNN-based model achieves high and consistent AUC scores in cross-validation and demonstrates quick convergence and resilience during the training process, thus providing a solid foundation for reliable reproductive toxicity prediction of new compounds.

To provide a comprehensive benchmark of predictive methods on our reproductive toxicity dataset, we first evaluated eleven classical machine learning algorithms alongside the ReproTox-CMPNN model. As summarized in Table 2, ReproTox-CMPNN achieved an AUC of 0.946, substantially outperforming the next-best methods Random Forest (0.825) and Linear SVM (0.823). This clear margin underscores ReproTox-CMPNN’s superior ability to capture intricate molecular topology and chemical context, leading to markedly improved discrimination between toxic and non-toxic compounds.

In terms of overall classification metrics, ReproTox-CMPNN attained a mean accuracy of 0.857 with a standard deviation of 0.019 and a balanced accuracy (BA) of 0.856 with a standard deviation of 0.018, significantly higher than those of other models (e.g., Random Forest 0.825).

These results indicate that ReproTox-CMPNN delivers not only excellent overall predictive power but also maintains high balance between positive (toxic) and negative (non-toxic) classes, effectively mitigating the biases arising from class imbalance. Furthermore, its sensitivity (mean 0.823 and standard deviation 0.076) and specificity (mean 0.890 and standard deviation 0.085) both exceed 0.80, demonstrating reliable recall of toxicants and exclusion of non-toxicants.

Examining more stringent agreement measures, ReproTox-CMPNN’s Cohen’s Kappa (mean 0.713 and standard deviation 0.037) and Matthews Correlation Coefficient (mean 0.721 and standard deviation 0.033) are well above those of traditional models (most of which fall between 0.50 and 0.65), highlighting strong consistency and correlation with true labels. The F1 score of mean of 0.846 and standard deviation of 0.018, further reflects a balanced trade-off between precision and recall, particularly excelling in identifying toxic compounds. In contrast, simpler approaches such as K-Nearest Neighbors (AUC = 0.717) or Naive Bayes (AUC = 0.783) exhibit inferior discrimination and stability.

In summary, the ReproTox-CMPNN model, through its advanced representation of molecular structures, significantly surpasses various conventional machine learning algorithms, offering a powerful and robust framework for reproductive toxicity prediction with promising applicability in risk assessment pipelines.

As shown in Table 3, in recent work regarding reproductive toxicity prediction, Jiang et al. (2019) used a much more extensive dataset (1,823 chemicals) compared to previous publications, multiple endpoints such as sperm reduction and infertility, and six machine learning methods to develop more reliable models. Their study recommended the SVM model using Molecular Access System Keys Fingerprints (MACCSFP), which generated an AUC of 0.900 and an accuracy score of 0.836. As previous work mainly focused on frequentist methods, Zhang et al. (2020) investigated Naive Bayes (NB) model together with six molecular descriptors and ten types of fingerprints. Their best model resulted in an AUC of 0.888 and an accuracy score of 0.830. Feng et al. (2021) used three machine learning methods with nine molecular fingerprints to build more ensemble models. The model they recommended had an AUC of 0.920 and an accuracy score of 0.844. Ren et al. (2024) developed a deep learning fragment-based graph transformer network (FGTN) model to predict reproductive toxicity, taking pre-generated fragments as nodes and bonds between fragments as edges with a super-molecule-level node to connect all fragment nodes. The FGTN model showed an AUC of 0.914 and an accuracy score of 0.861.

Due to differences in datasets and data splitting, a direct comparison between the aforementioned models and the ReproTox-CMPNN model is not possible. However, a simple comparison can still provide some valuable information. ReproTox-CMPNN remarkably outperformed other models in AUC and was very close to the FGTN model in ACC and MCC, exceeding the remaining models. These results suggest that the ReproTox-CMPNN model is both reliable and robust (Table 3).

CMPNN’s strong performance across multiple molecular property benchmarks is attributed to its enhanced messaging framework: unlike traditional MPNNs that focus solely on node-to-node communication, CMPNN utilizes a communicative kernel to reinforce interactions between both node and edge features, complemented by a message booster, generating richer molecular representations. As it is inherently based on 2D molecular graphs, CMPNN may miss stereochemistry and spatial relationships crucial for predicting properties like enantiomer-specific activity. Furthermore, because our training set mainly comprises organic environmental pollutants, our model’s generalizability to other chemical domains such as metal complexes or inorganic salts remains to be demonstrated. Finally, since our current evaluation focuses on specific endpoints, validating CMPNN across a broader spectrum of toxicity readouts such as DNA damage, endocrine disruption will be needed to further illuminate its applicability.

In this study, we built prediction models on a comprehensive reproductive toxicity dataset with 2,154 chemicals and compared the performance of the deep learning ReproTox-CMPNN model with other methods. Our results show that ReproTox-CMPNN learning results exceed the current best ML models in both embedding quality and predictive accuracy, making it a new state-of-the-art method in this field. ReproTox-CMPNN’s deep capture of multi-level molecular relationships offers an efficient and reliable computational tool for rapid chemical safety screening and risk assessment.

Future work includes training ReproTox-CMPNN to predict multiple toxicity endpoints or reproductive toxicity of chemical mixtures as well as combining ReproTox-CMPNN with mask language model (MLM) for protein embedding to further improve predictions.