Prediction of adverse biological effects of chemicals using knowledge graph embeddings

Abstract

We have created a knowledge graph based on major data sources used in ecotoxicological risk assessment. We have applied this knowledge graph to an important task in risk assessment, namely chemical effect prediction. We have evaluated nine knowledge graph embedding models from a selection of geometric, decomposition, and convolutional models on this prediction task. We show that using knowledge graph embeddings can increase the accuracy of effect prediction with neural networks. Furthermore, we have implemented a fine-tuning architecture which adapts the knowledge graph embeddings to the effect prediction task and leads to a better performance. Finally, we evaluate certain characteristics of the knowledge graph embedding models to shed light on the individual model performance.

1. Introduction

Ecotoxicology is a multidisciplinary field that studies the potentially adverse toxicological effects of chemicals on organisms, starting at molecular level to individuals, sub-populations, communities and ecosystems. One major societal contribution of ecotoxicology is ecological risk assessments, which compare environmental concentrations of chemicals with existing laboratory effect data to evaluate the ecosystem health status. While laboratory experiments are thus crucial, they are both labour intensive and result in a high number of animal testing. Therefore, the development of modelling techniques for extrapolating from existing laboratory effect data is a major effort in the field of ecotoxicology.

A very important challenge in ecotoxicology risk assessment is the interoperability of the disparate data sources, formats and vocabularies. The use of Semantic Web technologies and (RDF-based) knowledge graphs [6] can address this challenge and facilitate the orchestration of these datasets. Hence, extrapolation or prediction models can benefit from an integrated view of the data and the background knowledge provided by a knowledge graph. The use of knowledge graphs also enables the use of the available infrastructure to perform automated reasoning, explore the data via semantic queries, and compute semantic embeddings for machine learning prediction.

In this work we have created the Toxicological Effect and Risk Assessment Knowledge Graph (TERA) and implemented a prediction model over this knowledge graph to extrapolate adverse biological effects of chemicals on organisms. Here, we limit ourselves to binary effect prediction of mortality (shortened to effect prediction), i.e., where there is a chance that a chemical can affect a species in a lethal way. The work and evaluation conducted in this paper is driven by the following research question: does the use of contextual information in the form of knowledge graph embeddings brings added value in the prediction of adverse biological effects?

Our contributions can be summarized as follows:

TERA aims at consolidating the relevant information to the ecological risk assessment domain. TERA integrates several disparate datasets and enables a unified (semantic) access. The formats of these data sources vary from tabular, to RDF files and SPARQL endpoints over public linked data. We have exploited external resources (e.g., Wikidata [76]) and ontology alignment methods (e.g., LogMap [33]) to discover equivalences between the data sources.

We have designed and implemented a model tailored to binary lethal chemical effect prediction. This model relies on TERA and builds upon existing knowledge graph embedding models. Moreover, it supplies the knowledge graph embedding models with additional information. This is used to tailor the embeddings to this specific task.

We have evaluated nine knowledge graph embedding (KGE) models, together with a naive baseline on the binary chemical effect prediction task. This evaluation includes four data sampling strategies which highlight the different settings of chemical effect prediction (i.e., the test data contains unseen chemical-organism pairs where: (a) the chemical and the organism may be known (but not in previously seen pairs), (b) the chemical is unknown, (c) the organism is unknown, and (d) both the chemical and the organism are unknown).

These contributions are openly shared. A snapshot of the TERA knowledge graph is available on Zenodo [53] (https://doi.org/10.5281/zenodo.3559865) and the source scripts for creating TERA are available on GitHub (https://github.com/NIVA-Knowledge-Graph/TERA). Finally, the scripts to reproduce the conducted evaluation in this paper are also available on GitHub (https://github.com/NIVA-Knowledge-Graph/KGs_and_Effect_Prediction_2020).

This paper extends our preliminary work presented in the In-Use Track of the 18th International Semantic Web Conference [51]. We have (i) extended TERA with new sources (Encyclopedia of Life (EOL), MeSH, and a larger part of ChEMBL) and provided detailed steps about its creation; (ii) created a more robust prediction model with nine (up from three) embedding algorithms supported and a task-specific embedding fine-tuning strategy; and (iii) conducted a more comprehensive evaluation with all combinations of KGE models and sampling strategies totalling 648 data points (324 for each prediction model).

The rest of the paper is organized as follows. Section 2 introduces essential concepts to the subsequent sections. Section 3 introduces the use case where the knowledge graph and prediction models are applied. Section 4 introduces related work. The creation of the knowledge graph is described in Section 5. Section 6 introduces the prediction models, while Section 7 presents the evaluation of these models. Section 8 elaborates on the contributions and discusses future directions of research. Finally, the Appendix gives an overview of the knowledge graph embedding models used in this work.

2. Preliminaries

In this section we introduce important background concepts that will be used throughout the paper. Table 1 contain the most important symbols.

Table 1

Key symbols and acronyms used throughout the paper

Symbol	Definition
RDF	Resource Description Framework
OWL	Web Ontology Language
SPARQL	SPARQL Protocol and RDF Query Language
KG	Knowledge graph
KGE	Knowledge graph embedding
t	A triple
$\mathit{sb}$	The subject of a triple
$\mathit{ob}$	The object of a triple
p, r	The predicate/relation of a triple
e	A KG entity
$\mathcal{T}$	The set of KG triples
$\mathcal{E}$	The set of KG entities
$\mathcal{R}$	The set of KG relations
$\mathcal{L}$	The set of literal values
e	The vector representation of an entity or relation
k	The dimension of a vector
$\mathit{SF}$	The scoring function of a KGE model
$\mathit{PT}$	Pre-trained KGE-based model
$\mathit{FT}$	Fine-tuning KGE-based model
s	A species
c	A chemical
S	Refers to species
C	Refers to chemicals
κ	Chemical concentration

2.1. Ecotoxicological terminology

Taxonomy in this work refers to a species classification hierarchy. Any node in a taxonomy is called a taxon. Species is a taxon which is also a leaf node in the taxonomy. An Organism denotes an individual living organism which is an instance of a species. Chemicals or compounds are unique isotopes of substances consisting of two or more atoms. Effect, used in this work as short form for chemical effect, refers to the response of an organism (or population) to a chemical at a specific concentration. Endpoint1

¹

Not to be confused with SPARQL endpoint.

denotes a measured effect on the test population at a certain time; e.g., lethal concentration to 50% of test population (LC50) measured at 48 hours. Note that, an experiment can have several endpoints, e.g., LC50 at 48 hours and LC100 at 96 hours (lethal concentration for all test organisms). See Table 2 for the most common endpoints.

2.2. Ontology-enhanced knowledge graphs

In this work we consider the most broadly accepted notion of knowledge graph within the Semantic Web: an ontology enhanced RDF-based knowledge graph (KG) [32]. This kind of knowledge graph enables the use of the available Semantic Web infrastructure, including SPARQL engines and OWL reasoners.2

²

RDF, RDFS, OWL and SPARQL are standards defined by the W3C: https://www.w3.org/standards/semanticweb/.

Thus, in our setting, KGs are composed by RDF triples in the form of $⟨\mathit{sb}\mathrm{,}\mathit{p}\mathrm{,}\mathit{ob}⟩\in \mathcal{E}\times \mathcal{R}\times \mathcal{E}\cup \mathcal{L}$ ,3

³

$\mathcal{E}$ is the set of all classes and instances, $\mathcal{R}$ is the set of all properties, while $\mathcal{L}$ represents the set of all literal values.

where $\mathit{sb}$ represents a subject (a class or an instance), p represents a predicate (a property) and $\mathit{ob}$ represents an object (a class, an instance or a literal). KG entities (i.e., $\mathcal{E}\cup \mathcal{R}$ : classes, properties and instances) are represented by an URI (Uniform Resource Identifier).

An (ontology-enhanced) KG can be split into a TBox (terminology) and an ABox (assertions). The TBox is composed by triples using RDF Schema (RDFS) constructors like class subsumptions and property domain and range; and OWL constructors like disjointness, equivalence and property inverses.4

⁴

Note that the Web Ontology Language (OWL) [27] also enables the creation of complex axioms that are translated/serialized into more than one triple: https://www.w3.org/TR/owl2-mapping-to-rdf/.

The ABox contains assertions among instances, including OWL equality and inequality, and semantic type definitions. Table 5 shows several examples of TBox and ABox triples.

2.3. Ontology alignment

Ontology alignment is the process of finding mappings or correspondences between a source and a target ontology or knowledge graph [23,66]. These mappings typically represent equivalences or broader/narrower relationships among the entities of the input ontologies. In the ontology matching community [61], mappings are exchanged using the RDF Alignment format [18]; but they can also be interpreted as standard OWL axioms (e.g., [24,35]). In this work we treat ontology alignments as OWL axioms (e.g., triple ${\mathit{t}}_{13}$ in Table 5). An ontology matching system (e.g., LogMap [34]) is a program that, given as input two ontologies or knowledge graphs, generates as output a set of mappings (i.e., an alignment) M.

2.4. Embedding models

Knowledge graph embedding (KGE) [63,78] plays a key role in link prediction problems where it is applied to knowledge graphs to resolve missing facts in largely connected knowledge graphs, such as DBpedia [44]. Biomedical link prediction is another area where embedding models have been applied successfully (e.g., [1,5]).

The embeddings of the entities in a KG are commonly learned by (i) defining a scoring function over a triple, which is typically proportional to the probability of the existence of that triple in the KG,5

⁵

For the embedding process, we focus on triples where $\mathit{o}\in \mathcal{E}$ is a class or an instance.

i.e., $\mathit{SF}:\mathcal{E}\times \mathcal{R}\times \mathcal{E}\to \mathbb{R}$ , $\mathit{SF}\propto \mathit{P}\mathrm{(}⟨\mathit{sb}\mathrm{,}\mathit{p}\mathrm{,}\mathit{ob}⟩\in \mathrm{KG}\mathrm{)}$ ; and (ii) minimizing a loss function (i.e., deviation of the prediction of the scoring function with respect to the truth available in the KG). More specifically, KGE models (i) initialize the entities in a triple $⟨\mathit{sb}\mathrm{,}\mathit{p}\mathrm{,}\mathit{ob}⟩$ into a vector representation ${\mathbf{e}}_{\mathit{sb}}\mathrm{,}{\mathbf{e}}_{\mathit{p}}\mathrm{,}{\mathbf{e}}_{\mathit{ob}}\in {\mathbb{R}}^{\mathit{k}}\text{ or }{\mathbb{C}}^{\mathit{k}}$ , where k is the dimension of the vector; (ii) apply a scoring function to $\mathrm{(}{\mathbf{e}}_{\mathit{sb}}\mathrm{,}{\mathbf{e}}_{\mathit{p}}\mathrm{,}{\mathbf{e}}_{\mathit{ob}}\mathrm{)}$ ; and (iii) adapt the vector representations to improve the scoring and minimize the loss.

Several knowledge graph embedding models have been proposed. In this work, we used models of three major categories: decomposition models, geometric models, and convolutional models.6

⁶

The interested reader please refer to [63] for a comprehensive survey.

The decomposition models represent the triples of the KG into a one-hot 3-order tensor and apply matrix decomposition to learn entity vectors. Geometric models, also known as translational, try to learn embeddings by defining a scoring function where the predicate in the triple act as a geometric translation (e.g., rotation) from subject to object. Convolutional models, unlike previous models, learn entity embedding with non-linear scoring functions via convolutional layers.

3. Ecotoxicological risk assessment and adverse biological effect prediction

The task of ecotoxicological risk assessment is to study the potential hazardous effects of chemicals on organisms from individuals to ecosystems. In this context, risk is the result of the intrinsic hazards of a substance on species, populations or ecosystems, combined with an estimate of the environmental exposure, i.e., the product of exposure and effect (hazard).

Fig. 1.

Simplified ecological risk assessment pipeline.

Figure 1 shows a simplified risk assessment pipeline. Exposure data is gathered from analysis of environmental concentrations of one or more chemicals, while effects (hazards) are characterized for a number of species in the laboratory as a proxy for more ecologically relevant organisms. These two data sources are used to calculate the so-called risk quotient (RQ; ratio between exposure and effects). The RQ for one chemical or the mixture of many chemicals is used to identify chemicals with the highest RQs (risk drivers), identify relevant modes of action7

⁷

The mode of action describes the molecular pathway by which a chemical causes physiological change in an organism.

(MoA) and characterize detailed toxicity mechanisms for one or more species (or taxa). Results from these predictions can generate a number of new hypotheses that can be investigated in the laboratory or studied in the environment. Note that, this risk assessment pipeline is a simplified version of the one in use at the Norwegian Institute for Water Research,8

⁸

NIVA: https://www.niva.no/en.

however, similar methodologies are used across regulatory risk assessment pipelines.

Table 2

The most frequent endpoints in ECOTOX [74] chemical effect data

Endpoint	Frequency	Description
NR	0.21	Not reported
NOEL	0.17	No-observable-effect-level
LC50	0.16	Lethal concentration for $50\mathrm{\%}$ of test population
LOEL	0.14	Lowest-observable-effect-level
NOEC	0.05	No-observable-effect-concentration
EC50	0.05	Effective concentration for $50\mathrm{\%}$ of test population
LOEC	0.04	Lowest observable effect concentration
BCF	0.03	Bioconcentration factor
NR-LETH	0.02	Lethal to $100\mathrm{\%}$ of test population
LD50	0.02	Lethal dose for $50\mathrm{\%}$ of test population
Other	0.11

The chemical effect data is gathered during laboratory experiments, where a sub-population of a single species is exposed to an increasing concentration of a toxic chemical. The endpoints of the experiments are recorded at chemical concentrations and time after exposure. These endpoints are categorized into several categories, e.g., lethality rate of test population (see Table 2).

Ecological risk assessment methods require a large amount of these experimental data to give an accurate depiction of the long term risk to an ecosystem. The data must cover the relevant chemicals and species present in the ecosystem, e.g., an ecological risk assessment of agricultural runoff in Norway will mostly concern pesticides and waterflees, copepods, and frogs, among other species [42]. Just with a few relevant chemicals and species the search space becomes immense and performing laboratory experiments becomes unfeasible. Thus, it is essential to develop in silico methods to extrapolate new chemical-species effects from known combinations. We differentiate among two types complementary strategies: (i) highly specialized (restricted in chemical and species domains) models to predict chemical concentrations that will have an effect on a test species, and (ii) models that produce rankings of highly representative chemical-species pair hypothesis which can be used by a laboratory to perform targeted experiments. In this paper we focus on the latter strategy, using a method based on knowledge graph embeddings. Methods that fall into the first strategy are introduced in Section 4.1.

4. Related work

This section will cover related work from ecotoxicology and knowledge graph based prediction.

4.1. Toxicity extrapolation

There are two main research areas in toxicology to extrapolate chemical effects, i.e., Quantitative Structure-Activity Relationship (QSAR) and read-across. QSAR modelling try to find a relationship between the structure of a chemical and the chemical’s biological activity (cf. reviews [22,26]). This relationship is described using derived chemical features. Some features are simple, e.g., octanol-water partition coefficient or logP, others concern the entire chemical, e.g., chemical fingerprints. The basis of the QSAR relationship is usually modeled as polynomial equations. Parthasarathi and Dhawan [59] take this further by using the logarithm of chemical concentration to achieve a polynomial relationship: $log\mathrm{(}1\mathrm{/}\mathit{\kappa }\mathrm{)}=\mathit{f}\mathrm{(}\mathit{\pi }\mathrm{)}+\mathit{g}\mathrm{(}\mathit{\sigma }\mathrm{)}$ , $\mathit{f}\in {\mathit{P}}_2$ and $\mathit{g}\in {\mathit{P}}_1$ ( ${\mathit{P}}_{\mathit{n}}$ is a polynomial of nth degree), where κ is the chemical concentration while π and σ denote the derived chemical features hydrophobicity9

⁹

Measure of the absence of attraction to water.

and electronic effects in the molecule, respectively. The drawback of these models is the applicability domains. Usually, a QSAR model considers a small set of chemicals (10ths to 100ths) and one single species. This means that new features and relationships need to be developed for each species and each chemical group.

The read-across methods try to mitigate these drawbacks, mainly by considering extrapolation of the effect at the chemical and species levels. Similar to QSAR models, read-across of chemicals use the chemical features to create similarity measures between chemicals to justify the read-across of chemical effects. The read-across in the species domain is harder. Species do not tend to have easily derived features. Therefore, genetic similarity has emerged as a viable option. Sequence Alignment to Predict Across Species Susceptibility (SeqAPASS), developed by the United States Environmental Protection Agency (U.S. EPA), is an example of such an approach [20,41]. SeqAPASS uses a large amount of data available for humans, mice, rats, and zebrafish to extrapolate to areas with lower coverage.

4.2. Embedding models

In this work, we use nine KGE models across three categories of models. Here, we will give a brief introduction to the models, while a more extended explanation of the models is found in the Appendix. The interested reader please refer to [63] for a comprehensive survey.

The three categories of models are decomposition, geometric, and convolutional [63]. The decomposition models are DistMult, ComplEx, and HolE. DistMult models the score of a triple as the vector multiplication of the representation of each subject, predicate and object [83]. ComplEx uses the same scoring function as DistMult, however, in a complex vector space, such that it can handle inverse relations [73]. HolE is based on holographic embeddings [56], however, it has been shown that HolE is equivalent to ComplEx [30].

The geometric models are TransE, RotatE, pRotatE, and HAKE. TransE is the base of a whole family of models and scores triples based on the translation from subject to object using the representation of the predicate [10]. RotatE is similar to TransE, however, the translation using the predicate is done by rotating it (via Euler’s identity) [70]. Furthermore, pRotatE is a baseline for RotatE where the modulus in Euler’s identity is ignored [70]. Finally, the hierarchical-aware model, HAKE, where entities at each level in the hierarchy is at equal distance from the origin and relations at a level is modeled as rotation [86].

The convolutional models take a deep learning approach to the task of KGE. We use ConvKB [55] and ConvE [19], which are similar with slightly different architectures. They have shown good performance given the relative small number of parameters.

Although quite a few KGE models have been proposed, the adopted ones are either classic models or can achieve state-of-the-art performance in some benchmarks. They are representative of mainstream techniques, and have been widely adopted in KGE research and applications [63]. Thus, the benefits and shortcomings of the KGE models analysed in this study provide good evidence of the general performance of this type of models in a complex prediction task, i.e., adverse biological effect of chemicals on organisms.

4.3. Using KGE for prediction

Our focus to use KGE models is to predict if a chemical has a lethal effect on an organism. KGE models have been explored in the biomedical domain to solve similar predictions tasks (e.g., finding relationships between diseases, drugs, genes, and treatments). Several works have shown improvements in results by using KGE models for prediction, e.g., [1,5,46]. Chen et al. [15] used random walks over networks to perform drug-target predictions. The ChEMBL and DrugBank KGs have also been used to predict chemical mode of action (MoA) of anticancer drugs with high performance on benchmark datasets [82].

Opa2vec [68] and Blagec et al. [8] have developed embedding models to improve similarity-based prediction in the biomedical domain, while OpenBioLink [12] has created a framework for evaluating models in the biomedical domain.

EL Embeddings [40] and opa2vec [68] present new semantic embedding methods for KGs with expressive logic expressions (i.e., OWL ontologies) to predict protein interaction. The former utilizes complex geometric structures to model the logic relationships between entities, while the later learns a language model from a corpus extracted from the ontology. OWL2Vec* [13] also learns a language model from an ontology and applies the computed embeddings into two prediction tasks: class subsumption and class membership. OWL2Vec* has also been used to predict the plausibility of ontology alignments [14].

To the best of our knowledge there is no work using link prediction or KGE models to support ecotoxicological effect prediction. This study will give novel insights and empirical results of KGE models in this new domain.

5. TERA knowledge graph

One major challenge in ecological risk assessment processes is the interoperability of data. In this section, we introduce the Toxicological Effect and Risk Assessment (TERA), an ontology-enhanced RDF-based knowledge graph that aims at providing an integrated view of the relevant data sources for risk assessment.10

¹⁰

Resources to create and access TERA: https://github.com/NIVA-Knowledge-Graph/TERA.

The initial inspiration for TERA was the aid of ecotoxicological effect prediction where access to disparate resources was required (see Section 5.3). However, by integrating these sources into a KG, we were also able to directly apply TERA into the prediction process by leveraging knowledge graph embedding models (see Section 5.4).

The data sources integrated into TERA vary from tabular and RDF files to SPARQL endpoints over public linked data. The sources currently integrated into TERA are: (i) biological: NCBI Taxonomy, Encyclopedia of Life, and Wikidata mappings (∼500k species); (ii) chemical: PubChem, ChEMBL, MeSH, and Wikidata mappings (∼110M compounds); and (iii) biological effects: ECOTOXicology Knowledgebase (∼1M results, ∼12k compounds, ∼13k species), and system-generated mappings. These three distinct parts make up the sub-KGs of TERA, i.e., (i) the Taxonomy sub-KG ( ${\mathrm{KG}}_{\mathit{S}}$ ), (ii) the Chemical sub-KG ( ${\mathrm{KG}}_{\mathit{C}}$ ), and (iii) the Effects sub-KG ( ${\mathrm{KG}}_{\mathit{E}}$ ). The different processes to transform and integrate these sources into TERA are shown in Fig. 2.

Fig. 2.

Data sources and processes to create the TERA knowledge graph.

A snapshot of TERA is available on Zenodo [53], where licenses permit.11

¹¹

EOL: Various Creative commons (CC), NCBI: Creative Commons CC0 1.0 Universal (CC0 1.0), ECOTOX: No restrictions, PubChem: Open Data Commons Open Database License, ChEMBL: CC Attribution, MeSH: Open, Courtesy of the U.S. National Library of Medicine, Wikidata: CC0 1.0.

PubChem and ChEMBL are not included in the snapshot due to size constraints; these can be downloaded from the National Institutes of Health12

¹²

ftp://ftp.ncbi.nlm.nih.gov/pubchem/RDF/

and European Bioinformatics Institute,13

¹³

ftp://ftp.ebi.ac.uk/pub/databases/chembl/ChEMBL-RDF/

respectively. The subgraph of TERA used for prediction is available alongside the chemical effect prediction models in our GitHub repository.14

¹⁴

https://github.com/NIVA-Knowledge-Graph/KGs_and_Effect_Prediction_2020

Table 5 shows several examples of RDF triples from TERA.15

¹⁵

Prefixes associated to the URI namespaces of entities in TERA: et: (ECOTOXicology knowledgebase), ncbi: (NCBI taxonomy), eol: (Encyclopedia of Life), mesh: (Medical Subject Heading), compound: (PubChem compound), descr: (PubChem descriptors), vocab: (PubChem vocabulary), inchikey: (InChIKey identifiers), envo: (Environment Ontology) cheminf: (Chemical information ontology), chembl: (ChEMBL), chembl_m: (ChEMBL molecule subset), chembl_t: (ChEMBL target subset), wd: (WikiData entities), wdt: (Wikidata properties), qudt: (Quantities, Units, Dimensions and Types Catalog), snomedct: (SNOMED CT ontology), and bp: (Biological PAthway eXchange ontology). owl:, rdfs:, rdf: and xsd: are prefixes referring to W3C standard vocabularies.

5.1. Dataset overview

TERA, as mentioned above, is constructed by gathering a number of sources about chemicals, species and chemical toxicity, with a diverse set of formats including tabular data, RDF dumps and SPARQL endpoints.

Biological effect data of chemicals. The largest publicly available repository of effect data is the ECOTOXicology knowledgebase (ECOTOX) developed by the US Environmental Protection Agency [74]. This data is gathered from published toxicological studies and limited internal experiments. The dataset consists of $1\mathit{M}$ experiments covering $12\mathit{k}$ chemicals and $13\mathit{k}$ species,16

¹⁶

Version dated Sep. 15, 2020.

implying a chemical–species pair converge of maximum $\sim 0.6\mathrm{\%}$ . The resulting endpoint from an experiment is categorised in one of a plethora of predefined endpoints (see Table 2 above).

Table 3

ECOTOX database tests example

test_id	reference_number	test_cas	species_number	organism_habitat
1147366	12448	134623 (diethyltoluamide)	1 (Pimephales promelas)	Water

Table 4

ECOTOX database results example

result_id	test_id	endpoint	effect	conc1_mean	conc1_unit
102570	1147366	LC50	MOR	110000	μg/L

Tables 3 and 4 contain an excerpt of the ECOTOX database. ECOTOX includes information about the chemicals and species used in the tests. This information, however, is limited and additional (external) resources are required to complement ECOTOX.

Chemicals. The ECOTOX database uses an identifier called CAS Registry Number assigned by the Chemical Abstracts Service to identify chemicals. The CAS numbers are proprietary, however, Wikidata [76] (indirectly) encodes mappings between CAS numbers and open identifiers like InChIKey, a 27-character hash of the International Chemical Identifier (InChI) which encodes chemical information uniquely [31].17

¹⁷

While InChI is unique, InChiKey is not, and collisions have greater than zero probability [79].

Wikidata also provides mappings to well known databases like PubChem, ChEMBL and MeSH, which include relevant chemical information such as chemical structure, structural classification and functional classification.

Taxonomy. ECOTOX contains a taxonomy18

¹⁸

In the context of the paper “taxonomy” typically refers to a classification of organisms.

(of species), however, this only considers the species represented in the ECOTOX effect data. Hence, to enable extrapolation of effects across a larger taxonomic domain, we include the NCBI Taxonomy [64]. This taxonomy data source consists of a number of database dump files, which contains a hierarchy for all sequenced species, which equates to around $10\mathrm{\%}$ of the currently known life on Earth and is one of the most comprehensive taxonomic resources. For each of the taxa (species and classes), the taxonomy defines a handful of labels, the most commonly used of which are the scientific and common names. However, labels such as authority can be used to see the citation where the species was first mentioned, while synonym is a alternate scientific name, that may be used in the literature.

Species traits. As an analog to chemical features, we use species traits to expand the coverage of the knowledge graph. Apart from taxonomic classifications, traits are the most important information to identify species and will be of great importance when predicting the effect on the species.

The traits we have included in the knowledge graph are the habitat, endemic regions, and presence (and classifications of these). This data is gathered from the Encyclopedia of Life (EOL) [57], which is available as a property graph. Moreover, EOL uses external definitions of certain concepts, and mappings to these sources are available as glossary files. In addition to traits, researchers may be interested in species that have different conservation statuses, e.g., if the population is stable or declining, etc. This data can also be extracted from EOL.

5.2. Dataset preprocessing

In this section we present the different steps to extract, transform and integrate the source datasets into the main TERA components and sub-KGs. All data is transformed using custom mappings (scripts) from the sources to RDF triples. Table 5 shows an excerpt of the triples in TERA.

Table 5

Example triples from the TERA knowledge graph. For space reasons, we have added the full id or label for some of the entities using footnote marks where ¹inchikey:MMOXZBCLCQITDF-UHFFFAOYSA-N, ²Pimephales, ³Cyprinidae, ⁴Headwater, ⁵Benzamides, ⁶Insect Repellents, ⁷CHRNA3, ⁸CHRNB4, ⁹DETA-20, ¹⁰DETA Epichlorohydrin, ¹¹Has component, ¹²Triclocarban, ¹³Trichlorocarbanilide-containing product, ¹⁴Similar to, ¹⁵3-Chloromethyl-N,N-diethylbenzamide

#	subject	predicate	object
Effects sub-KG
${\mathit{t}}_1$	et:test/1147366	et:compound	et:chemical/134623
${\mathit{t}}_2$	et:test/1147366	et:species	et:taxon/1
${\mathit{t}}_3$	et:test/1147366	et:hasResult	et:result/102570
${\mathit{t}}_4$	et:result/102570	et:endpoint	et:endpoint/LC50
${\mathit{t}}_5$	et:result/102570	et:effect	et:effect/Mortality
${\mathit{t}}_6$	et:taxon/1	rdf:type	et:taxon/Pimephales
${\mathit{t}}_7$	et:taxon/Pimephales	rdfs:subClassOf	et:taxon/Cyprinidae
${\mathit{t}}_8$	et:taxon/1	et:latinName	“Pimephales promelas”
${\mathit{t}}_9$	et:taxon/1	et:commonName	“Fathead Minnow”
${\mathit{t}}_{10}$	et:taxon/1	et:speciesGroup	et:group/Fish
${\mathit{t}}_{11}$	et:taxon/1	et:rank	et:rank/species
${\mathit{t}}_{12}$	et:chemical/134623	rdfs:label	“diethyltoluamide”
Entity Mappings
${\mathit{t}}_{13}$	et:taxon/1	owl:sameAs	ncbi:taxon/90988
${\mathit{t}}_{14}$	ncbi:taxon/90988	owl:sameAs	wd:Q2700010
${\mathit{t}}_{15}$	wd:Q2700010	owl:sameAs	eol:211492
${\mathit{t}}_{16}$	et:chemical/134623	owl:sameAs	wd:Q408389
${\mathit{t}}_{17}$	wd:Q408389	owl:sameAs	chembl_m:CHEMBL1453317
${\mathit{t}}_{18}$	wd:Q408389	owl:sameAs	compound:CID4284
${\mathit{t}}_{19}$	wd:Q408389	owl:sameAs	mesh:D003671
${\mathit{t}}_{20}$	wd:Q408389	owl:sameAs	inchikey:MMOXZBCLC… 1
Taxonomy sub-KG
${\mathit{t}}_{21}$	ncbi:taxon/90988	rdf:type	ncbi:taxon/51137 2
${\mathit{t}}_{22}$	ncbi:taxon/90988	rdf:type	ncbi:division/10
${\mathit{t}}_{23}$	ncbi:taxon/90988	ncbi:scientific_name	“Pimephales promelas”
${\mathit{t}}_{24}$	ncbi:taxon/90988	ncbi:rank	ncbi:species
${\mathit{t}}_{25}$	ncbi:taxon/51137	rdfs:subClassOf	ncbi:taxon/7953 3
${\mathit{t}}_{26}$	ncbi:division/10	rdfs:label	“Vertebrates”
${\mathit{t}}_{27}$	ncbi:division/10	owl:disjointWith	ncbi:division/1
${\mathit{t}}_{28}$	ncbi:division/1	rdfs:label	“Invertebrates”
${\mathit{t}}_{29}$	eol:211492	eol:habitat	envo:00000153 4
Chemical sub-KG
${\mathit{t}}_{30}$	mesh:D003671	mesh:broaderDescriptor	mesh:D001549 5
${\mathit{t}}_{31}$	mesh:D003671	mesh:pharmacologicalAction	mesh:D007302 6
${\mathit{t}}_{32}$	chembl_m:CHEMBL1453317	chembl:hasTarget	chembl_t:CHEMBL1907594 7
${\mathit{t}}_{33}$	chembl_t:CHEMBL1907594	chembl:relSubsetOf	chembl_t:CHEMBL3137273 8
${\mathit{t}}_{34}$	compound:CID89845769 9	vocab:hasParentCompound	compound:CID4284
${\mathit{t}}_{35}$	compound:CID131721069 10	cheminf:CHEMINF_000478 11	compound:CID4284
${\mathit{t}}_{36}$	compound:CID131721069	rdf:type	bp:SmallMolecule
${\mathit{t}}_{37}$	compound:CID7547 12	vocab:is_active_ingredient_of	snomedct:411346009 13
${\mathit{t}}_{38}$	compound:CID131721069	cheminf:CHEMINF_000480 14	compound:CID10751691 15

Fig. 3.

Example of an ECOTOX test and related triples.

5.2.1. Effects sub-KG construction

The effect data in ECOTOX consist of two parts, i.e., test definitions and results associated with the test definitions (see Tables 3 and 4, respectively). The important columns of a test are the chemical and the species used. Other columns include metadata, but these are optional and often empty. Each result is composed by an endpoint, an effect, and a concentration (with a unit) at which the endpoint and effect are recorded.

This tabular data in ECOTOX is transformed into triples that form the effects sub-KG in TERA ( ${\mathrm{KG}}_{\mathit{E}}$ ). Note that a test can have multiple results. A subset of the effect triples are listed in Table 5 (see triples ${\mathit{t}}_1$ – ${\mathit{t}}_{12}$ ). A graphical representation for an effect test and its result is also shown in Fig. 3.

ECOTOX contains metadata about the species and chemicals used in the experiments. This metadata is also included in TERA to facilitate the alignment with other resources (see Section 5.2.2).

• The ECOTOX metadata file species.txt includes common and Latin names, along with a (species) ECOTOX group (see triples ${\mathit{t}}_8$ – ${\mathit{t}}_{10}$ in Table 5). This group is a categorization of the species based on ECOTOX use cases. Prefixes and abbreviations like sp., var. are removed from the label names.

• The full hierarchical lineage19

  ¹⁹

  As defined by U.S. EPA. Note that species hierarchies are contested among researchers.

  is also available in the metadata file species.txt. Each column represents a taxonomic level, e.g., genus or family. If a column is empty, we construct an intermediate classification; for example, Daphnia magna has no genus classification in the data, then its classification is set to Daphniidae genus (family name + genus, actually called Daphnia). We construct these classifications to ensure the number of levels in the taxonomy is consistent (see triples ${\mathit{t}}_6$ and ${\mathit{t}}_7$ in Table 5). Note that when adding triples such as ${\mathit{t}}_{11}$ in Table 5, we also add a taxonomic rank to facilitate the querying for a specific taxonomic level.

• The ECOTOX source file chemicals.txt includes chemical metadata and it is handled similarly to species.txt. The file includes chemical name (see  ${\mathit{t}}_{12}$ in Table 5) and a (chemical) ECOTOX group.

For the units in the effect data, e.g., chemical concentrations (mg/L, mol/L, mg/kg, etc.), we reuse the QUDT 1.120

²⁰

QUDT 1.1: http://linkedmodel.org/catalog/qudt/1.1/

ontologies. When an unit such as mg/L is not defined, we define it according to Listing 1.

Listing 1.

Unit definition of mg/L using QUDT

5.2.2. Alignment with state-of-the-art tools

ECOTOX database provides proprietary chemical identifiers (i.e., CAS numbers) and internal ECOTOX ids for species. In order to extrapolate effects across a larger set of chemicals and species than those available in ECOTOX, TERA integrates taxonomy and trait data from NCBI and EOL, and chemical data from PubChem, ChEMBL and MeSH.

Alignment between ECOTOX and the NCBI Taxonomy. There does not exist a complete and public alignment between the 23,439 ECOTOX species and the 1,830,312 the NCBI Taxonomy species.21

²¹

There are a total of 27,133 and 2,246,074 taxa in ECOTOX and NCBI, respectively. However, we focus on species, i.e., instances.

We have used three methods, two state-of-art ontology alignments systems and a baseline, to align ECOTOX and the NCBI Taxonomy: (i) LogMap [33,34], (ii) AgreementMakerLight (AML) [25], and (iii) a string matching algorithm based on Levenshtein distance [45]. LogMap and AML were chosen since they have performed well across many datasets in the Ontology Alignment Evaluation Initiative (e.g., [2,3,61]). Most mappings in our setting are expected to be lexical, therefore, we also selected a purely lexical matcher to evaluate if more sophisticated systems like LogMap and AML bring an additional value.

Due to the large size of the NCBI Taxonomy, we needed to split NCBI into manageable chunks to enable the use of ontology alignment systems. Fortunately, this can be easily done by considering the species division, e.g., mammal or invertebrate. This divides the NCBI Taxonomy into 11 distinct parts, which can be aligned to the taxonomy in ECOTOX.

Table 6

Alignment results for ECOTOX-NCBI. #M: number of mappings (at instance level), R: Recall, ${\mathrm{P}}^{\approx }$ : estimated precision

Method	1-to-1 mappings
	#M	R	${\mathrm{P}}^{\approx }$
LogMap	20,585	0.81	0.87
AML	14,148	0.77	0.94
String similarity ( $\mathrm{>}0.8$ )	20,423	0.76	0.87
Consensus ( $\mathrm{LogMap}\cap \mathrm{AML}$ )	12,740	0.76	0.98
$\mathbf{LogMap}\mathbf{\cup }\mathbf{AML}$	21,145	0.83	0.86

Note that it is expected an entity from ECOTOX to match to a single entity in the NCBI Taxonomy, and vice-versa. Hence, 1-to-N and N-to-1 alignments were filtered according to the system computed confidence. A partial mapping curated by experts can be obtained through the ECOTOX Web.22

²²

ECOTOX interface: https://cfpub.epa.gov/ecotox/search.cfm.

We have gathered a total of 2,321 mappings for validation purposes. Table 6 shows the alignment results over the ground truth samples for the 1-to-1 (filtered) system mappings. We report number of mappings (#M), Recall (R) and estimated precision ( ${\mathrm{P}}^{\approx }$ ) with respect to the known entities in the incomplete ground truth, assuming only 1-to-1 mappings are valid. ${\mathrm{P}}^{\approx }$ is calculated as $$\begin{array}{ll}\text{(1)}& \hfill {\mathit{P}}^{\approx }=& |{\mathit{M}}^{\approx }\cap {\mathit{M}}_{\mathrm{ref}}|\mathrm{/}\left|{\mathit{M}}^{\approx }\right|\mathrm{,}\hfill {\mathit{M}}^{\approx }=& \{⟨{\mathit{e}}_{\mathit{e}}\mathrm{,}\mathtt{owl}:\mathtt{sameAs}\mathrm{,}{\mathit{e}}_{\mathit{n}}⟩\in \mathit{M}\mid \\ \text{(2)}& \hfill & {\mathit{e}}_{\mathit{e}}\in {\mathcal{E}}_{\mathit{e}}^{\mathrm{ref}}\vee {\mathit{e}}_{\mathit{n}}\in {\mathcal{E}}_{\mathit{n}}^{\mathrm{ref}}\}\mathrm{,}\end{array}$$ where ${\mathit{M}}_{\mathrm{ref}}$ is the (incomplete) reference mapping set and M is the set of generated mappings between entities ${\mathit{e}}_{\mathit{e}}\in {\mathcal{E}}_{\mathit{e}}$ from ECOTOX and entities ${\mathit{e}}_{\mathit{n}}\in {\mathcal{E}}_{\mathit{n}}$ from the NCBI Taxonomy, ${\mathcal{E}}_{\mathit{e}}^{\mathrm{ref}}\subseteq {\mathcal{E}}_{\mathit{e}}$ and ${\mathcal{E}}_{\mathit{n}}^{\mathrm{ref}}\subseteq {\mathcal{E}}_{\mathit{n}}$ are the sets of entities that appear in the reference mappings. Thus, ${\mathit{M}}^{\approx }$ is defined as a subset of mappings from M involving entities in the reference mapping set ${\mathit{M}}_{\mathrm{ref}}$ . Recall is defined in the standard way as $$\begin{array}{}\text{(3)}& \mathit{R}=|\mathit{M}\cap {\mathit{M}}_{\mathrm{ref}}|\mathrm{/}|{\mathit{M}}_{\mathrm{ref}}|.\end{array}$$ Note that, the recall will be the same for M and ${\mathit{M}}^{\approx }$ .

We have selected the union of the 1-to-1 equivalence23

²³

There is no need for more complex mappings in this use case.

mappings computed by AML and LogMap to be integrated within TERA, as they represent the mapping set with the best recall with a reasonable estimated precision. This choice was made by considering the large uncertainty of downstream applications (effect prediction and risk assessment), where we prefer a larger coverage of the domain. See triple ${\mathit{t}}_{13}$ in Table 5 for an example of a system computed mapping between ECOTOX and the NCBI Taxonomy.

Listing 2.

Construct taxon mapping between Wikidata and, NCBI and EOL. wd:Q16521 is the class of all taxa, while wdt:P31, wdt:P685 and wdt:P830 are the relations instance of, NCBI Taxonomy ID and Encyclopedia of Life ID, respectively

We use Wikidata as source of alignments between the NCBI Taxonomy and EOL, and among the used chemical datasets. Alignments are extracted via Wikidata’s query interface (i.e., SPARQL endpoint).24

²⁴

Wikidata endpoint: https://query.wikidata.org/sparql.

The data in Wikidata concerning species and chemicals are in large parts manually curated [77] and will have a low error rate, comparatively to using the automated ontology alignment systems.

Alignment between the NCBI Taxonomy and EOL. In order to include in TERA trait data from EOL, we need to establish an alignment between EOL and the NCBI Taxonomy. We have constructed equivalence triples between the NCBI Taxonomy and EOL identifiers using Wikidata. The species identifiers are available as literals in Wikidata. Therefore, we concatenate them with the appropriate namespace. Listing 2 represents the SPARQL CONSTRUCT query used against the Wikidata endpoint. Here, we query Wikidata for instances of taxa, thereafter adding optional triple patterns for NCBI Taxonomy and EOL identifiers which are added as owl:sameAs triples to TERA.

Examples of resulting mapping triples are shown in ${\mathit{t}}_{14}$ – ${\mathit{t}}_{15}$ in Table 5. The proportion of species in Wikidata where this mapping exists is $49\mathrm{\%}$ .

Alignment between chemical entities. The mapping between ECOTOX chemical identifiers (CAS Registry Numbers) to Wikidata entities enables the alignment to a vast set of chemical datasets, e.g., PubChem, ChEBI, KEGG, ChemSpider, MeSH, UMLS, to name a few. The construction of equivalence triples between CAS, ChEMBL, MeSH, PubChem and Wikidata identifiers is shown in Listing 3. As for the case of species identifiers, the literal representing a chemical identifier is concatenated with the corresponding namespace. For the CAS Registry Numbers we also remove the hyphens to match ECOTOX notation. Examples of resulting mapping triples are shown in ${\mathit{t}}_{16}$ – ${\mathit{t}}_{20}$ in Table 5.

Listing 3.

Construct chemical mapping between Wikidata and ECOTOX, ChEMBL, MeSH and PubChem. wdt:P31 is the predicate for instance of and wd:Q11173 is the class of all chemical compounds. wdt:P231, wdt:P592, wdt:P486, wdt:P662 and wdt:P235 are the relations for CAS Registry Number, ChEMBL ID, MeSH ID, PubChem CID and InChIKey, respectively

These mappings are not complete, but for some the coverage is large. Out of the chemicals used in ECOTOX, $73\mathrm{\%}$ have an equivalence in Wikidata (through the CAS registry numbers). Moreover, Wikidata chemicals has $4\mathrm{\%}$ ChEMBL identifiers, $0.5\mathrm{\%}$ MeSH identifiers, $55\mathrm{\%}$ PubChem identifiers, and $95\mathrm{\%}$ InChiKey identifiers.

5.2.3. Taxonomy sub-KG construction

The Taxonomy sub-KG ( ${\mathrm{KG}}_{\mathit{S}}$ ) integrates data from the NCBI Taxonomy and the EOL trait data. The integration of the NCBI Taxonomy into the TERA knowledge graph is split into several sub-tasks.

• We load the hierarchical structure included in the NCBI Taxonomy file nodes.dmp. The columns of interest are the taxon identifiers of the child and parent taxon, along with the rank of the child taxon and the division where the taxon belongs. We use this to create triples like ${\mathit{t}}_{21}$ – ${\mathit{t}}_{22}$ and ${\mathit{t}}_{24}$ – ${\mathit{t}}_{25}$ in Table 5.

• To aid alignment between the NCBI Taxonomy and the ECOTOX identifiers, we add the synonyms found in names.dmp. Here, the taxon identifier, its name and name type are used to create triples like ${\mathit{t}}_{23}$ in Table 5. Note that a taxon in the NCBI Taxonomy can have several synonyms while a taxon in ECOTOX usually has two, i.e., common name and scientific name.

• Finally, we add the labels of the divisions found in divisions.dmp (see triples ${\mathit{t}}_{26}$ and ${\mathit{t}}_{28}$ ). We also add disjointness axioms among unrelated divisions, e.g., triple ${\mathit{t}}_{27}$ in Table 5.

We use the TraitBank from EOL [58] to add species traits to TERA. The TraitBank is modeled as a property graph and can be accessed as a neo4j database or via a set of tabular files. To integrate the TraitBank into TERA we validate the identifiers used in EOL and convert to URIs. If an identifier is not a valid URI, we replace invalid symbols. A trait example is shown as triple ${\mathit{t}}_{29}$ in Table 5. The EOL TraitBank also includes subsumption definitions (i.e., via rdfs:subClassOf) for a large portion of traits. These subsumptions can be downloaded separately and are added to TERA in a similar way as mentioned above.

5.2.4. Chemical sub-KG construction

The Chemical sub-KG ( ${\mathrm{KG}}_{\mathit{C}}$ ) is created from PubChem [38], ChEMBL [29], and MeSH [47]. These datasets are available for download as RDF triples. In addition, ChEMBL and MeSH can be accessed through the EBI and MeSH SPARQL endpoints, respectively.

The chemical subset of PubChem is used since information about chemicals is standardized in PubChem, while information about substances is not. In this subset we use: (i) component information, i.e., what are the building blocks of the chemical or parts of a mixture; (ii) type assertions, which either link to ChEBI or describe the type of molecule, e.g., small or large; (iii) role assertions, which describe additional attributes or relationships of the chemical, e.g., FDAApprovedDrug; and (iv) drug products, which link to the clinical data in SNOMED CT [7]. Examples of these can be seen in triples ${\mathit{t}}_{35}$ , ${\mathit{t}}_{36}$ and ${\mathit{t}}_{37}$ in Table 5.

Parent chemical data in PubChem is limited to permutations e.g., bonds, polarity, and part of mixtures axioms (triple ${\mathit{t}}_{34}$ in Table 5). Therefore, we use the hierarchical data about chemicals from MeSH. In addition to this data, we create similarity triples between chemicals. This is impractical to download, but can be calculated on demand. We add similarity triples to TERA where the Tanimoto (Jaccard) distance between the chemical fingerprints (gathered using PubChemPy [71]) is $⩾0.9$ ,25

²⁵

Default value used in PubChem [37].

see triple  ${\mathit{t}}_{38}$ in Table 5.

ChEMBL contains facts about bioactivity of chemicals. This contributes in assessing the danger of a chemical. In TERA, we use the mode of action (MoA) and target (receptor targeted by MoA; triple ${\mathit{t}}_{32}$ in Table 5). These targets are organized in a hierarchy using chembl:relSubsetOf relations (see triple ${\mathit{t}}_{33}$ ). The receptors will link to which organism it belongs to, however, we leave the inclusion of this information for future work.

We use the entire MeSH dataset in TERA. MeSH is organised as several hierarchies. The most prominent classifications are based on chemical groups and the intended use of the chemicals. Triples ${\mathit{t}}_{30}$ and ${\mathit{t}}_{31}$ in Table 5 show examples of chemical group and functional classifications.

Listing 4.

Query to select all species, chemicals, concentrations and units, where the species is endemic to the Oslofjord

5.3. TERA for data access

TERA covers knowledge and data relevant to the ecotoxicological domain and enables an integrated semantic access across data sets. In addition, the adoption of an RDF-based knowledge graph enables the use of an extensive range of Semantic Web infrastructure (e.g., reasoning engines, ontology alignment systems, SPARQL query engines).

The data integration efforts and the construction of TERA go in line with the vision in the computational risk assessment communities (e.g., Norwegian Institute for Water Research’s Computational Toxicology Program (NCTP)), where increasing the availability and accessibility of knowledge enables optimal decision making.

The knowledge in TERA can be accessed via predefined queries26

²⁶

Predefined queries are typically abstractions of SPARQL queries.

(e.g., classification, sibling, and name queries, and fuzzy queries over the species names) and arbitrary SPARQL queries. The (final) output is flexible to the task, and can be given either as a graph or in tabular format. Listing 4 shows an example query to extract the chemicals and concentrations, at which, the species in the Oslofjord experience lethal effects.

5.4. TERA for effect prediction

TERA is used as background knowledge in combination with machine learning models for chemical effect prediction. TERA’s sub-KGs play different roles in effect prediction. The rich semantics of the species and chemical entities in the Taxonomy sub-KG ( ${\mathrm{KG}}_{\mathit{S}}$ ) and the Chemical sub-KG ( ${\mathrm{KG}}_{\mathit{C}}$ ), respectively, are embedded into low-dimensional vectors; while the Effects sub-KG ( ${\mathrm{KG}}_{\mathit{E}}$ ) provides the training samples for the prediction model. Each sample is composed of a chemical, a species, a chemical concentration, and the outcome or endpoint of the experiment. More details are given in Section 6, where the effect prediction model is built upon state-of-the-art knowledge graph embedding models.

Table 7

Densities and entropies of benchmark datasets. TERA ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ are the chemical and species parts of TERA, while ${\mathrm{KG}}_{\mathit{C}}^{\prime }$ and ${\mathrm{KG}}_{\mathit{S}}^{\prime }$ denote the parts of TERA used in prediction in Section 7

Dataset	RD	ED	RE	EE	AD
TERA ${\mathrm{KG}}_{\mathit{C}}$	$2.3\times {10}^5$	5.5	3.0	24	$4.6\times {10}^{-7}$
TERA ${\mathrm{KG}}_{\mathit{S}}$	$6.6\times {10}^4$	5.1	2.7	23	$3.7\times {10}^{-7}$
TERA ${\mathrm{KG}}_{\mathit{C}}^{\prime }$	$6.9\times {10}^3$	8.6	2.3	17	$7.7\times {10}^{-5}$
TERA ${\mathrm{KG}}_{\mathit{S}}^{\prime }$	$3.8\times {10}^2$	15	2.3	14	$8.9\times {10}^{-4}$
YAGO3-10	$2.9\times {10}^4$	18	2.0	20	$7.1\times {10}^{-5}$
FB15k-237	$1.3\times {10}^3$	43	4.5	16	$1.3\times {10}^{-3}$
WN18	$8.4\times {10}^3$	7.4	2.1	16	$9.0\times {10}^{-5}$
WN18RR	$8.5\times {10}^3$	4.5	1.5	19	$5.5\times {10}^{-5}$

Table 7 shows the sparsity-related measures of common benchmark datasets27

²⁷

YAGO3-10 [69], FB15k-237 [9], WN18 [48] and WN18RR [19].

and TERA’s ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ (triples involving literals are removed). We follow Pujara et al. [62] and calculate the relational density, $\mathit{RD}=|\mathcal{T}|\mathrm{/}|\mathcal{R}|$ , and entity density, $\mathit{ED}=2|\mathcal{T}|\mathrm{/}|\mathcal{E}|$ , where $\mathcal{T}$ , $\mathcal{R}$ , and $\mathcal{E}$ are the sets of triples, relations, and entities in the knowledge graph, respectively. The entity entropy (EE) and the relation entropy (RE) indicate whether there are biases (the lower EE or RE, the larger bias) in the triples in the KG [62], and are calculated as $$\begin{array}{}\text{(4)}& \mathit{P}\mathrm{(}\mathit{r}\mathrm{)}={\displaystyle \frac{|\mathit{t}.\mathit{p}=\mathit{r}|}{|\mathit{T}|}}\mathrm{,}\text{(5)}& \mathit{P}\mathrm{(}\mathit{e}\mathrm{)}={\displaystyle \frac{|\mathit{t}.\mathit{s}\mathit{b}=\mathit{e}|+|\mathit{t}.\mathit{o}\mathit{b}=\mathit{e}|}{|\mathit{T}|}}\mathrm{,}\text{(6)}& \mathrm{RE}={\displaystyle \sum }_{\mathit{r}\in \mathcal{R}}-\mathit{P}\mathrm{(}\mathit{r}\mathrm{)}log\mathrm{(}\mathit{P}\mathrm{(}\mathit{r}\mathrm{)}\mathrm{)}\mathrm{,}\text{(7)}& \mathrm{EE}={\displaystyle \sum }_{\mathit{e}\in \mathcal{E}}-\mathit{P}\mathrm{(}\mathit{e}\mathrm{)}log\mathrm{(}\mathit{P}\mathrm{(}\mathit{e}\mathrm{)}\mathrm{)}\mathrm{,}\end{array}$$ where $|\mathit{t}.\mathit{p}=\mathit{r}|$ is the number of triples with r as predicate, and $|\mathit{t}.\mathit{s}\mathit{b}=\mathit{e}|+|\mathit{t}.\mathit{o}\mathit{b}=\mathit{e}|$ is the number triples with e as subject or object.

In addition, we calculate the absolute density of the graph, which is  $\mathit{AD}=|\mathcal{T}|\mathrm{/}\mathrm{(}|\mathcal{E}|\mathrm{(}|\mathcal{E}|-1\mathrm{)}\mathrm{)}$ . This is the ratio of edges to the maximum number of edges possible in a simple directed graph [17].

High RD and low RE typically lead to a worse performance, while high ED and low EE often lead to better link prediction performance (e.g., [19]). In Table 7 we can see that the density and entropy values are in between those for YAGO3-10 and FB15k-237, which typically lead to worse and better predictive performance, respectively [19]. This shows that TERA is a suitable background knowledge to extrapolate effect data and, at the same time, an interesting dataset to benchmark state-of-the-art knowledge graph embedding models. Note that using the full TERA (i.e., ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ ), according to RD, will be more challenging than using the reduced TERA fragments (i.e., ${\mathrm{KG}}_{\mathit{C}}^{\prime }$ and ${\mathrm{KG}}_{\mathit{S}}^{\prime }$ ) for prediction. Full details of the construction of ${\mathrm{KG}}_{\mathit{C}}^{\prime }$ and ${\mathrm{KG}}_{\mathit{S}}^{\prime }$ are given in Section 7.1.1.

6. Adverse biological effect prediction

The aim of chemical effect prediction is to extrapolate exiting data to new combinations of (possibly unknown) chemicals and species. In this section we present three classification models used to predict the adverse biological effect of chemicals on species: (i) a multilayer perceptron (MLP) model (our baseline), (ii) the baseline model fed with pre-trained KG embeddings, (iii) a model that simultaneously trains the baseline model and the KGE models (i.e., it fine-tunes the KG embeddings). A MLP was chosen as baseline as it is a basic model where additional components and penalties can be easily added and assessed as we do in our third model (see Section 6.3).

The models have three inputs, namely a chemical c, a species s, and a chemical concentration κ (denoted ${\mathit{x}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}$ ). The output is a binary value that represents whether the chemical at the given concentration has a lethal effect on the species: $$\begin{array}{}\text{(8)}& \hfill {\mathit{y}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}=\left\{\begin{array}{cc}1\hfill & \mathit{c}\text{ is lethal to }\mathit{s}\text{ at }\mathit{\kappa }\mathrm{,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}\right.& \end{array}$$ Note that the effect can have a more fine-grained categorization (endpoints $\mathrm{LC}\mathit{x}$ , $\mathrm{LD}\mathit{x}$ , $\mathrm{EC}\mathit{x}$ ,28

²⁸

If effect is mortality (e.g., see Table 4).

and NR-LETH in Table 2). Without losing the generality in introducing and evaluating our effect prediction methods, we simplify the effect into two cases: “lethal” and “non-lethal”.

Notation. Throughout this section we use bold lower case letters to denote vectors while matrices are denoted as bold upper case letters. The vector representation of an entity and a relation are noted as ${\mathbf{e}}_{\mathit{e}}$ and ${\mathbf{e}}_{\mathit{p}}$ , respectively. These vectors are either in ${\mathbb{R}}^{\mathit{k}}$ or ${\mathbb{C}}^{\mathit{k}}$ , where k is the embedding dimension.

Fig. 4.

Baseline model. Inputs: c, s, κ as in Equation (9); Outputs: $\hat{\mathit{y}}$ as in Equation (15).

6.1. Baseline model

Our baseline prediction model is a multilayer perceptron (MLP) with multiple hidden layers. ${\mathit{n}}_{\mathit{c}}$ hidden layers are appended to the embedding ${\mathbf{e}}_{\mathit{c}}$ of the chemical c, ${\mathit{n}}_{\mathit{s}}$ hidden layers are appended to the embedding ${\mathbf{e}}_{\mathit{s}}$ of species s, and ${\mathit{n}}_{\mathit{\kappa }}$ hidden layers appended to the real valued chemical concentration κ. Thereafter, n hidden layers are further appended to the output of the previous hidden layers concatenated. Specifically, the model can be expressed by the following equations (with ${\mathit{x}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}$ as input): $$\begin{array}{}\text{(9)}& {\mathbf{y}}_{\mathit{c}}^0={\mathbf{e}}_{\mathit{c}}\mathrm{,}{\mathbf{y}}_{\mathit{s}}^0={\mathbf{e}}_{\mathit{s}}\mathrm{,}{\mathit{y}}_{\mathit{\kappa }}^0=\mathit{\kappa }\text{(10)}& {\mathbf{y}}_{\mathit{c}}^{\mathit{h}}=\mathit{ReLu}\mathrm{(}{\mathbf{y}}_{\mathit{c}}^{\mathit{h}-1}{\mathbf{W}}_{\mathit{c}}^{\mathit{h}}+{\mathbf{b}}_{\mathit{c}}^{\mathit{h}}\mathrm{)}\mathrm{,}\mathit{h}\in \{0\mathrm{,}\dots \mathrm{,}{\mathit{n}}_{\mathit{c}}\}\text{(11)}& {\mathbf{y}}_{\mathit{s}}^{\mathit{h}}=\mathit{ReLu}\mathrm{(}{\mathbf{y}}_{\mathit{s}}^{\mathit{h}-1}{\mathbf{W}}_{\mathit{s}}^{\mathit{h}}+{\mathbf{b}}_{\mathit{s}}^{\mathit{h}}\mathrm{)}\mathrm{,}\mathit{h}\in \{0\mathrm{,}\dots \mathrm{,}{\mathit{n}}_{\mathit{s}}\}\text{(12)}& {\mathbf{y}}_{\mathit{\kappa }}^{\mathit{h}}=\mathit{ReLu}\mathrm{(}{\mathbf{y}}_{\mathit{\kappa }}^{\mathit{h}-1}{\mathbf{W}}_{\mathit{\kappa }}^{\mathit{h}}+{\mathbf{b}}_{\mathit{\kappa }}^{\mathit{h}}\mathrm{)}\mathrm{,}\mathit{h}\in \{0\mathrm{,}\dots \mathrm{,}{\mathit{n}}_{\mathit{\kappa }}\}\text{(13)}& {\mathbf{y}}^0=[{\mathbf{y}}_{\mathit{c}}^{{\mathit{n}}_{\mathit{c}}}\mathrm{,}{\mathbf{y}}_{\mathit{s}}^{{\mathit{n}}_{\mathit{s}}}\mathrm{,}{\mathbf{y}}_{\mathit{\kappa }}^{{\mathit{n}}_{\mathit{\kappa }}}]\text{(14)}& {\mathbf{y}}^{\mathit{h}}=\mathit{ReLu}\mathrm{(}{\mathbf{y}}^{\mathit{h}-1}{\mathbf{W}}^{\mathit{h}}+{\mathbf{b}}^{\mathit{h}}\mathrm{)}\mathrm{,}\mathit{h}\in \{1\mathrm{,}\dots \mathrm{,}\mathit{n}\}\text{(15)}& \hat{\mathit{y}}=\mathit{\sigma }\mathrm{(}{\mathbf{y}}^{\mathit{n}}{\mathbf{W}}^{\mathit{n}}+{\mathbf{b}}^{\mathit{n}}\mathrm{)}\end{array}$$ ${\mathbf{e}}_{\mathit{c}}\mathrm{,}{\mathbf{e}}_{\mathit{s}}\in {\mathbb{R}}^{\mathit{k}}$ in (9) denote the embeddings of c and s respectively, and are calculated as $$\begin{array}{}\text{(16)}& {\mathbf{e}}_{\mathit{c}}={\mathit{\delta }}_{\mathit{c}}{\mathbf{W}}_{\mathit{c}}\mathrm{,}{\mathbf{e}}_{\mathit{s}}={\mathit{\delta }}_{\mathit{s}}{\mathbf{W}}_{\mathit{s}}\end{array}$$ where ${\mathit{\delta }}_{\mathit{c}}$ and ${\mathit{\delta }}_{\mathit{s}}$ denote the one-hot encoding vectors of the chemical entity c (w.r.t. all the entities in ${\mathcal{E}}_{\mathit{C}}$ from ${\mathrm{KG}}_{\mathit{C}}$ ) and the species entity s (w.r.t. all the entities in ${\mathcal{E}}_{\mathit{S}}$ from ${\mathrm{KG}}_{\mathit{S}}$ ), respectively;29

²⁹

${\mathit{\delta }}_{\mathit{c}}\in {\mathbb{R}}^{|{\mathcal{E}}_{\mathit{C}}|}$ , where ${\mathit{\delta }}_{\mathit{c}}^{\mathit{i}}=1$ if c is the ith chemical in ${\mathcal{E}}_{\mathit{C}}$ , else 0. ${\mathit{\delta }}_{\mathit{s}}$ is defined similarly.

${\mathbf{W}}_{\mathit{c}}\in {\mathbb{R}}^{|{\mathcal{E}}_{\mathit{C}}|\times \mathit{k}}$ and ${\mathbf{W}}_{\mathit{s}}\in {\mathbb{R}}^{|{\mathcal{E}}_{\mathit{S}}|\times \mathit{k}}$ are embedding transformation matrices to learn. (10), (11) and (14) represent the hidden layers, where $\mathit{ReLu}$ denotes the rectifier function (i.e., $\mathit{ReLu}\mathrm{(}\mathit{x}\mathrm{)}=max\mathrm{(}0\mathrm{,}\mathit{x}\mathrm{)}$ ), ${\mathbf{W}}_{\mathit{c}}^{\mathit{t}}$ , ${\mathbf{W}}_{\mathit{s}}^{\mathit{t}}$ and ${\mathbf{W}}^{\mathit{t}}$ denote the weights, ${\mathbf{b}}_{\mathit{c}}^{\mathit{t}}$ , ${\mathbf{b}}_{\mathit{s}}^{\mathit{t}}$ and ${\mathbf{b}}^{\mathit{t}}$ denote the biases. $[·\mathrm{,}·]$ in (13) denotes vector concatenation. σ in (15) denotes the sigmoid function (i.e., $\mathit{\sigma }\mathrm{(}\mathit{x}\mathrm{)}=1\mathrm{/}\mathrm{(}1+exp\mathrm{(}-\mathit{x}\mathrm{)}\mathrm{)}$ ). Note that a dropout and a normalization layer is stacked after each hidden layer for regularization.

We differentiate between two settings of the baseline model (see Fig. 4):

• Simple setting. Figure 4a shows the model without embedding transformation layers, i.e., ${\mathit{n}}_{\mathit{s}}={\mathit{n}}_{\mathit{c}}={\mathit{n}}_{\mathit{\kappa }}=0$ , and $\mathit{n}=1$ .

• Complex setting. The complex model shown in Fig. 4b introduces transformation layers on the embeddings and chemical concentration input. These transformations aim at extracting the important information in the inputs and disregard the redundant information based on the output.

In the experiments we refer to the baseline models as Simple one-hot and Complex one-hot, depending on the selected MLP setting.

6.2. Baseline model with pre-trained KG embeddings

This models relies on pre-trained embeddings of chemicals and species computed using state-of-the-art KGE models (see Section 4.2 and the Appendix for an overview). A (different) KGE model is applied to the chemicals ${\mathrm{KG}}_{\mathit{C}}$ and the species ${\mathrm{KG}}_{\mathit{S}}$ .

These pre-trained KG embeddings are then given as input instead of the one-hot encoding vectors in the baseline model. We replace the trainable matrices ${\mathbf{W}}_{\mathit{c}}$ and ${\mathbf{W}}_{\mathit{s}}$ in Equation (16) by the matrices composed of embeddings by the respective KGE models. Namely ${\mathbf{W}}_{\mathit{c}}$ is set to $[{\mathbf{e}}_{\mathit{c}\mathrm{,}1};{\mathbf{e}}_{\mathit{c}\mathrm{,}2};\dots ;{\mathbf{e}}_{\mathit{c}\mathrm{,}|{\mathcal{E}}_{\mathit{C}}|}]$ , ${\mathbf{W}}_{\mathit{s}}$ is set to $[{\mathbf{e}}_{\mathit{s}\mathrm{,}1};{\mathbf{e}}_{\mathit{s}\mathrm{,}2};\dots ;{\mathbf{e}}_{\mathit{s}\mathrm{,}|{\mathcal{E}}_{\mathit{S}}|}]$ , where $[·;·]$ denotes stacking vectors, ${\mathbf{e}}_{\mathit{c}\mathrm{,}\mathit{i}}$ denotes the embedding of the ith chemical in the chemicals ${\mathrm{KG}}_{\mathit{C}}$ , ${\mathbf{e}}_{\mathit{s}\mathrm{,}\mathit{i}}$ denotes the embedding of the ith species in the species ${\mathrm{KG}}_{\mathit{S}}$ .

In the experiments we refer to these models as Simple PT ${\mathrm{KGE}}_{\mathit{C}}$ - ${\mathrm{KGE}}_{\mathit{S}}$ and Complex PT ${\mathrm{KGE}}_{\mathit{C}}$ - ${\mathrm{KGE}}_{\mathit{S}}$ , depending on the selected MLP setting, where PT stands for pre-trained, and ${\mathrm{KGE}}_{\mathit{C}}$ and ${\mathrm{KGE}}_{\mathit{S}}$ are the KGE models used for the chemicals KG and the species KG, respectively (e.g., Complex PT DistMult-HAKE). For simplicity, we also refer to these models as PT-based models.

6.3. Fine-tuning optimization model

This model improves upon the pre-trained KG embeddings with fine-tuning based on the effect prediction data. This is done by simultaneously training the (selected) KGE models and the MLP-based baseline model. Such that the ${\mathbf{W}}_{\mathit{C}}$ and ${\mathbf{W}}_{\mathit{S}}$ , and the MLP weights ( ${\mathbf{W}}_{\mathit{x}}$ and ${\mathbf{b}}_{\mathit{x}}$ in Equations (10), (11), (14) and (15)) are optimized simultaneously. Note that we initialize the KGE models with the previously pre-trained embeddings.

Fig. 5.

Fine-tuning optimization model. In addition to variables described in Figs 4a and 4b, ${\mathit{t}}_{\mathit{C}}=\mathrm{(}{\mathit{sb}}_{\mathit{C}}\mathrm{,}{\mathit{p}}_{\mathit{C}}\mathrm{,}{\mathit{ob}}_{\mathit{C}}\mathrm{)}\in {\mathrm{KG}}_{\mathit{C}}\cup {\overline{\mathrm{KG}}}_{\mathit{C}}$ , ${\mathit{t}}_{\mathit{S}}=\mathrm{(}{\mathit{sb}}_{\mathit{S}}\mathrm{,}{\mathit{p}}_{\mathit{S}}\mathrm{,}{\mathit{ob}}_{\mathit{S}}\mathrm{)}\in {\mathrm{KG}}_{\mathit{S}}\cup {\overline{\mathrm{KG}}}_{\mathit{S}}$ . Entity lookups transform an entity into a vector (see Equation (16)). ${\mathit{SF}}_{{\mathit{KGE}}_{\mathit{C}}}$ and ${\mathit{SF}}_{{\mathit{KGE}}_{\mathit{S}}}$ are the triple scoring functions implemented by the selected KGE model (see the Appendix). ${\mathit{SF}}_{{\mathit{t}}_{\mathit{C}}}$ and ${\mathit{SF}}_{{\mathit{t}}_{\mathit{S}}}$ are the scores for a chemicals and species triple, respectively. ${\mathit{x}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}$ is the prediction input and ${\mathit{y}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}$ is described in Equation (8). ${\mathit{l}}_{{\mathit{t}}_{\mathit{C}}}$ and ${\mathit{l}}_{{\mathit{t}}_{\mathit{S}}}$ are the triple labels (i.e., True or False). $\mathit{BCE}$ is the binary cross-entropy loss function (from Equation (18)). The summation of the losses is described in Equation (17), that is the loss used by the optimizer to apply changes to model weights.

The model architecture is shown in Fig. 5 and the overall loss to minimize is $$\begin{array}{}\text{(17)}& \mathit{L}={\mathit{\alpha }}_{\mathit{C}}{\mathit{L}}_{{\mathit{KGE}}_{\mathit{C}}}+{\mathit{\alpha }}_{\mathit{S}}{\mathit{L}}_{{\mathit{KGE}}_{\mathit{S}}}+{\mathit{\alpha }}_{\mathrm{MLP}}{\mathit{L}}_{\mathrm{MLP}}\end{array}$$ where ${\mathit{L}}_{{\mathit{KGE}}_{\mathit{C}}}$ and ${\mathit{L}}_{{\mathit{KGE}}_{\mathit{S}}}$ respectively denote the loss of the chemical ${\mathrm{KG}}_{\mathit{C}}$ and the species ${\mathrm{KG}}_{\mathit{S}}$ when a specific KGE model is used,30

³⁰

Appendix A.5 introduces the used loss-functions in this work. The selection of the loss function for a KGE model will be via a hyper-parameter.

${\mathit{\alpha }}_{\mathit{C}}$ and ${\mathit{\alpha }}_{\mathit{S}}$ denote their weights respectively, ${\mathit{L}}_{\mathrm{MLP}}$ and ${\mathit{\alpha }}_{\mathrm{MLP}}$ denote the loss of the MLP and its weight. Specifically, we use binary cross-entropy (BCE) as the loss for the classification. ${\mathit{L}}_{\mathrm{MLP}}$ is calculated as $$\begin{array}{}\text{(18)}& {\mathit{L}}_{\mathrm{MLP}}=-{\displaystyle \frac{1}{\mathit{N}}}{\displaystyle \sum }_{\mathit{i}}^{\mathit{N}}{\mathit{y}}_{\mathit{i}}log\mathrm{(}{\hat{\mathit{y}}}_{\mathit{i}}\mathrm{)}+\mathrm{(}1-{\mathit{y}}_{\mathit{i}}\mathrm{)}log\mathrm{(}1-{\hat{\mathit{y}}}_{\mathit{i}}\mathrm{)}\end{array}$$ where N denotes the size of training samples, ${\mathit{y}}_{\mathit{i}}$ and ${\hat{\mathit{y}}}_{\mathit{i}}$ denote the sample label and the MLP output, respectively (as in Equation (8)). With the overall loss, gradient-based learning algorithms such as Adam optimizer [39] can be adopted to jointly training the embeddings of both KGEs and the MLP.

Figure 5 shows the full simultaneous fine-tuning model and the optimization process. The initial state of the entity lookups is the pre-trained embeddings. The full training procedure is summarised as follows:

Select N triples from ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ , where N is the length of the effects training set.31

³¹

Section 7.1 describes how the known effect data extracted from ECOTOX is split into training, validation and test sets.

Generate negative knowledge graph triples (see Appendix A.5 for details) from the extracted subsets of triples from ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ , these negative KGs triples are referred to as ${\overline{\mathrm{KG}}}_{\mathit{C}}$ and ${\overline{\mathrm{KG}}}_{\mathit{S}}$ .

Feed-forward the input through the model and calculate loss for each model component and combine according the loss weights.

Optimize the KG entity and relation embeddings, and the MLP layers.

These steps are repeated until the loss (only ${\mathit{L}}_{\mathrm{MLP}}$ ) over the validation set stops improving.

In the experiments we refer to these models as Simple FT ${\mathrm{KGE}}_{\mathit{C}}$ - ${\mathrm{KGE}}_{\mathit{S}}$ and Complex FT ${\mathrm{KGE}}_{\mathit{C}}$ - ${\mathrm{KGE}}_{\mathit{S}}$ , depending on the selected MLP setting, where FT stands for fine-tuning, and ${\mathrm{KGE}}_{\mathit{C}}$ and ${\mathrm{KGE}}_{\mathit{S}}$ are the KGE models used for the chemicals KG and the species KG, respectively (e.g., Simple FT HAKE-HAKE). For simplicity, we also refer to these models as FT-based models.

7. Results

7.1. Experimental setup

All models are implemented using Keras [16] and the model codes are available in our GitHub repository, alongside all data preparation and analysis scripts.32

³²

https://github.com/NIVA-Knowledge-Graph/KGs_and_Effect_Prediction_2020

7.1.1. Preparation of TERA for prediction

As shown earlier, TERA consists of three sub-KGs. These are the basis for the chemical effect prediction.33

³³

All data used to create TERA was downloaded on the 14th of May 2020.

We process the sub-KGs further to limit their size by removing irrelevant triples for prediction. This is necessary to scale up the training of the KGE models. The reduction of TERA’s sub-KGs is performed according to the following steps:

• Effect data. For prediction purposes, the effect data in ${\mathrm{KG}}_{\mathit{E}}$ is limited to four features, namely, chemical, species, chemical concentration, and effect. The chemical concentrations (κ, converted to mg/L) are log-normalized to remove the large discrepancy in scales. As mentioned, we separate the effects into two categories for simplicity, lethal and non-lethal effects. This reduces the possibility of ambiguity among the effects that does not cause death in the test species. We label lethal effects as 1 and non-lethal effects as 0

• ${\mathrm{KG}}_{\mathit{C}}$ . For each chemical in the effect data, we extract all triples connected to them using a directed crawl. This reduces the size of ${\mathrm{KG}}_{\mathit{C}}$ to a manageable size for the KGE models. Moreover, we do not deem triples not directly connected to the effect data relevant for the prediction task, and may introduce unnecessary noise. As mentioned before, PubChem contains similarities between chemicals based on chemical fingerprints, however, for our use-case it is unpractical to query them from the PubChem RDF data, therefore, we calculate similarity triples based on queried PubChem fingerprints. We use the same similarity threshold as PubChem, i.e., 0.9 [37].

• ${\mathrm{KG}}_{\mathit{S}}$ . The same steps as for ${\mathrm{KG}}_{\mathit{C}}$ are conducted for all species in the effect data.

A simple directed crawl over all predicates is sufficient to gather the interesting data in this setting as both ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ are primarily hierarchical and we start the crawls at the leaf nodes.

These steps reduce ${\mathrm{KG}}_{\mathit{C}}$ to 241,442 triples and ${\mathrm{KG}}_{\mathit{S}}$ to 59,673 triples. Some statistics of ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ , and the reduced fragments ${\mathrm{KG}}_{\mathit{C}}^{\prime }$ and ${\mathrm{KG}}_{\mathit{S}}^{\prime }$ , are given in Table 7 (Section 5.4). In the rest of the paper were refer to TERA’s reduced sub-KGs simply as ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ .

The transformation from TERA’s ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ to model input is done by first dropping literals, thereafter assigning each entity an unique integer identifier which corresponds to the index of a column vector in matrices ${\mathbf{W}}_{\mathit{c}}$ or ${\mathbf{W}}_{\mathit{s}}$ in Equation (16), depending on which sub-KG is transformed.34

³⁴

$\mathit{i}\in [0\mathrm{,}|{\mathcal{E}}_{\mathit{C}}|-1]$ for ${\mathrm{KG}}_{\mathit{C}}$ and $\mathit{i}\in [0\mathrm{,}|{\mathcal{E}}_{\mathit{S}}|-1]$ for ${\mathrm{KG}}_{\mathit{S}}$ .

Relations are treated similarly.

7.1.2. Sampling

We use four sampling strategies of the effect data to analyze how the proposed classification models behave by varying the data parts that are used for training and testing. Note that, we only consider effect data where the chemical and species have mappings to external sources (e.g., NCBI Taxonomy and Wikidata, cf. Section 5.2.2) so that there is additional contextual information that can be used by the KGE models. For each of the strategies, the validation and test sets contain unseen chemical-organism pairs with respect to the training set. The strategies, however, differ with respect to the individual organism and chemical as follows:

Random $70\mathrm{\%}\mathrm{/}15\mathrm{\%}\mathrm{/}15\mathrm{\%}$ training/validation/test split on the entire dataset (i.e., the chemicals and the organisms in the validation and test will most probably be known).

Training/validation/test split where there is no overlap between chemicals in the three sets (i.e., the chemicals in the validation and test sets are unknown). This resulted on a $77\mathrm{\%}\mathrm{/}14\mathrm{\%}\mathrm{/}9\mathrm{\%}$  split.

Training/validation/test split where there is no overlap between species in the three sets (i.e., the species in the validation and test sets are unknown). This resulted on a $77\mathrm{\%}\mathrm{/}14\mathrm{\%}\mathrm{/}9\mathrm{\%}$ split.

Training/validation/test split with no chemicals or species overlap in the three sets (i.e., both the chemicals and the organisms in the validation and test sets are unknown). This resulted on a $72\mathrm{\%}\mathrm{/}14\mathrm{\%}\mathrm{/}14\mathrm{\%}$ split.

Note that since we use the species and chemicals as groups to divide the data rather than the samples, the splits can vary. For strategies (i)–(iii) there is a total of 14,377 effect data samples while for strategy (iv) the total number samples is 5,621. As above, this discrepancy is down to the way we split the data. We do not split across samples, but across chemicals and species. For example, some chemicals are used on (close to) all species, therefore, these chemicals are discarded in the sampling strategy (iv), affecting the final number of samples.

There were originally 57,560 samples, however, this includes experiment duplicates, i.e., same chemical, species, and endpoint, with different chemical concentrations. This is down to large discrepancies in laboratory testing variance, therefore, we use the median concentration across the duplicates. The prior probability is approximately $0.16\mathrm{/}0.84$ (i.e., $\approx 16\mathrm{\%}$ of samples are labelled as non-lethal and $\approx 84\mathrm{\%}$ of samples are labelled as lethal) across all sampling methods. We solve this when training by randomly oversampling the minority class until the prior probabilities are $0.5\mathrm{/}0.5$ in the training set. In this case, the oversampling is performed by adding duplicates samples labelled as non-lethal. Oversampling is a well established technique used in many classification problems to remove bias during learning [11].

Table 8

Hyper-parameter choices for the models. Please refer to the Equations (9)–(15) in Section 6.1 for the prediction hyper-parameters

KGE hyper-parameters	Search space
Loss function	$\{{\mathit{L}}_{{\mathit{H}}_1}\mathrm{,}{\mathit{L}}_{{\mathit{H}}_2}\mathrm{,}{\mathit{L}}_{{\mathit{L}}_1}\mathrm{,}{\mathit{L}}_{{\mathit{L}}_2}\}$
Margin (only hinge loss)	$\{1\mathrm{,}2\mathrm{,}\dots \mathrm{,}10\}$
Bias (only geometric models)	$\{0\mathrm{,}1\mathrm{,}\dots \mathrm{,}20\}$
Embedding dimension	$\{100\mathrm{,}101\mathrm{,}\dots \mathrm{,}400\}$
Negative samples	$\{10\mathrm{,}11\mathrm{,}\dots \mathrm{,}100\}$

Prediction hyper-parameters	Search space
${\mathit{n}}_{\mathit{c}}$ (10), ${\mathit{n}}_{\mathit{s}}$ (11), ${\mathit{n}}_{\mathit{\kappa }}$ (12), n (14)	$\{0\mathrm{,}1\mathrm{,}2\mathrm{,}3\}$
# units (10), (11), (14)	$\{2^{\mathit{u}}\text{with}\mathit{u}\in \{4\mathrm{,}5\mathrm{,}\dots \mathrm{,}10\}\}$
# units (12)	$\{2^{\mathit{u}}\text{with}\mathit{u}\in \{2\mathrm{,}3\mathrm{,}4\mathrm{,}5\}\}$

Table 9

Best hyper-parameters for KGE models. The two values before and after / are for the embeddings of ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ , respectively

Model	Loss function	Margin	Bias	Embedding dimension	Negative samples
DistMult	${\mathit{L}}_{{\mathit{L}}_2}$  /  ${\mathit{L}}_{{\mathit{H}}_2}$	– / 2	–	143 / 383	28 / 43
ComplEx	${\mathit{L}}_{{\mathit{L}}_2}$  /  ${\mathit{L}}_{{\mathit{H}}_2}$	– / 4	–	163 / 372	27 / 42
HolE	${\mathit{L}}_{{\mathit{H}}_2}$  /  ${\mathit{L}}_{{\mathit{L}}_2}$	6 / –	–	188 / 376	30 / 100
TransE	${\mathit{L}}_{{\mathit{H}}_2}$  /  ${\mathit{L}}_{{\mathit{H}}_1}$	4 / 7	14 / 20	226 / 196	23 / 57
RotatE	${\mathit{L}}_{{\mathit{H}}_2}$  /  ${\mathit{L}}_{{\mathit{H}}_2}$	5 / 2	16 / 6	271 / 398	75 / 22
pRotatE	${\mathit{L}}_{{\mathit{L}}_2}$  /  ${\mathit{L}}_{{\mathit{L}}_2}$	– / –	14 / 16	164 / 210	34 / 82
HAKE	${\mathit{L}}_{{\mathit{L}}_2}$  /  ${\mathit{L}}_{{\mathit{L}}_2}$	– / –	12 / 10	108 / 359	56 / 13
ConvKB	${\mathit{L}}_{{\mathit{L}}_2}$  /  ${\mathit{L}}_{{\mathit{H}}_2}$	– / 5	–	248 / 276	18 / 90
ConvE	${\mathit{L}}_{{\mathit{H}}_1}$  /  ${\mathit{L}}_{{\mathit{H}}_1}$	7 / 3	–	228 / 196	68 / 40

7.1.3. Hyper-parameters

To optimize the hyper-parameters for the KGE and classification models we use random search over the parameter ranges. We conduct 20 trials per model. Tables 8 and 9 contain the best hyper-parameters and can be used to reproduce the top performing models.

To find the best hyper-parameters for the KGE models, we use the loss as a proxy for performance, normalized by the initial loss, $\mathit{R}{\mathit{L}}_{\mathrm{ep}}={\mathit{L}}_{\mathrm{ep}}\mathrm{/}{\mathit{L}}_0$ , where ${\mathit{L}}_{\mathrm{ep}}$ is the training loss at epoch $\mathrm{ep}$ , ${\mathit{L}}_0$ is the loss with the initial weights.

We use validation loss to select the best hyper-parameter setting for the classification models presented in Section 6. The best prediction models are refitted and evaluated 10 times to reduce the influence of initial conditions on the metrics. The average and standard deviation of the metrics are presented in Section 7.2.

The hyper-parameter ranges for the KGE models are shown in Table 8 based on common values used in the literature. We conduct 20 trials of random hyper-parameters choices and validate over the validation data. In Table 9 we show the best hyper-parameters.

Table 10

Number of units in the hidden layers in the (complex) one-hot model and the top-1 prediction models with pre-trained KG embeddings. The same parameters are used for the fine-tuning models. Organized as follows: $\mathrm{(}|{\mathbf{b}}_{\mathit{c}}^1|\mathrm{,}\dots \mathrm{,}|{\mathbf{b}}_{\mathit{c}}^{{\mathit{n}}_{\mathit{c}}}|\mathrm{)}\mathrm{/}\mathrm{(}|{\mathbf{b}}_{\mathit{s}}^1|\mathrm{,}\dots \mathrm{,}|{\mathbf{b}}_{\mathit{s}}^{{\mathit{n}}_{\mathit{s}}}|\mathrm{)}\mathrm{/}\mathrm{(}|{\mathbf{b}}_{\mathit{\kappa }}^1|\mathrm{,}\dots \mathrm{,}|{\mathbf{b}}_{\mathit{\kappa }}^{{\mathit{n}}_{\mathit{\kappa }}}|\mathrm{)}\mathrm{/}\mathrm{(}|{\mathbf{b}}^1|\mathrm{,}\dots \mathrm{,}|{\mathbf{b}}^{\mathit{n}}|\mathrm{)}$ as in Equations (10), (11), (12), and (14)). − denotes no hidden layers. e.g., $\mathrm{(}128\mathrm{)}\mathrm{/}\mathrm{(}256\mathrm{)}\mathrm{/}\mathrm{(}8\mathrm{,}8\mathrm{)}\mathrm{/}-$ denotes ${\mathit{n}}_{\mathit{c}}=1$ , ${\mathit{n}}_{\mathit{s}}=1$ , ${\mathit{n}}_{\mathit{\kappa }}=2$ , $\mathit{n}=0$ and $|\mathbf{b}{|}_{\mathit{c}}^1=128$ , $|{\mathbf{b}}_{\mathit{s}}^1|=256$ , $|{\mathbf{b}}_{\mathit{\kappa }}^1|=8$ and $|{\mathbf{b}}_{\mathit{\kappa }}^2|=8$

Model	Sampling	# units
Complex one-hot	(i)	$\mathrm{(}128\mathrm{)}\mathrm{/}\mathrm{(}128\mathrm{)}\mathrm{/}\text{–}\mathrm{/}\text{–}$
	(ii)	$\mathrm{(}128\mathrm{)}\mathrm{/}\mathrm{(}256\mathrm{)}\mathrm{/}\mathrm{(}8\mathrm{,}8\mathrm{)}\mathrm{/}\text{–}$
	(iii)	$\mathrm{(}256\mathrm{,}128\mathrm{)}\mathrm{/}\mathrm{(}128\mathrm{)}\mathrm{/}\mathrm{(}4\mathrm{,}4\mathrm{,}4\mathrm{)}\mathrm{/}\text{–}$
	(iv)	$\mathrm{(}256\mathrm{,}256\mathrm{)}\mathrm{/}\mathrm{(}128\mathrm{)}\mathrm{/}\mathrm{(}8\mathrm{,}8\mathrm{)}\mathrm{/}\mathrm{(}128\mathrm{)}$
Complex PT DistMult-HAKE (top-1 in (i))	(i)	$\mathrm{(}256\mathrm{,}256\mathrm{)}\mathrm{/}\mathrm{(}256\mathrm{)}\mathrm{/}\mathrm{(}16\mathrm{,}4\mathrm{)}\mathrm{/}\mathrm{(}512\mathrm{,}64\mathrm{)}$
Complex PT HolE-ConvKB (top-1 in (ii))	(ii)	$\mathrm{(}512\mathrm{,}128\mathrm{,}128\mathrm{)}\mathrm{/}\mathrm{(}512\mathrm{)}\mathrm{/}\text{–}\mathrm{/}\mathrm{(}64\mathrm{)}$
Complex PT HAKE-DistMult (top-1 in (iii), (iv))	(iii)	$\mathrm{(}64\mathrm{)}\mathrm{/}\mathrm{(}512\mathrm{)}\mathrm{/}\mathrm{(}16\mathrm{,}32\mathrm{)}\mathrm{/}\mathrm{(}16\mathrm{)}$
	(iv)	$\mathrm{(}128\mathrm{)}\mathrm{/}\text{–}\mathrm{/}\mathrm{(}4\mathrm{,}8\mathrm{,}8\mathrm{)}\mathrm{/}\mathrm{(}256\mathrm{,}128\mathrm{)}$

We can see in Table 9 that the decomposition models have similar hyper-parameters for ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ . As shown in Section 5.4, the major difference between ${\mathrm{KG}}_{\mathit{C}}$ and ${\mathrm{KG}}_{\mathit{S}}$ is the relational density. Therefore, it is reasonable to believe that a lower relational density KG requires more parameters to have an equivalent representation in the embedding space. We can get the same observation for the geometric models except for TransE, where the embedding dimensions are similar. ConvE is more efficient in embedding dimension than ConvKB, however, since ConvE is slightly more complex than ConvKB this is expected. The difference in negative samples could be down to our implementation of ConvE, which varies from the original. Our implementation of all models relies on 1-to-1 scoring of triples, while the implementation of ConvE originally used 1-to- $|\mathcal{E}|$ scoring, where $|\mathcal{E}|$ is the number of entities in the KG [19].

The fine-tuning optimization model (Section 6.3), in order to save on intensive computation, reuses the same hyper-parameters found for the KGE models. Depending on the optimizer choice, the choice of loss weights, ${\mathit{\alpha }}_{\mathit{C}}$ , ${\mathit{\alpha }}_{\mathit{S}}$ , and ${\mathit{\alpha }}_{\mathrm{MLP}}$ , is important. However, our optimizer choice has dynamic learning rates per variable, and therefore, will adapt regardless of the loss weights and we can set ${\mathit{\alpha }}_{\mathit{C}}={\mathit{\alpha }}_{\mathit{S}}={\mathit{\alpha }}_{\mathrm{MLP}}=1$ . Had we used, e.g., stochastic gradient descent, these variables would needed to be tuned.

7.1.4. Initialization of the fine-tuning optimization models

As presented in Section 6.3, we simultaneously train the KGE models and the MLP-based baseline model. This is done by initializing the model with (i) the weights learned in the correspondent baseline model with pre-trained embeddings, and (ii) the KG embeddings learned with the respective KGE models. For example, the Complex FT DistMult-HAKE model is initialized with the learned weights with the Complex PT DistMult-HAKE model and the pre-trained KG embeddings using DistMult and HAKE models. Then the model is further trained with a small learning rate. We found that reducing the learning rate by a factor of 100 worked well. Using this learning rate we optimize the model until convergence.

7.1.5. Simple and complex settings

As presented in Section 6.1, we use two settings in our classification models: simple and complex. This will help us isolate the effects of the KG embeddings versus the power of the MLP model. The simple setting uses no branching layers, i.e., ${\mathit{n}}_{\mathit{C}}={\mathit{n}}_{\mathit{S}}={\mathit{n}}_{\mathit{\kappa }}=0$ and $\mathit{n}=1$ as in Equations (10), (11), (12) and (14) with 128 units in the hidden dense layer. For the complex models we use random search (20 trials) to find the optimal number of layers and units out of the ranges shown in Table 8. The optimal choices for the top performing models (using one-hot and pre-trained embeddings) are shown in Table 10.

Looking at the increasing complexity of the layer configuration of the one-hot models in Table 10 we can see a correlation from the simplest sampling strategy (i.e., (i)) through the most challenging one (i.e., (iv)). The same can be seen for PT HAKE-DisMult from strategy (iii) to (iv), where the number of layers increase. Overall we can see that the layer configurations of the chemical branch is more complex than for the species branch. This indicates that the KGE models are better at representing ${\mathrm{KG}}_{\mathit{S}}$ than  ${\mathrm{KG}}_{\mathit{C}}$ .

7.2. Prediction results

In this section we present a summary of the conducted chemical effect prediction evaluation. Complete results are available at the project repository.35

³⁵

https://github.com/NIVA-Knowledge-Graph/KGs_and_Effect_Prediction_2020

The default decision threshold is set to 0.5. That is, if a model predicts $\hat{\mathit{y}}\mathrm{>}0.5$ for an input ${\mathit{x}}_{\mathit{c}\mathrm{,}\mathit{s}\mathrm{,}\mathit{\kappa }}$ then the chemical c is considered lethal to s at a concentration κ.36

³⁶

We set the decision threshold $\hat{\mathit{y}}\mathrm{>}0.5$ since the model output bias (cf. Equation (15)) will be (close to) 0.5 after training. Recall that we have oversampled the classes to reach a $0.5\mathrm{/}0.5$ prior probability during training (cf. Section 7.1.2).

We use several metrics to compare the different prediction models. These are Sensitivity (i.e., recall), Specificity, and Youden’s index ( $\mathrm{YI}$ ) [85]. Precision and F-score were also considered as metrics. However, they were not representative for the performance with respect to non-harmful chemicals. This is attributed to the larger number of positive samples (i.e., harmful chemicals) than negative samples (i.e., non-harmful chemicals) in the test data.

Sensitivity and Specificity are defined as $$\begin{array}{}\text{(19)}& \text{Sensitivity}={\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}\mathrm{,}\text{(20)}& \text{Specificity}={\displaystyle \frac{\mathrm{TN}}{\mathrm{FP}+\mathrm{TN}}}\mathrm{,}\end{array}$$ where TP, FN, TN, and FP are true positives, false negatives, true negatives and false positives, respectively. YI is defined as $$\begin{array}{}\text{(21)}& \mathrm{YI}=\text{Sensitivity}+\text{Specificity}-1.\end{array}$$ We also present the maximized Youden’s index ( ${\mathrm{YI}}_{max}$ ), this is defined as $$\begin{array}{}\text{(22)}& {\mathrm{YI}}_{max}=\underset{\mathit{\tau }}{max}\text{Sensitivity}+\text{Specificity}-1\mathrm{,}\end{array}$$ i.e., we maximize Youden’s index based on the decision threshold (τ), we call this optimal threshold ${\mathit{\tau }}_{max}$ . This metric is equivalent to the maximum of the Receiver operating characteristic (ROC) curve over a random model and can be used to select the optimal decision threshold in a production environment (based on validation data). We do not present ROC (or area under ROC, AUC) as a metric as it correlates ( $\mathrm{>}0.99$ ) with ${\mathrm{YI}}_{max}$ in our case.

In our setting, sensitivity is a measure on how well the models identify harmful chemicals while specificity measures models’ ability to identify non-harmful chemicals. Youden’s index is used to capture the usefulness of a diagnostic test (or in our case, a toxicity test). A useless test will have $\mathrm{YI}=0$ while with $\mathrm{YI}\mathrm{>}0$ a test is useful. $\mathrm{YI}$ is also thought of as how well informed a decision might be. Note that, $\mathrm{YI}$ can be less than 0, but this is solved by swapping labeled classes. Similarly to how negative correlation is still useful.

Table 11

Prediction results (mean and standard deviation over 10 runs) for sampling strategy (i). Bold denotes best mean result and underline denotes within one standard deviation of best result. PT prefix denotes pre-trained and FT denotes fine-tuning. Simple denotes ${\mathit{n}}_{\mathit{C}}={\mathit{n}}_{\mathit{S}}={\mathit{n}}_{\mathit{\kappa }}=0$ and $\mathit{n}=1$ while in complex, ${\mathit{n}}_{\mathit{C}}$ , ${\mathit{n}}_{\mathit{S}}$ , ${\mathit{n}}_{\mathit{\kappa }}$ and n are hyper-parameters in Equations (10), (11), (12) and (14)

Model	Sensitivity	Specificity	YI	${\mathrm{YI}}_{max}$	${\mathit{\tau }}_{max}$
Simple one-hot	$0.939\pm 0.009$	$0.657\pm 0.018$	$0.595\pm 0.015$	$0.666\pm 0.011$	$0.809\pm 0.049$
Simple PT HAKE-HAKE	$0.912\pm 0.006$	$0.773\pm 0.018$	$0.685\pm 0.016$	$0.719\pm 0.012$	$0.707\pm 0.044$
Simple PT pRotatE-HAKE	$0.934\pm 0.005$	$0.749\pm 0.044$	$0.683\pm 0.04$	$0.718\pm 0.02$	$0.665\pm 0.082$
Simple PT ConvE-HAKE	$0.937\pm 0.006$	$0.738\pm 0.006$	$0.674\pm 0.004$	$0.724\pm 0.007$	$0.721\pm 0.054$
Simple PT pRotatE-ConvE	$0.924\pm 0.029$	$0.436\pm 0.155$	$0.36\pm 0.182$	$0.469\pm 0.196$	$0.784\pm 0.052$
Simple PT RotatE-ConvE	$\mathbf{0.997}\pm \mathbf{0.003}$	$0.024\pm 0.035$	$0.021\pm 0.035$	$0.195\pm 0.111$	$0.812\pm 0.086$
Simple FT HAKE-HAKE	$0.921\pm 0.005$	$0.814\pm 0.009$	$\mathbf{0.734}\pm \mathbf{0.006}$	$0.743\pm 0.007$	$0.547\pm 0.074$
Simple FT pRotatE-HAKE	$0.92\pm 0.005$	$0.808\pm 0.013$	$\underset{\_}{0.728\pm 0.011}$	$0.738\pm 0.007$	$0.56\pm 0.107$
Simple FT ConvE-HAKE	$0.942\pm 0.003$	$0.733\pm 0.019$	$0.675\pm 0.019$	$0.729\pm 0.007$	$0.864\pm 0.053$
Simple FT pRotatE-ConvE	$0.949\pm 0.003$	$0.766\pm 0.017$	$0.715\pm 0.016$	$\mathbf{0.765}\pm \mathbf{0.006}$	$0.842\pm 0.064$
Simple FT RotatE-ConvE	$0.928\pm 0.015$	$0.797\pm 0.036$	$\underset{\_}{0.726\pm 0.022}$	$\underset{\_}{0.761\pm 0.01}$	$0.722\pm 0.069$
Complex one-hot	$0.937\pm 0.004$	$0.748\pm 0.016$	$0.685\pm 0.015$	$0.728\pm 0.009$	$0.769\pm 0.094$
Complex PT DistMult-HAKE	$0.895\pm 0.008$	$0.817\pm 0.008$	$0.713\pm 0.007$	$0.723\pm 0.008$	$0.456\pm 0.088$
Complex PT HAKE-ConvKB	$0.927\pm 0.006$	$0.784\pm 0.017$	$0.711\pm 0.013$	$0.739\pm 0.009$	$0.686\pm 0.109$
Complex PT HolE-ConvKB	$0.932\pm 0.013$	$0.779\pm 0.024$	$0.711\pm 0.013$	$0.729\pm 0.009$	$0.676\pm 0.104$
Complex PT ComplEx-DistMult	$0.96\pm 0.006$	$0.584\pm 0.04$	$0.543\pm 0.039$	$0.664\pm 0.024$	$0.838\pm 0.048$
Complex PT HolE-pRotatE	$\underset{\_}{0.996\pm 0.006}$	$0.011\pm 0.02$	$0.006\pm 0.014$	$0.182\pm 0.041$	$0.804\pm 0.071$
Complex FT DistMult-HAKE	$0.903\pm 0.009$	$0.816\pm 0.015$	$0.719\pm 0.008$	$0.729\pm 0.005$	$0.597\pm 0.098$
Complex FT HAKE-ConvKB	$0.935\pm 0.006$	$0.791\pm 0.021$	$\underset{\_}{0.726\pm 0.018}$	$0.754\pm 0.008$	$0.776\pm 0.109$
Complex FT HolE-ConvKB	$0.895\pm 0.01$	$\mathbf{0.835}\pm \mathbf{0.016}$	$\underset{\_}{0.73\pm 0.01}$	$0.739\pm 0.011$	$0.61\pm 0.123$
Complex FT ComplEx-DistMult	$0.927\pm 0.005$	$0.78\pm 0.018$	$0.707\pm 0.016$	$0.742\pm 0.011$	$0.797\pm 0.093$
Complex FT HolE-pRotatE	$0.913\pm 0.008$	$0.795\pm 0.017$	$0.708\pm 0.012$	$0.734\pm 0.008$	$0.777\pm 0.049$

Table 12

Prediction results for sampling strategy (ii). Same notation as Table 11

Model	Sensitivity	Specificity	YI	${\mathrm{YI}}_{max}$	${\mathit{\tau }}_{max}$
Simple one-hot	$0.88\pm 0.022$	$0.628\pm 0.048$	$0.508\pm 0.057$	$0.556\pm 0.051$	$0.713\pm 0.13$
Simple PT HAKE-ConvKB	$0.926\pm 0.007$	$0.823\pm 0.016$	$0.748\pm 0.017$	$0.775\pm 0.013$	$0.623\pm 0.064$
Simple PT HAKE-HAKE	$0.908\pm 0.007$	$0.829\pm 0.014$	$0.738\pm 0.012$	$0.759\pm 0.01$	$0.613\pm 0.132$
Simple PT pRotatE-HAKE	$0.924\pm 0.003$	$0.802\pm 0.009$	$0.726\pm 0.008$	$0.76\pm 0.006$	$0.79\pm 0.084$
Simple PT RotatE-ConvKB	$0.972\pm 0.021$	$0.42\pm 0.255$	$0.392\pm 0.236$	$0.62\pm 0.111$	$0.814\pm 0.06$
Simple PT RotatE-ConvE	$\mathbf{0.997}\pm \mathbf{0.004}$	$0.021\pm 0.057$	$0.018\pm 0.054$	$0.22\pm 0.088$	$0.824\pm 0.095$
Simple FT HAKE-ConvKB	$0.909\pm 0.003$	$0.883\pm 0.006$	$\mathbf{0.792}\pm \mathbf{0.006}$	$0.803\pm 0.004$	$0.556\pm 0.138$
Simple FT HAKE-HAKE	$0.897\pm 0.007$	$0.86\pm 0.01$	$0.757\pm 0.012$	$0.769\pm 0.006$	$0.61\pm 0.134$
Simple FT pRotatE-HAKE	$0.905\pm 0.004$	$0.859\pm 0.012$	$0.764\pm 0.012$	$0.775\pm 0.011$	$0.544\pm 0.099$
Simple FT RotatE-ConvKB	$0.93\pm 0.007$	$0.853\pm 0.013$	$\underset{\_}{0.784\pm 0.008}$	$\mathbf{0.81}\pm \mathbf{0.008}$	$0.732\pm 0.119$
Simple FT RotatE-ConvE	$0.912\pm 0.02$	$0.821\pm 0.028$	$0.733\pm 0.01$	$0.753\pm 0.005$	$0.735\pm 0.17$
Complex one-hot	$0.875\pm 0.014$	$0.859\pm 0.015$	$0.734\pm 0.012$	$0.749\pm 0.009$	$0.448\pm 0.2$
Complex PT HolE-ConvKB	$0.894\pm 0.006$	$0.889\pm 0.014$	$\underset{\_}{0.783\pm 0.014}$	$0.793\pm 0.01$	$0.489\pm 0.035$
Complex PT pRotatE-ConvKB	$0.901\pm 0.012$	$0.875\pm 0.027$	$\underset{\_}{0.776\pm 0.024}$	$0.79\pm 0.018$	$0.592\pm 0.081$
Complex PT TransE-ConvKB	$0.906\pm 0.008$	$0.868\pm 0.021$	$\underset{\_}{0.774\pm 0.019}$	$0.787\pm 0.012$	$0.588\pm 0.112$
Complex PT ComplEx-ConvE	$0.928\pm 0.006$	$0.768\pm 0.015$	$0.696\pm 0.015$	$0.731\pm 0.008$	$0.689\pm 0.095$
Complex PT ConvKB-pRotatE	$\underset{\_}{0.995\pm 0.005}$	$0.011\pm 0.012$	$0.007\pm 0.008$	$0.265\pm 0.054$	$0.77\pm 0.089$
Complex FT HolE-ConvKB	$0.871\pm 0.007$	$0.906\pm 0.007$	$0.778\pm 0.007$	$0.791\pm 0.005$	$0.441\pm 0.07$
Complex FT pRotatE-ConvKB	$0.869\pm 0.008$	$\mathbf{0.914}\pm \mathbf{0.011}$	$0.783\pm 0.007$	$0.794\pm 0.006$	$0.483\pm 0.083$
Complex FT TransE-ConvKB	$0.878\pm 0.008$	$0.895\pm 0.011$	$0.772\pm 0.008$	$0.792\pm 0.006$	$0.511\pm 0.133$
Complex FT ComplEx-ConvE	$0.916\pm 0.009$	$0.83\pm 0.021$	$0.746\pm 0.016$	$0.76\pm 0.011$	$0.596\pm 0.151$
Complex FT ConvKB-pRotatE	$0.9\pm 0.013$	$0.794\pm 0.026$	$0.694\pm 0.018$	$0.723\pm 0.014$	$0.785\pm 0.111$

Table 13

Prediction results for sampling strategy (iii). Same notation as Table 11

Model	Sensitivity	Specificity	YI	${\mathrm{YI}}_{max}$	${\mathit{\tau }}_{max}$
Simple one-hot	$0.822\pm 0.058$	$0.439\pm 0.054$	$0.261\pm 0.058$	$0.31\pm 0.047$	$0.597\pm 0.182$
Simple PT ConvKB-DistMult	$0.966\pm 0.007$	$0.626\pm 0.047$	$\underset{\_}{0.591\pm 0.045}$	$\underset{\_}{0.623\pm 0.049}$	$0.67\pm 0.058$
Simple PT HAKE-DistMult	$0.958\pm 0.023$	$0.628\pm 0.026$	$0.586\pm 0.033$	$\underset{\_}{0.626\pm 0.045}$	$0.613\pm 0.092$
Simple PT ConvKB-TransE	$0.969\pm 0.009$	$0.614\pm 0.048$	$\underset{\_}{0.583\pm 0.04}$	$\underset{\_}{0.642\pm 0.01}$	$0.643\pm 0.059$
Simple PT ConvE-RotatE	$0.934\pm 0.055$	$0.276\pm 0.026$	$0.209\pm 0.043$	$0.273\pm 0.071$	$0.596\pm 0.13$
Simple PT HolE-HAKE	$0.88\pm 0.089$	$0.115\pm 0.083$	$-0.005\pm 0.075$	$0.077\pm 0.057$	$0.783\pm 0.18$
Simple FT ConvKB-DistMult	$0.947\pm 0.014$	$0.667\pm 0.02$	$\underset{\_}{0.614\pm 0.013}$	$\underset{\_}{0.645\pm 0.011}$	$0.736\pm 0.087$
Simple FT HAKE-DistMult	$0.947\pm 0.012$	$\underset{\_}{0.662\pm 0.035}$	$\underset{\_}{0.609\pm 0.031}$	$\underset{\_}{0.634\pm 0.026}$	$0.701\pm 0.132$
Simple FT ConvKB-TransE	$0.934\pm 0.009$	$\underset{\_}{0.68\pm 0.018}$	$\underset{\_}{0.615\pm 0.014}$	$\underset{\_}{0.642\pm 0.015}$	$0.687\pm 0.065$
Simple FT ConvE-RotatE	$0.915\pm 0.013$	$0.454\pm 0.028$	$0.369\pm 0.027$	$0.402\pm 0.028$	$0.658\pm 0.083$
Simple FT HolE-HAKE	$0.931\pm 0.009$	$0.118\pm 0.036$	$0.049\pm 0.038$	$0.171\pm 0.038$	$0.882\pm 0.127$
Complex one-hot	$0.796\pm 0.028$	$0.571\pm 0.041$	$0.367\pm 0.054$	$0.398\pm 0.043$	$0.526\pm 0.076$
Complex PT HAKE-DistMult	$0.969\pm 0.016$	$\underset{\_}{0.642\pm 0.044}$	$\underset{\_}{0.61\pm 0.034}$	$\underset{\_}{0.643\pm 0.026}$	$0.675\pm 0.105$
Complex PT pRotatE-ComplEx	$0.929\pm 0.024$	$\underset{\_}{0.668\pm 0.048}$	$\underset{\_}{0.597\pm 0.048}$	$\underset{\_}{0.62\pm 0.046}$	$0.526\pm 0.145$
Complex PT ConvKB-DistMult	$0.965\pm 0.013$	$\underset{\_}{0.631\pm 0.078}$	$\underset{\_}{0.597\pm 0.07}$	$\underset{\_}{0.627\pm 0.039}$	$0.597\pm 0.149$
Complex PT ComplEx-HolE	$\mathbf{0.991}\pm \mathbf{0.01}$	$0.237\pm 0.106$	$0.228\pm 0.098$	$0.45\pm 0.028$	$0.721\pm 0.047$
Complex PT ComplEx-HAKE	$0.9\pm 0.055$	$0.097\pm 0.047$	$-0.003\pm 0.064$	$0.133\pm 0.081$	$0.696\pm 0.22$
Complex FT HAKE-DistMult	$0.932\pm 0.011$	$\mathbf{0.69}\pm \mathbf{0.024}$	$\mathbf{0.622}\pm \mathbf{0.023}$	$\mathbf{0.652}\pm \mathbf{0.022}$	$0.706\pm 0.134$
Complex FT pRotatE-ComplEx	$0.931\pm 0.025$	$\underset{\_}{0.672\pm 0.042}$	$\underset{\_}{0.602\pm 0.045}$	$\underset{\_}{0.631\pm 0.037}$	$0.627\pm 0.157$
Complex FT ConvKB-DistMult	$0.953\pm 0.008$	$0.642\pm 0.027$	$\underset{\_}{0.596\pm 0.027}$	$\underset{\_}{0.625\pm 0.028}$	$0.753\pm 0.138$
Complex FT ComplEx-HolE	$0.898\pm 0.035$	$0.591\pm 0.064$	$0.489\pm 0.042$	$0.521\pm 0.027$	$0.612\pm 0.156$
Complex FT ComplEx-HAKE	$0.88\pm 0.032$	$0.255\pm 0.026$	$0.135\pm 0.034$	$0.204\pm 0.06$	$0.775\pm 0.268$

Table 14

Prediction results sampling strategy (iv). Same notation as Table 11

Model	Sensitivity	Specificity	YI	${\mathrm{YI}}_{max}$	${\mathit{\tau }}_{max}$
Simple one-hot	$0.612\pm 0.096$	$0.421\pm 0.107$	$0.033\pm 0.14$	$0.113\pm 0.076$	$0.555\pm 0.306$
Simple PT HAKE-ComplEx	$\underset{\_}{0.971\pm 0.011}$	$0.361\pm 0.065$	$0.332\pm 0.056$	$\underset{\_}{0.546\pm 0.031}$	$0.89\pm 0.042$
Simple PT pRotatE-ComplEx	$\mathbf{0.972}\pm \mathbf{0.008}$	$0.36\pm 0.079$	$0.332\pm 0.074$	$\underset{\_}{0.527\pm 0.045}$	$0.852\pm 0.04$
Simple PT HolE-ComplEx	$\underset{\_}{0.965\pm 0.032}$	$0.363\pm 0.068$	$0.328\pm 0.063$	$\underset{\_}{0.549\pm 0.075}$	$0.856\pm 0.077$
Simple PT pRotatE-RotatE	$0.917\pm 0.01$	$0.168\pm 0.016$	$0.084\pm 0.013$	$0.151\pm 0.021$	$0.779\pm 0.182$
Simple PT HAKE-HAKE	$0.8\pm 0.095$	$0.128\pm 0.066$	$-0.072\pm 0.07$	$0.033\pm 0.027$	$0.736\pm 0.321$
Simple FT HAKE-ComplEx	$\underset{\_}{0.963\pm 0.01}$	$0.423\pm 0.102$	$\underset{\_}{0.386\pm 0.096}$	$\underset{\_}{0.57\pm 0.03}$	$0.875\pm 0.079$
Simple FT pRotatE-ComplEx	$0.954\pm 0.009$	$\underset{\_}{0.5\pm 0.058}$	$\underset{\_}{0.454\pm 0.052}$	$\underset{\_}{0.569\pm 0.024}$	$0.854\pm 0.073$
Simple FT HolE-ComplEx	$\underset{\_}{0.965\pm 0.007}$	$0.418\pm 0.058$	$0.383\pm 0.053$	$\mathbf{0.571}\pm \mathbf{0.042}$	$0.9\pm 0.046$
Simple FT pRotatE-RotatE	$0.806\pm 0.039$	$0.229\pm 0.027$	$0.035\pm 0.016$	$0.131\pm 0.032$	$0.782\pm 0.157$
Simple FT HAKE-HAKE	$0.893\pm 0.046$	$0.104\pm 0.051$	$-0.003\pm 0.031$	$0.037\pm 0.033$	$0.588\pm 0.332$
Complex one-hot	$0.656\pm 0.069$	$0.422\pm 0.075$	$0.078\pm 0.053$	$0.124\pm 0.036$	$0.645\pm 0.178$
Complex PT HAKE-DistMult	$0.923\pm 0.013$	$0.434\pm 0.059$	$0.357\pm 0.052$	$0.488\pm 0.074$	$0.808\pm 0.07$
Complex PT HolE-DistMult	$0.949\pm 0.016$	$0.38\pm 0.084$	$0.33\pm 0.076$	$0.443\pm 0.089$	$0.805\pm 0.07$
Complex PT ConvKB-DistMult	$0.942\pm 0.01$	$0.387\pm 0.038$	$0.329\pm 0.039$	$0.484\pm 0.066$	$0.817\pm 0.052$
Complex PT HolE-RotatE	$0.932\pm 0.014$	$0.15\pm 0.018$	$0.082\pm 0.023$	$0.168\pm 0.015$	$0.861\pm 0.064$
Complex PT TransE-HAKE	$0.756\pm 0.047$	$0.19\pm 0.077$	$-0.054\pm 0.089$	$0.057\pm 0.046$	$0.742\pm 0.253$
Complex FT HAKE-DistMult	$0.925\pm 0.021$	$\underset{\_}{0.513\pm 0.064}$	$\underset{\_}{0.437\pm 0.058}$	$0.522\pm 0.034$	$0.83\pm 0.09$
Complex FT HolE-DistMult	$0.926\pm 0.015$	$\mathbf{0.536}\pm \mathbf{0.03}$	$\mathbf{0.462}\pm \mathbf{0.03}$	$\underset{\_}{0.543\pm 0.039}$	$0.81\pm 0.084$
Complex FT ConvKB-DistMult	$0.933\pm 0.01$	$\underset{\_}{0.525\pm 0.065}$	$\underset{\_}{0.459\pm 0.063}$	$\underset{\_}{0.55\pm 0.04}$	$0.746\pm 0.122$
Complex FT HolE-RotatE	$0.863\pm 0.057$	$0.194\pm 0.053$	$0.057\pm 0.015$	$0.11\pm 0.021$	$0.81\pm 0.278$
Complex FT TransE-HAKE	$0.892\pm 0.027$	$0.075\pm 0.043$	$-0.033\pm 0.049$	$0.072\pm 0.048$	$0.958\pm 0.077$

Tables 11–14 show the results for each of the data sampling strategies (i)–(iv), respectively. The tables include the three best models (based on $\mathrm{YI}$ ) for the baseline model using one-hot and pre-trained (PT) KG embeddings, and the fine-tuning (FT) models using the same combination of KGE models as the selected PT-based models. We have also included a model with middling performance (i.e., 40 out of 81 models) and the worst performing model. Note that for the PT- and FT-based models we have evaluated 81 combinations ${\mathrm{KGE}}_{\mathit{C}}$ - ${\mathrm{KGE}}_{\mathit{S}}$ of KGE models. All models were evaluated using the simple and complex MLP settings. For example, the model Complex FT DistMult-HolE denotes that fine-tuning was used together with the complex MLP setting, and DistMult was selected to embed the chemicals ${\mathrm{KG}}_{\mathit{C}}$ while HolE was used to embed the species ${\mathrm{KG}}_{\mathit{S}}$ . We present the mean and standard deviation over 10 evaluation runs, i.e., we re-initialize and re-train the models 10 times. Results highlighted in bold are the best mean results of the corresponding metrics. Underlined results are where there is a $⩾32\mathrm{\%}$ chance that a single run outperforms the best mean (i.e., one standard deviation contains $68\mathrm{\%}$ of results, assuming normally distribute results).37

³⁷

Note that we only consider the best mean result and not the standard deviation in both directions.

Overall, models with the complex setting and fine-tuning are needed as the data sampling strategies become more challenging. Moreover, all models favour sensitivity over specificity at default decision threshold (0.5). This is down to the imbalance in the data. We can see the imbalance by ${\mathit{\tau }}_{max}$ , it is $\mathrm{>}0.5$ for most models. As we use a log-loss instead of a discrete loss, this is to be expected for imbalanced data.

For settings (iii) and (iv) the performance drops and the standard deviation increases compared to the other strategies. This large standard deviation leads to large overlaps in quantiles among top-3 models in all categories, such that, by chance, one of these models could perform best in one individual evaluation.

7.2.1. One-hot baseline models

For the sampling strategy (i) the one-hot baseline models perform well, especially, with the complex one-hot model. This complex model is equivalent in terms of $\mathrm{YI}$ as the best simple pre-trained model. The story is largely the same in setting (ii), where the complex one-hot model performs within $1.5\mathrm{\%}$ of the best simple pre-trained models. With strategies (iii) and (iv) the one-hot models degrade, especially in strategy (iv) where the Youden’s index is near zero ( $\mathrm{<}0.1$ ). This is expected as the one-hot baseline models lack important background information about the entities, specially for unseen chemicals and species, that the KG embedding models aim at capturing.

7.2.2. Baseline with pre-trained KG embeddings

We can see that the PT-based models do not lead to an important improvement with respect to ${\mathrm{YI}}_{max}$ in sampling strategy (i). The top-1 complex PT model, however, yields a better balance between sensitivity and specificity leading to an improved $\mathrm{YI}$ over the complex one-hot models. The two middling performing models, Simple PT pRotatE-ConvE and Complex PT ComplEx-DistMult, still retain a decent level of performance.

The results with the strategy (ii) are similar to strategy (i), the delta in $\mathrm{YI}$ between the simple and the complex PT-based models are about $5\mathrm{\%}$ . This slight improvement is due to the increased balance between sensitivity and specificity which in turn leads to a higher $\mathrm{YI}$ .

In the sampling strategy (iii) we can observe that the improvement of the PT-based models over the one-hot models increases. The increase is up to $25\mathrm{\%}$ in $\mathrm{YI}$ of the best PT-based model over the best one-hot model. In addition, we observe in this strategy that the standard deviation increases, especially in specificity, leading to a large portion of the models that are within one standard deviation of the best model in terms of  $\mathrm{YI}$ .

Finally, the impact of using a PT-based models is strengthen in strategy (iv). The delta between the one-hot and PT-based models is up to $40\mathrm{\%}$ in $\mathrm{YI}$ , and larger for ${\mathrm{YI}}_{max}$ . We see that all models struggle with specificity in this setting, this is down to the difficulty of predicting true negatives. This also leads to a larger variation, with certain models yielding standard deviation in the same order of magnitude as the metric (e.g., Simple FT HAKE-ComplEx).

7.2.3. Fine-tuning optimization model

The FT-based models, with some exceptions, improve the results over the PT-based models, most notably in sampling strategies (iii) and (iv). For example, the FT-based models Complex FT HolE-DistMult and Simple FT HolE-ComplEx are the best models in terms of $\mathrm{YI}$ and ${\mathrm{YI}}_{max}$ in strategy (iv), respectively. We can also see in strategies (i) and (ii) that the FT-based models improve middling and worst performing PT-based models, e.g., Simple FT RotatE-ConvE in strategy (i) improves from $\mathrm{YI}=0.021$ to $\mathrm{YI}=0.726$ using fine-tuning of the KG embeddings. The results are expected as the fine-tuned KG embeddings are tailored to the effect prediction task.

7.3. KG embedding analysis

In this section we look at correlations between KGE model choices and prediction performance. KGE models are designed to capture certain structures in the data, and this can give some explanation of which parts of the KGs are important for prediction.

First, in Table 15 we show how many times a KGE model is used when regarding the top 10 performing combinations (out of the total 81 possible combinations). We focus on the choices when using the simple MLP setting to reduce the influence of the non-linear transforms on the embeddings.

Table 15

Usage of KGE models for each sampling strategy in simple MLP setting in top-10 performing combinations. Note that, there is one model for the ${\mathrm{KG}}_{\mathit{C}}$ and one for ${\mathrm{KG}}_{\mathit{S}}$ , such that there is a total of 20 models per sampling strategy. Notation: ‘used in ${\mathrm{KG}}_{\mathit{C}}$ / used in ${\mathrm{KG}}_{\mathit{S}}$ ’, e.g., HAKE, $2\mathrm{/}8$ in sampling strategy (i), indicates that HAKE is used to embed ${\mathrm{KG}}_{\mathit{C}}$ 2 out of top-10 combinations and it is used to embed ${\mathrm{KG}}_{\mathit{S}}$ 8 out of top-10 combinations

KGE model	# uses (i)	# uses (ii)	# uses (iii)	# uses (iv)
DistMult	$1\mathrm{/}0$	$0\mathrm{/}1$	$1\mathrm{/}7$	$0\mathrm{/}4$
ComplEx	$1\mathrm{/}1$	$1\mathrm{/}3$	$2\mathrm{/}1$	$1\mathrm{/}5$
HolE	$2\mathrm{/}0$	$1\mathrm{/}0$	$1\mathrm{/}0$	$1\mathrm{/}0$
Total decomposition	$4\mathrm{/}1$	$2\mathrm{/}4$	$4\mathrm{/}8$	$2\mathrm{/}9$
TransE	$1\mathrm{/}0$	$2\mathrm{/}0$	$1\mathrm{/}2$	$0\mathrm{/}0$
RotatE	$0\mathrm{/}0$	$0\mathrm{/}0$	$0\mathrm{/}0$	$1\mathrm{/}0$
pRotatE	$1\mathrm{/}0$	$1\mathrm{/}0$	$1\mathrm{/}0$	$3\mathrm{/}0$
HAKE	$2\mathrm{/}8$	$3\mathrm{/}5$	$1\mathrm{/}0$	$2\mathrm{/}0$
Total geometric	$4\mathrm{/}8$	$6\mathrm{/}5$	$3\mathrm{/}2$	$5\mathrm{/}0$
ConvKB	$1\mathrm{/}1$	$0\mathrm{/}1$	$2\mathrm{/}0$	$0\mathrm{/}1$
ConvE	$1\mathrm{/}0$	$2\mathrm{/}0$	$1\mathrm{/}0$	$2\mathrm{/}0$
Total convolutional	$2\mathrm{/}1$	$2\mathrm{/}1$	$3\mathrm{/}0$	$2\mathrm{/}1$

Looking at Table 15 we can see that the KGE models used to embed the chemicals ${\mathrm{KG}}_{\mathit{C}}$ in the best performing models is distributed evenly across most models and settings. This indicates that the performance of the prediction models is not highly correlated with the use of a KGE model on ${\mathrm{KG}}_{\mathit{C}}$ . Referencing Table 7, the high relational density in ${\mathrm{KG}}_{\mathit{C}}$ can contribute to worse performance [62] and therefore equal distribution of models in Table 15. This is different for ${\mathrm{KG}}_{\mathit{S}}$ . For sampling strategies (i) and (ii), HAKE is extensively used in the top models to embed ${\mathrm{KG}}_{\mathit{S}}$ . HAKE is designed to embed hierarchies. Therefore, this indicates that in strategies (i) and (ii) the hierarchical structure of ${\mathrm{KG}}_{\mathit{S}}$ dwarfs the rest of the KG. ${\mathrm{KG}}_{\mathit{S}}$ has a higher entity density and lower entity entropy (Table 7) than ${\mathrm{KG}}_{\mathit{C}}$ . This should lead to higher performance generally, but might also lead to larger discrepancies between models as seen in Table 15.

The use of the decomposition models increase in strategies (iii) and (iv) for the embedding of ${\mathrm{KG}}_{\mathit{S}}$ , which indicates that KG structures, other than the hierarchy, are important. Overall, DistMult and ComplEx can be used to great effect in strategies (iii) and (iv) while the geometric model, HAKE, is more successful in the less challenging strategies (i) and (ii).

7.3.1. Explained variance

Explained variance is a measure of how many principal components are required to describe all components.38

³⁸

We use the scikit-learn implementation [60] based on [72].

In Fig. 6, we present how the $\mathrm{YI}$ metric depends on the explained variance of the top-10 principal components (i.e., ${\sum }_{\mathit{i}=1}^{10}{\mathit{pca}}_{\mathit{i}}$ ). We show all (81 per sampling strategy) PT-based prediction model results, simple MLP setting in Fig. 6a and complex setting in Fig. 6b. For example, in Fig. 6a, the best model in the strategy (iv), Simple PT pRotatE-ComplEx have a explained variance of 0.49 compared to the worst model, Simple PT HAKE-HAKE, with explained variance of 0.34. Coincidentally, these two points does not follow the trend lines in these figures which indicate negative correlation between $\mathrm{YI}$ and explained variance. The trend lines can be interpreted in two ways. First, it is counter-intuitive as we would expect more descriptive embeddings, i.e., larger explained variance, to perform better. On the other hand, the top-10 principal components may not be representative enough to capture the semantics of the KG embeddings, and thus, a large explained variance does not necessarily correlate with a high performance.

Fig. 6.

Relation between explained variance using 10 principal components and model performance represented as $\mathrm{YI}$ .

Fig. 7.

Relation between explained variance using 10 principal components and model performance represented as sensitivity.

Figure 7 represents the explained variance against sensitivity. We can see that the trend is flat for strategy (iv), but positive for strategies (i)-(iii). This means that the trends in Fig. 6 are explained by specificity rather than sensitivity. By balancing sensitivity and specificity, i.e., ${\mathrm{YI}}_{max}$ as seen in Fig. 8, the rate of change is reduced compared to $\mathrm{YI}$ in Fig. 6.

Fig. 8.

Relation between explained variance using 10 principal components and model performance represented as ${\mathrm{YI}}_{max}$ .

7.4. Example predictions

Table 16 shows a few examples of correct (TP and TN) and incorrect predictions (FN and FP).

Table 16

Example predictions by complex FT HolE-DistMult (best model) for sampling strategy (iv)

Chemical	Species	$log\mathrm{(}\mathit{\kappa }\mathrm{)}$	Predicted	Lethal	Classification
D001556 (hexachlorocyclohexane)	59899 (walking catfish)	−3.4	0.97	1 (yes)	TP
C037925 (benthiocarb)	7965 (sea urchins)	0.9	0.2	0 (no)	TN
D026023 (permethrin)	378420 (bivalves)	0.7	0.96	1 (yes)	TP
D011189 (potassium chloride)	938113 (megacyclops viridis)	6.7	0.27	1 (yes)	FN
C427526 (carfentrazone-ethyl)	208866 (eudicots)	−0.9	0.82	0 (no)	FP
D010278 (parathion)	201691 (green sunfish)	−0.9	0.86	0 (no)	FP

Benthiocarb and permethrin are both biocides with different targets: benthiocarb is a herbicide and permethrin is an insecticide. It is therefore not surprising that benthiocarb has a low predicted effect on sea urchins, while permethrin has a severe effect on bivalves.

There are several possible explanations for the failed predictions. A wrong prediction of potassium chloride toxicity to a marine copepod (Megacyclops viridis) could be due to the prediction model not being accurate enough for metal salts, or the copepod species being particularly sensitive to changes in osmolarity due to salt content. The wrong prediction of lack of herbicide toxicity (i.e., carfentrazone-ethyl) to a flower (i.e., eudicots) could be due to the fact that flowers, and plants in general, are severely underrepresented in the available effect prediction data.

8. Discussion

We have introduced the Toxicological Effect and Risk Assessment (TERA) knowledge graph and shown how we can directly use it in chemical effect prediction. The use of TERA improves the PT-based prediction models over the one-hot baselines. In the most challenging data sampling strategies, we have also seen the benefits of creating tailored (i.e., fine-tuned) KG embeddings in the FT-based prediction models.

8.1. TERA knowledge graph

The constructed knowledge graph consists of several sources from the ecotoxicological domain. There are three major parts in TERA: the effects data, the chemical data, and the species taxonomic data. Integrating each part has different challenges. The chemical and pharmacological communities have come a long way in annotating their data as knowledge graphs and ontologies. Here, selecting the correct subsets to work with the chemical effect prediction data was a major challenge. This had to be done based on mappings between effect data and chemical data that were extracted from Wikidata. We selected a relatively small subset of the chemical sub-KG to facilitate faster model training, however, still larger than the extracted fragment from the species sub-KG. The species sub-KG was created from tabular data and cleaned by removing several annotation labels with redundant information. This sub-KG was aligned using ontology alignment systems to the species taxonomy in the effects sub-KG. This required pre-processing of the KG, where it was divided into smaller parts such that the selected systems could perform the alignment. We used several standard ontologies to facilitate the transformation of the effect data into a knowledge graph. This involved not only automatic processes, but also an important amount of manual work.

Integrating more data into TERA involves the creation of mappings to the existing data. This is possible for a large amount of chemical datasets as Wikidata links multiple datasets, e.g., the chemical compound diethyltoluamide (wd:Q408389) has $\sim 35$ distinct identifiers. Biological data, both taxonomic and effects, might be harder to align to TERA as these mappings are not available in Wikidata. Here, ontology alignment systems play an important role to fill this gap.

The additional integrated data will give larger coverage of the domain, and thereby, improve model performance. However, adding more data will also increase the memory and time requirements of KGE models. This was bypassed in this work by reducing TERA to only relevant parts.

Adding additional domain knowledge is also critical in other applications, such as using TERA for data access.

8.2. Performance of prediction models

We have shown that the ability to embed some structure types of different KGE models largely impact the prediction models. We see that some KGE models fail to capture the semantics of the chemicals and the species, which leads to similar performance to the one-hot baselines. Moreover, in a few isolated cases the performance is reduced further which leads us to believe that the embeddings collapse in one or some dimensions, making it impossible to distinguish among entities.

We suspect that the even distribution of KGE models to embed ${\mathrm{KG}}_{\mathit{C}}$ (Table 15) in most settings is likely down to the structure of ${\mathrm{KG}}_{\mathit{C}}$ . This sub-KG has, unlike ${\mathrm{KG}}_{\mathit{S}}$ ’s tree structure, a forest structure, and models that can deal with trees (as in ${\mathrm{KG}}_{\mathit{S}}$ ) fail here, e.g., an entity in ${\mathrm{KG}}_{\mathit{C}}$ can have multiple parents, but only one grand-parent. In this case, some models may create very similar or the same embeddings for the parent nodes.

9. Conclusions and future work

TERA is a novel knowledge graph which includes large amounts of data required by ecological risk assessment. We have conducted an extensive evaluation of KGE embedding models in a novel and very challenging application domain. Moreover, we have shown the value of using TERA in an ecotoxicological effect prediction task. The fine-tuning optimization model architecture to adapt the KG embeddings to the prediction task has, to our knowledge, not been applied elsewhere.

9.1. Value for the ecotoxicology community

The creation of TERA is of great importance to future effect modelling and computational risk assessment approaches within ecotoxicology. Where the strategic goal is designing and developing prediction models to assess the hazard and risks of chemicals and their mixtures where traditional laboratory data cannot easily be acquired.

A great effort in the hazard and risk assessment of chemicals is the reduction of regulatory-mandated animal testing. Wide-scale predictive approaches, as described here, answer a direct and current need for generalized prediction frameworks. These can aid in identifying especially sensitive species and toxic chemicals. At the Norwegian Institute for Water Research (NIVA), TERA will be used in this regard and will support several research projects.

In environmental risk assessment it is often unfeasible to assess the hazard and risk a chemical poses to a local species in the environment. These species may not be suitable for lab testing, or may even be endangered and thus are protected by national or international legislation. The currently presented work provides an in silico approach to predict the hazard to such species based on the taxonomic position of the species within the tree of life.

From an economic perspective, TERA and the prediction models are useful tools to evaluate new industrial chemicals during the synthetic in silico stage. Candidate chemicals can be evaluated for their potential environmental hazard, which is in line with the Green Chemistry initiatives by authorities such as the European Parliament or the US Environmental Protection Agency.

The effect prediction using TERA is also in line with a larger shift in ecological risk assessment towards the use of artificial intelligence [80]. We also believe the development of TERA contributes to a methodological change in the community, and encourages others to make their data interoperable.

9.2. TERA as background knowledge

As mentioned, in this work we use TERA directly in prediction models. However, TERA could be used as background knowledge to improve many emerging techniques for toxicity prediction (e.g., [65]). These methods often use chemical features, images, fingerprints and so on as input, and machine learning methods such as Convolutional Neural Networks and Random Forests as prediction models [81,84]. These models are often uninterpretable, and the predictions lack domain explanations. TERA can also provide context for machine learning tasks such as pre-processing, feature extraction, transfer and zero/few-shot learning. Furthermore, the knowledge graph is a possible source for the (semantic) explanation of the predictions (e.g., [43]).

9.3. Benchmarking KG embedding models

We have shown that embedding TERA brings new challenges to state-of-the-art KGE models with respect to capturing the semantics of the chemicals and the species. Furthermore, as shown in Section 5.4 the sparsity-related measures indicate that TERA represent an interesting KG. KGE models could be benchmarked in a standard KG completion task or in a specific task such as the chemical effect prediction.

9.4. Value to the ontology alignment community

As mentioned in Section 5.2, there does not exist a complete and public alignment between ECOTOX species and the NCBI Taxonomy. Therefore the computed mappings can also be seen as a very relevant resource to the ecotoxicology community. The used alignment techniques achieve high scores for recall over the available (incomplete) reference mappings. However, aligning such large and challenging datasets requires preprocessing before ontology alignment systems can cope with them. We removed all nodes which did not share a word (or shared only a stop word) in labels across the two taxonomies. This quartered the size of ECOTOX and reduced NCBI Taxonomy 50 fold. However, the possible alignment between entities without labels is lost when reducing the dataset size. Thus, the alignment of ECOTOX and NCBI Taxonomy has the potential of becoming a new track of the Ontology Alignment Evaluation Initiative (OAEI) [52] to push the limits of large scale ontology alignment tools. Furthermore, the output of the different OAEI participants could be merged into a rich consensus alignment (e.g., as done in the phenotype-disease domain [28]) that could become the reference alignment to integrate ECOTOX and NCBI Taxonomy.

9.5. Future work

We plan to extend TERA to include a larger part of ChEBI (which ChEMBL is a part of). ChEBI includes relevant data on the interaction between chemicals and species at a cellular level, which may be very important for chemical effect prediction. In this work we only consider effect data from ECOTOX as this is the largest data set available, however, the inclusion of e.g., TOXCAST [75] is in our interest. New sources will always bring more coverage of the domain and will improve TERA for prediction, as background knowledge, and for data access.

We plan to evaluate the effect prediction under different parts of TERA, i.e., which sources in TERA provide value and which do not contribute in terms of the effect prediction. A similar effort in exploring different KG crawling techniques has been explored in [67]. In a similar vain, we plan to evaluate how materialization, via OWL reasoning, of TERA’s implicit triples affects prediction performance.

Finally, as mentioned already, some KGE models cannot deal with parts of the structure of TERA. An in-depth analysis of this is an interesting direction for future research. This could be solved by embedding the hierarchy separately, e.g., [50], or imposing restrictions on the embeddings, such as a minimum distance constraint.

9.6. Resources

We encourage feedback from domain researchers on extensions to TERA and associated tools.

A snapshot of TERA is available at

https://doi.org/10.5281/zenodo.3559865

This snapshot does not include data that is impractical to re-share (i.e., partial ${\mathrm{KG}}_{\mathit{C}}$ as described in Section 5). However, we include the full ${\mathrm{KG}}_{\mathit{E}}$ and ${\mathrm{KG}}_{\mathit{S}}$ .

All the material related to this project is available at

https://github.com/NIVA-Knowledge-Graph/

Source codes to create TERA are available in the TERA GitHub repository. The prediction models and data used for prediction can be found in the KGs_and_Effect_Prediction_2020 GitHub repository. The prediction models require the implementation of the KGE models from the KGE-Keras GitHub repository.