In silico toxicology: computational methods for the prediction of chemical toxicity.

Determining the toxicity of chemicals is necessary to identify their harmful effects on humans, animals, plants, or the environment. It is also one of the main steps in drug design. Animal models have been used for a long time for toxicity testing. However, in vivo animal tests are constrained by time, ethical considerations, and financial burden. Therefore, computational methods for estimating the toxicity of chemicals are considered useful. In silico toxicology is one type of toxicity assessment that uses computational methods to analyze, simulate, visualize, or predict the toxicity of chemicals. In silico toxicology aims to complement existing toxicity tests to predict toxicity, prioritize chemicals, guide toxicity tests, and minimize late-stage failures in drugs design. There are various methods for generating models to predict toxicity endpoints. We provide a comprehensive overview, explain, and compare the strengths and weaknesses of the existing modeling methods and algorithms for toxicity prediction with a particular (but not exclusive) emphasis on computational tools that can implement these methods and refer to expert systems that deploy the prediction models. Finally, we briefly review a number of new research directions in in silico toxicology and provide recommendations for designing in silico models. WIREs Comput Mol Sci 2016, 6:147-172. doi: 10.1002/wcms.1240 For further resources related to this article, please visit the WIREs website.

INTRODUCTION

Toxicity is a measure of any undesirable or adverse effect of chemicals. Specific types of these adverse effects are called toxicity endpoints, such as carcinogenicity or genotoxicity, and can be quantitative (e.g., LD50: lethal dose to 50% of tested individuals)1 or qualitative, such as binary (e.g., toxic or non‐toxic) or ordinary (e.g., low, moderate, or high toxicity).2 Toxicity tests aim to identify harmful effects caused by substances on humans, animals, plants, or the environment through acute‐exposure (single dose) or multiple‐exposure (multiple doses).3 Several factors determine the toxicity of chemicals, such as route of exposure (e.g., oral, dermal, inhalation), dose (amount of the chemical), frequency of exposure (e.g., single versus multiple exposure), duration of exposure (e.g., 96 h), ADME properties (absorption, distribution, metabolism, and excretion/elimination), biological properties (e.g., age, gender), and chemical properties.4

Animal models have been used for a long time for toxicity testing.3 However, in vitro toxicity tests became plausible due to the advances in high throughput screening.3 In silico toxicology (computational toxicology) is one type of toxicity assessment that uses computational resources (i.e., methods, algorithms, software, data, etc.) to organize, analyze, model, simulate, visualize, or predict toxicity of chemicals.5, 6 It is intertwined with in silico pharmacology, which uses information from computational tools to analyze beneficial or adverse effects of drugs for therapeutic purposes.5, 6

Computational methods aim to complement in vitro and in vivo toxicity tests to potentially minimize the need for animal testing, reduce the cost and time of toxicity tests, and improve toxicity prediction and safety assessment. In addition, computational methods have a unique advantage of being able to estimate chemicals for toxicity even before they are synthesized.7 In silico toxicology encompasses a wide variety of computational tools (Figure 1): (A) databases for storing data about chemicals, their toxicity, and chemical properties; (B) software for generating molecular descriptors; (C) simulation tools for systems biology and molecular dynamics; (D) modeling methods for toxicity prediction; (E) modeling tools such as statistical packages and software for generating prediction models; (F) expert systems that include pre‐built models in web servers or standalone applications for predicting toxicity; and (G) visualization tools.

Figure 1

In silico toxicology tools, steps to generate prediction models, and categories of prediction models.

The purpose of this study is to provide a comprehensive overview of existing modeling methods and algorithms for toxicity prediction (element D above), with a particular (but not exclusive) emphasis on computational tools that can implement these methods (element E), and expert systems that deploy the prediction models (element F). Due to the nature of this expanding field, this study cannot provide an exhaustive overview of all the seven in silico components mentioned above. Therefore, the reader is encouraged to refer to existing literature to get more information about toxicity databases,6, 8, 9, 10, 11 molecular descriptors generation software,12 toxicology simulation tools,13, 14 statistical modeling packages,12 expert systems6, 9, 11, 12, 15, 16, 17, and visualization tools.18

Generally, modeling methods include five major steps while developing prediction models19 (Figure 1): (1) gathering biological data that contain associations between chemicals and toxicity endpoints, (2) calculating molecular descriptors of the chemicals, (3) generating a prediction model, (4) evaluating the accuracy of the model, and (5) interpreting the model.

The scope of this review covers the third step, generating prediction models. We focus on using computational methods to predict toxicity of different types of substances such as drugs, other chemicals, mixtures, and nanomaterials both quantitatively and qualitatively. There are various methods to solve such problems, and each method has its strengths, limitations, scope of application, and specificity of interpretation. The goal is to find the most effective method for the problem at hand. However, all the five steps mentioned above are inter‐related. Therefore, we discuss the remaining steps whenever necessary. We divided this review into four main sections. First, we provide an explanation and discussion of the existing in silico methods. Then, we briefly discuss two special types of chemicals: mixtures and nanomaterials. Subsequently, we offer recommendations on how to develop and apply toxicity prediction models. Finally, we provide an overview of 21st century toxicology.

IN SILICO MODELING METHODS

Many in silico methods have been developed to predict the toxicity of chemicals. The methods we discuss here are chosen either because they illustrate the historical development of in silico toxicology or they represent the state‐of‐the‐art method for predicting toxicity. For each method, we provide (if applicable) a mathematical description, discussion of strengths and limitations, recommendations about when and why to use the method, and existing tools that implement the method. Additionally, for the sake of clarity, we keep equations and visual representations of models as general as possible.

Structural Alerts and Rule‐based Models

Structural alerts (SAs)12, 20 (also called toxicophores/toxic fragments17) are chemical structures that indicate or associate to toxicity.6, 12 SAs can consist of only one atom or several connected atoms.21 A combination of SAs may contribute to toxicity more than a single SA.21 SAs are often used in rules defined in the form ‘if A is B then T,’ where A is an SA, B is the value of the SA, and T is the toxicity prediction with assigned certainty level,6 as illustrated in the following example:

IF (chemical_substructure) IS (present) THEN (skin_sensitizer IS certain)

There are two main types of rule‐based models that we will consider: human‐based rules (HBRs) and induction‐based rules (IBRs).12 HBRs are derived from human knowledge of field experts or from literature, but IBRs are derived computationally.6, 12 HBRs are more accurate but are limited to human knowledge that could be incomplete or biased.6, 12 Moreover, updating HBRs is often impractical as it requires detailed literature analysis.21 On the contrary, IBRs can be generated efficiently from large datasets. IBRs may propose hypotheses about associations between chemical structural properties (or their combinations) and toxicity endpoints, which may not be identified through human insights.6, 21 IBRs are implemented using probabilities to determine if SAs correspond to the toxic or non‐toxic class. It is possible to have hybrid‐based rules systems that contain IBRs and HBRs, with new rules being generated computationally.12

It is easy to interpret and implement SAs.15 They are useful in drug design to determine how drugs should be altered to reduce their toxicity. Using structure to predict toxicity allows identifying the structure of potential metabolites.10 However, SAs have a number of limitations. SAs use only binary features (e.g., chemical structures are either present or absent) and only qualitative endpoints (e.g., carcinogenic or non‐carcinogenic).12 SAs do not provide insights into the biological pathways of toxicity and may not be sufficient for predicting toxicity. Depending on the concurrent absence or presence of other chemical properties, toxicity may decrease or increase.15 The list of SAs and rules may be incomplete, which may cause a large number of false negatives (i.e., toxic chemicals predicted as non‐toxic) in predictions.12, 15, 20

The last point is particularly important. It is necessary to understand how to interpret the output of SA models. If a chemical does not include SAs or does not match any toxicity rules, this does not indicate non‐toxicity.6 This is especially true for HBRs that usually include SAs or rules that indicate toxicity but do not include SAs or rules that indicate non‐toxicity.6 Therefore, in developing such models, it is necessary to ensure that the list of SAs and rules are comprehensive and that they are refined when more experimental data becomes available. However, there should be a balance between the list of SAs and rules, their comprehensiveness, and predictive power. If SAs and rules are diverse, they can be applied to a large number of chemicals, but this may increase false positives (i.e., non‐toxic chemicals predicted as toxic). However, if they are too narrow, they can be applied only to a small group of chemicals, and this may increase false negatives (i.e., toxic chemicals predicted as non‐toxic).

An example of SA list for skin sensitization was published in 1982 by Dupuis and Benezra.22 Another SA list was proposed by Ashby and Tennant19, 21, 23 in 1988 to predict carcinogenicity and mutagenicity. One of the most developed lists of carcinogenic SAs was proposed by Benigni and Bossa19, 21, 24 in the Organisation for Economic Co‐operation and Development (OECD) Quantitative Structure‐Activity Relationship (OSAR) Toolbox25 and in Toxtree.19, 21, 26 Recently, Benigni and Bossa published a new list of non‐genotoxic carcinogenic SAs.27 Other lists of SAs and rule‐based models were developed for hepatotoxicity,28 cytotoxicity,29 irritation/corrosion of skin,30 and eye31 and skin sensitization.19, 32

Several systems (also called ‘expert systems’20) provide pre‐built rule‐based and SAs lists, for example, Oncologic Cancer Expert System (OCES),33 Toxtree,15, 26 Derek Nexus,15, 34 HazardExpert15, 35 and Meteor.36 Other tools that can extract SAs from datasets that contains toxic or non‐toxic chemicals are reviewed in21, such as computer assisted structure elucidation (CASE),37 prediction of activity spectra for substances (PASS)38, and categorical‐structure activity relationship (cat‐SAR).39

Additionally, there are several approaches for extracting the longest frequent molecular substructures such as Apriori (based on breadth‐first search) and pattern growth (based on depth‐first search).21 Examples of algorithms that implement the pattern growth approach are reviewed in21, such as molecular fragment miner (mofa),40 graph‐based substructure pattern mining (gSpan),41 fast frequent subgraph mining (FFSM),42 and Graph/Sequence/Tree Extraction (gaston).43 Significant substructures capable of discriminating between toxic and non‐toxic chemicals can be extracted using an emerging chemical pattern approach as explained in ref. 21.

Chemical Category, Read‐Across, and Trend Analysis

A chemical category44 is a group of chemicals whose properties and toxicity effects are similar or follow a similar pattern.11, 45 Chemicals in the category are also called source chemicals. The OECD Guidance On Grouping Of Chemicals lists several methods for grouping, such as chemical identity and composition, physicochemical and ADME properties, mechanism of action (MoA), and chemical/biological interactions.46 Structural similarity is described in the OECD guidelines as the starting point for grouping, but it is also criticized for lacking a ‘scientifically supportable basis’ for grouping, and it can be used if impurities or other constituents in the chemical composition would not change toxicity.46

Read‐across is a method of predicting unknown toxicity of a chemical using similar chemicals (called chemical analogs) with known toxicity from the same chemical category.9, 11, 12, 45, 47 Trend analysis is a method of predicting toxicity of a chemical by analyzing toxicity trends (increase, decrease, or constant) of tested chemicals.9 A hypothetical example of trend analysis shows that when carbon chain length (CCL) increases, acute aquatic toxicity increases (Figure 2).45

Figure 2

Different approaches of read‐across: analog versus category approaches, interpolation versus extrapolation, category boundary and outliers.

Here, we focus on the read‐across method. A summary of different parameters that must be considered when designing a read‐across model is depicted in Figure 3 and explained later. Note, however, that the points discussed are similarly applicable to trend analysis.

Figure 3

Different properties of read‐across models.

There are two ways to develop a read‐across method12, 45, 48: analog approach (AN) (called one‐to‐one), which uses one or few analogs, and a category approach (CA) (called many‐to‐one), which uses many analogs. AN may be sensitive to outliers because two analogs may have different toxicity profiles.12 Using many analogs for CA is useful to detect trends within a category and may increase confidence in the toxicity predictions.11, 45, 48 CA requires defining a category boundary to determine if a chemical belongs to the category45 and implementing a ‘combination of predictions’ method for analogs that have conflicting toxicity profiles. A combination of predictions can be done using (if applicable) minimum, maximum, mode, median, average, linear, quadratic, or other nonlinear combinations of the predictions.47

Read‐across can be qualitative if the toxicity endpoint is qualitative; otherwise, read‐across is quantitative.6, 9, 12 Also, interpolation using source chemicals surrounding the target chemical (see Figure 2) is better than extrapolation from one side.17 In Figure 2, interpolation is used with the chemical that has CCL of length 6, but extrapolation is used with a chemical that has CCL of length 12.

Identifying similar chemicals can be done in two steps: representing chemicals as feature vectors of chemical properties, and then calculating similarity of chemicals. The first step is implemented using either binary or holographic fingerprints. A binary fingerprint is a feature vector of binary bits representing presence (1) or absence (0) of a property (e.g. presence of a methyl group).44, 47 However, a holographic fingerprint uses frequency of properties (e.g. number of methyl groups). Continuous chemical properties (e.g., melting point) can be used as well. A hierarchy of categories and subcategories can be better than a single feature vector. At each level of the hierarchy, a property is applied for category formation. Subsequently, categories are divided using another property to generate subcategories and so on. The hierarchy can allow for investigating the significance of properties and can simplify model interpretation.47 An example of hierarchal categories is provided in ref. 47. Statistical similarity of two chemicals can be calculated using different types of distances, such as Hamming, Euclidean, Cosine, Mahalanobis, Tanimoto distance, or linear or nonlinear relationships of the features.45, 47

There are several advantages of read‐across. Read‐across is transparent,16 easy to interpret and implement.44 Read‐across can model quantitative and qualitative toxicity endpoints, and it allows for a wide range of types of descriptors and similarity measures to be used to express similarity between chemicals.47

However, there are also limitations. Statistical similarity measures do not provide biological insight of toxicity.47 Moreover, complex similarity measures may complicate model interpretation.47 In reality, read‐across uses small datasets compared to other approaches such as QSAR because there are usually only a few analogs for a given chemical.47 Additionally, accuracy depends on the number and choice of analogs, similarity metrics, strength in chemicals’ similarity, chemical properties, and category boundaries.47 These parameters are very subjective, mutually dependent, endpoint‐specific, and may require expert opinions.17, 44, 45, 46, 47 Moreover, this approach could be inapplicable or inaccurate if analogs have conflicting toxicity profiles11 or the number of analog chemicals is insufficient. In these cases, the QSAR approach can be used.11, 12, 44, 45

Read‐across was applied to predict carcinogenicity,49 hepatoxicity,28 aquatic toxicity,50 reproductive toxicity,51 skin sensitization,52 and environmental toxicity.53 Examples of tools implementing read‐across are The OECD QSAR Toolbox,25 Toxmatch,54 ToxTree,26 AMBIT,55 AmbitDiscovery,56 AIM,57 DSSTox,58 or ChemIDplus.59 A detailed explanation of some of these tools is available in refs. 6, 9, 11, 16, 17, 44, 45, 48.

Dose–Response and Time–Response Models

Dose–response (or time–response) models are relationships between doses (or time) and the incidence of a defined biological effect (e.g., toxicity or mortality).60 A dose is ‘the total quantity of a substance administered to, taken up, or absorbed by an organism, organ, or tissue and can be measured with in vitro or in vivo experiments.’60 Time can be the time to produce a response or the time for recovery.61 Exposure time can be continuous, intermittent, or random, and exposure can be acute, short‐term, sub‐chronic , and chronic exposure.60 Time–dose models describe the relationship between time and dose for a constant response.62 Figure 4 shows different types of dose/time–response models. These models that describe relationships between response versus dose or time can be generated by regression to fit the data.

Figure 4

Different types of relationships for dose–response models. Similar relationships can be generated for time–response models.

The first dose–response model relates concentration (C) and time (t) with response (K), which is Haber's law (law of toxicity)60, 63:

However, Haber's law does not hold in many situations,63 and it does not take detoxification into consideration. The law assumes that any combination of concentration and time that has the same C × t product should produce the same level of toxicity. However, in reality, this is not the case. Toxicity of some chemicals can be more dependent on concentration than time. Subsequently, Haber's law was generalized. Let C denote a threshold concentration, and n and m are constants. Several well‐known generalizations of Haber's law are shown below:

Ostwald: (C – C) = K that emphasizes concentration64;

Druckery: C × t = K that emphasizes time63; and

Miller et al.: (C – C) m = K that emphasizes both concentration and time.63

One of the frequently measured responses is mortality (the number of deceased individuals). The Bliss method (or Probit model) (Figure 5)61, 63, 65 transforms time–mortality and dose–mortality relationships into linear relationships. This transformation follows the next steps: (a) link mortality frequency (the number of deceased subjects) to dose or time; (b) convert frequency to percentages (percentage of deceased subjects); (c) transform percentages to probits (probability unites) and express dose or time in on logarithmic scale. Probits are inferred doses (or time) that correspond to a given mortality percentage. Bliss devised a special table called ‘probits table’ to calculate the probits as explained in ref. 65. This method takes into consideration the variation of an individual's susceptibility to toxic agents. For example, a certain dose (or time exposure) can cause the mortality of some individuals but not others.61, 65

Figure 5

Bliss method. (a) Plot mortality frequency (the number of dead subjects) versus dose or time. (b) Convert frequency to percentages (percentage of deceased subjects). (c) Transform percentages to probits and transform dose or time to logarithms.

There are many inherited differences between dose–mortality and time–mortality models. Time–mortality curves are based on the same individuals whose susceptibility is measured at specific intervals. The percentage of mortality at a given interval cannot be less than that of the preceding interval, and the susceptibility of individuals in successive time intervals are correlated. However, dose–mortality curves are based on different individuals for each dose. Therefore, susceptibility of individuals at successive doses is unrelated especially if there are individuals who have a high toxicity resistance.61

The effectiveness of time–mortality curves depends on the ‘whole’ distribution of susceptibilities and their relationship to the response. Time–mortality curves that measure the response time can be incomplete for small doses due to individuals who have a high resistance and fail to show the measured response. Similarly, time–response curves that measure the recovery time may be incomplete for large doses if some individuals fail to recover. Bliss explained how to estimate the truncated distribution of time–mortality models.61

Miller et al.63 proposed a three‐dimensional model for concentration–time–response that can reliably interpolate within the scope of experimental data, and they provided an estimation of error when extrapolating outside the scope. Recently, Brown and Foureman62 used a time–concentration–response model to generalize the concentration–response models using time as a parameter.

There are many advantages of time–response, dose–response, and dose–time–response models: ease of interpretation and implementation, consideration of dose and time of exposure, interpolation of effects between different doses of the same chemical within the range of experimental data61, 63 using dose–response models, and interpolation between different exposure times for the same toxicant and dose within the range of experimental data61, 63 using time–response models.

However, there are many limitations. The three models cannot extrapolate to other chemicals.60 Additionally, time–response models cannot extrapolate to other doses of the same chemical.61, 63 Time–response models require that tested individuals have uniform susceptibility levels,61 or these models may be unreliable if some individuals have an extremely low or high resistance. If time intervals are long, time–response models may overestimate or underestimate the response at a given moment.61 The three models do not take into consideration target tissue, biological process, ADME, toxicokinetics, toxicodynamics, detoxification, damage or repair, or chemical properties.60, 63

These time–response and dose–response models are complementary to one another and must be used together to achieve reliable conclusions. Several databases include dose–response data such as CEBS,6, 66 PubChem,18, 67 and ToxRefDB.8, 68 These models were used, for example, for modeling rectal cancer,69 mutagenicity,70 and developmental toxicity.71

Pharmacokinetic Models and Pharmacodynamic Models

Pharmacokinetic (PK) models relate chemical concentration in tissues to time, estimate the amount of chemicals in different parts of the body, and quantify ADME processes.13, 72 Toxicokinetic models are PK models used to relate chemical concentration in tissues to the time of toxic responses. PK models can be compartmental and non‐compartmental.60, 72 A compartment is the whole or part of an organism in which the concentration is uniform.73 Compartmental models consist of one or more compartments, and each compartment is usually represented by differential equations.60, 72, 74

One‐compartment models represent the whole body as a single compartment, assume rapid equilibrium of chemical concentration within the body after administration, and do not consider the time to distribute the chemical. The concentration C at a given time t is computed by72

where C is the initial concentration and k is the elimination constant. The plotting log of concentration versus time results in a straight line of slope (−k).72 However, these models do not consider the distribution time of chemicals. Additionally, concentrations in some organs reach equilibrium faster than in others. Two‐compartment models consist of two compartments: central (for rapidly‐perfused tissues e.g., liver or kidney) and peripheral (for slowly perfused tissues e.g., muscle or skin). Each compartment is represented by a differential equation similar to the one‐compartment models. After solving the coupled equations, the concentration is the sum of two exponential terms of time (interpreted as distribution phase with initial concentration C and slope ‐a and elimination phase with initial concentration C and slope ‐b). The concentration C based on this model is represented by72

These models, however, cannot extrapolate between species or provide a mechanistic insight.72 On the other hand, physiologically based pharmacokinetic (PBPK) models include, in addition to concentration and time, physiological descriptors of tissues and ADME processes such as volumes, blood flows, chemical binding/partitioning, metabolisms, or excretions.13, 60 PBPK models represent each organ as a compartment, represented by a differential equation that includes PK parameters.13, 60, 72, 74 An organ can be split into several compartments if there is a high variability in organ tissue. Also, one compartment can represent several similar organs.72 A general PBPK model to calculate plasma concentration (C) uses a feature vector of PK parameters (θ), time (t), and dose (X)74 as follows:

where f is a function that models the relationship. Because equation structure and the physiological parameters are tissue specific, PBPK models allow for interspecies extrapolation and provide a mechanistic basis of ADME.11, 60, 72 PBPK models can convert administered doses to tissue dosimetry,60 which is ‘the amount of chemical that is distributed to a tissue or part of a tissue,’13 and generate concentration versus time models.11

Pharmacodynamic (PD) models relate a biological response to the concentration of chemical in tissue.72 Toxicodynamic models are PD models that relate toxicity to the concentration of the chemical. PD models that are based on anatomy, physiology, biochemistry, and biology are called physiologically based pharmacodynamic (PBPD) models.75 Similar to dose–response models, PD models can be linear or nonlinear. Linear models should be used with caution because they do not consider the upper limit of responses and assume that responses always increase when concentrations increase.72 Similar to PBPK models, PBPD can be described by differential equations. A general PBPD model calculates the response (R) using a feature vector of PD parameters (θ), plasma concentration (C), which is calculated using the PBPK model given above, or biophase concentration (C), and chemical‐independent system parameters (Z)74 can be represented as:

where f is a function that models the relationship. PD models can be combined with PK models.72 The resulting model is called biologically based dose–response models (BBDR) and can be used to relate doses with responses.13, 60, 72, 75 In addition to PK and PD parameters, BBDR may include biological parameters such as cell division rates, mortality rates, or production rates of hormones.13, 60, 76 BBDR models are more powerful than dose–response models because the former consider time‐dependent changes of concentration and can extrapolate at low doses and between species.60, 75

There are many advantages for the models discussed in this section. Determining internal doses rather than administered doses and key metabolites allows for a more direct relationship with the response.60 Additionally, using ADME, PK, and PD properties permits route‐to‐route and species‐to‐species (e.g., animal‐to‐human) extrapolations and in vitro‐to‐in vivo extrapolation.60 BBDR is useful for extrapolating at low doses. Such low doses provide realistic estimates for human toxicity as human exposure to toxicants is at much lower doses than those tested on animals.75

However, there are a number of disadvantages. PK and PD parameters may be unavailable or inaccurate. In such cases, the parameters are estimated using in vitro‐to‐in vivo or species‐to‐species extrapolation.72 Otherwise, QSAR modeling could be more appropriate because it depends only on molecular descriptors.11, 60 Additionally, if biological data is not available, empirical dose–response models are used instead of BBPD. Using BBDR for extrapolation between species assumes that the relationship between dose and response in animals is the same in humans.76, 77 The same problem applies when using animal studies to estimate PK or PD parameters for modeling toxicity in humans.78 Although BBDR models have been proposed more than 20 years ago as a tool to minimize uncertainty for low‐dose and interspecies extrapolation, it was recently shown that BBDR has not progressed to reach such expectations due to uncertainty in modeled parameters and data, limited applicability of BBDR models to a small group of chemicals, or inherited complexity of BBDR models or toxicity mechanisms as discussed in ref. 76. Expert knowledge is required for defining MoA, toxicity pathways and chemical interactions that cause the response.

Different types of PK and PD models are reviewed in74 and summarized in Supplementary Table S1. An example of developing a PBPK model is available in ref. 79. Also, methods for estimating PK parameters are reviewed in ref. 80. Examples of PK and PD modeling tools are WinNonlin,81 Kinetica,82 and ADAPT 5.83 For example, PBPK was used for route‐to‐route extrapolation,84 toxicity and risk assessment,85 and carcinogenicity assessment.86

Uncertainty Factor Models

Uncertainty factors (UFs) (also called assessment/extrapolation/risk factors) are used for assessing risk from chemical exposure or the recommended daily intake of chemicals.87 A UF model is the simplest form of model for inter‐species extrapolation (e.g., from animals to humans), intra‐species extrapolation (e.g., from healthy people to special groups of the population such as elderly people, pregnant women, children, and fetuses), or exposure duration extrapolation (e.g., from short exposure to long exposure). It requires two main factors88: no observed adverse effect levels (NOAEL), which is the highest dose not exhibiting observable toxicity and a UF, which is a numerical value to account for variability in inter‐species, intra‐species, exposure duration, or exposed dose. Extrapolation is done by dividing NOAEL by UF.

However, there are two limitations for using NOAEL88: (1) the definition of NOAEL indicates the absence of the ‘appreciable risk’ of toxicity, but it does not indicate a zero‐effect threshold; and (2) NOAEL values are not constants and can vary depending on experimental designs such as the number of tested animals, number of doses, and toxicity endpoints. It was shown that low statistical power (e.g., a small number of tested animals or a small number of tested doses) would result in higher NOAEL. However, it is possible to use a least observable adverse effect level (LOAEL, which is the least dose or concentration that causes the observed effect) or to use a benchmark dose level (BMDL, which is ‘the lower statistical confidence limit of the dose resulting in a predetermined response’) if NOAEL is not available.87, 88

In addition to UFs, ‘modifying factors’ (MFs) are used to account for uncertainties in the data and the database. Additionally, ‘safety factors’ (SFs) are used for irreversible effects, such as teratogenicity and non‐genotoxic carcinogenicity. Although, existing UFs account for intra‐species variability, the use of additional factors for child safety is recommended. The values of MFs and SFs cannot exceed 10.87

UFs are necessary to estimate reference dose (RfD) and reference concentration (RfC). RfD or RfC ‘provide quantitative information for use in risk assessments for health effects known or assumed to be produced through a nonlinear (presumed threshold) mode of action.60’ The reference values are calculated as

where POD is the point of departure (e.g., NOAEL, LOAEL, or BMDL).60 A default UF of 100 was first proposed in 1954.87 However, this default value does not account for the quality of the database, the nature of the effect, the duration of the exposure, route‐to‐route extrapolation, and consideration of special groups of the population. Therefore, several factors have been calculated by different agencies as explained in87, 89 and summarized in Supplementary Table S2.

There are several advantages of UF models. It is easy to implement and understand them. They provide adequate safety levels for a single chemical and mixtures of chemicals.88 Additionally, they account for inter‐species and inter‐individual as well as PK and PD differences. However, there are some limitations of UF models. Default UFs or sub‐factors are not conservative nor do they assume the worst‐case scenario. Therefore, extrapolated safety levels of chemicals are not always below the realistic safety threshold for humans.88 These models cannot be used to extrapolate toxicity levels of genotoxic carcinogens as these chemicals always cause toxicity effects that are proportional to the dose, even at small doses.87

Quantitative Structure–Activity Relationship

Quantitative structure–activity relationship (QSAR) is a family of models that uses molecular descriptors to predict chemicals’ toxicity. It is assumed that chemicals that fit the same QSAR model may work through the same mechanism.10 A general QSAR model to predict toxicity (T) using a feature vector of chemical properties (θ) and a function f that calculates T given θ is

A local QSAR is generated from congeneric chemicals (i.e., similar chemicals); otherwise, the model is a global QSAR6 if it was made from diverse chemicals. Local QSARs are more accurate as they are customized for specific chemicals. However, there is an overhead to develop a local QSAR for each type of chemical. Therefore, global QSARs are more practical but may be less accurate. Local QSARs can also provide insight on the MoA of specific chemicals, which global QSARs may overlook.

Quantitative Structure Toxicity/Property Relationship (QSTR/QSPR) models are QSAR models that predict toxicity and chemical properties, respectively.10, 19 Structure activity relationships (SAR) are used for categorical endpoints.19 There are different types of models in the QSAR family as summarized in Supplementary Table S3.

There are two main steps to develop a QSAR model: generating molecular descriptors and then generating models to fit the data. Several types of molecular descriptors can be used to describe chemicals as summarized in Supplementary Table S3. Therefore, feature selection algorithms based on, for example, simulated annealing, genetic algorithm, or principal component analysis can be used.5, 19 If there are a small number of descriptors, using two‐dimensional scatter plots of each descriptor versus the biological activity can help identify significant descriptors19 (Figure 6).

Figure 6

2D scatter plots of molecular descriptors and toxicity levels. (a) no correlation between molecular descriptor 1 and the toxicity endpoint. (b) and (c) linear and nonlinear relationships between the molecular descriptors 2 and 3, respectively, with the toxicity endpoint. (b) and (c) can be modeled with linear and nonlinear algorithms, respectively.

There are several types of algorithms to generate QSAR models: linear models such as those based on linear regression analysis, multiple linear regression and partial least squares for continuous endpoints, and linear discriminant analysis for categorical endpoints5, 19; nonlinear models such as artificial neural networks or support vector machines5, 19; and data‐driven models such as those based on decision trees, clustering, Naïve Bayes, and K‐nearest neighbor.7

A comparison of different machine learning and regression models is provided in ref. 7. Linear models are simpler and, in general, require tuning fewer parameters than nonlinear models. However, many relationships between chemicals and toxicity are nonlinear. Therefore, nonlinear models are commonly used for developing QSARs. The two‐dimensional scatter plots can help identify the type of regression models19 as illustrated in Figure 6.

Additionally, SAR landscapes are three‐dimensional plots through which one can visualize structure–activity relationships. The X–Y plane represents the molecular descriptors, and the Z‐axis represents response. Figure 7 shows a hypothetical example of a SAR landscape. The smooth region corresponds to chemicals that have a similar structure and similar activity. However, the ragged region corresponds to chemicals that have a similar structure but different activity levels (also called activity cliffs). The activity cliffs are the most interesting part of the SAR landscape.18 They show that small structural changes correspond to huge changes in activity. Additionally, they affect the performance of machine learning models, either because these regions are discarded as outliers, cause over‐fitting, complicate the prediction models, or increase the prediction error while generating the model.

Figure 7

SAR landscapes.

SAR landscapes can be visualized using SAR maps. SAR maps are two‐dimensional plots of activity similarity versus structure similarity that characterize SAR landscapes through four regions18 as shown in Figure 8. Moreover, a structure activity landscape index (SALI) and a structure activity index (SARI) can be used to analyze SAR landscapes as explained in ref. 90.

Figure 8

Four regions in SAR maps are scaffold hops, smooth region, nondescript, and activity cliffs.

Historically, one of the early QSAR models was developed in 1962 by Hansch et al. 19 in which the log of chemical concentration (C) was estimated using the octanol/water partition coefficient (π) and the Hammett constant (σ)91:

If the coefficient of a descriptor is positive, there is a positive relationship between the toxicity endpoint and the descriptor; otherwise, there is a negative relationship.92

Examples of QSARs for predicting toxicity of aromatic nitro compounds, nitrobenzene compounds, cytotoxicity of TIBO derivatives, and carcinogenicity of sulfa drugs are explained in ref. 5. A discussion of the performance of QSAR models to predict carcinogenicity, mutagenicity, reproductive, and developmental toxicity endpoints are available in ref. 6. Case studies on applications of QSAR to skin sensitization and developmental toxicity are available in refs. 93 and 12, respectively.

There are many tools that provide pre‐built QSAR models such as OECD QSAR Toolbox,25 TopKat,94 Derek Nexus,34 HazardExpert,35 VEGA,95 and METEOR.36 Their characteristics are summarized together with those of other QSAR‐based tools in.10, 12, 17, 96 Case studies on combining the results of different prediction tools are available in refs. 15, 17. However, specialized software tools for generating QSAR models such as ADAPT and TOPKAT include databases for toxicity data and can calculate molecular descriptors.12 Additionally, several stand‐alone databases have been compiled to provide toxicity data and/or molecular descriptors as summarized in refs. 6, 7, 10.

There are several advantages of QSAR models. They are easy to interpret if the descriptors are meaningful. They can model categorical and continuous toxicity endpoints and molecular descriptors and toxic and non‐toxic chemicals. Using different types of descriptors allows for modeling complex endpoints.6, 11

However, QSARs may not be always applicable. QSARs require a large number of chemicals in model development to achieve statistical significance.5 Additionally, QSARs require using feature selection to identify the most significant and independent molecular descriptors, and a large number of descriptors makes the multidimensional space complex and fragmented.97 QSARs cannot be used for extrapolation between species, routes of exposure, or doses unless biological data is used. Moreover, QSARs may not be biologically interpretable, and QSARs do not take dose, duration, or metabolites into consideration.

A brief description of all the tools mentioned in the IN SILICO MODELING METHODS section is available in Table 1.

Table 1

Summary of In Silico Modeling Methods

Method	Definition	Approaches	Advantages	Limitations	Existing Software or Databases
Structural alerts (SAs) and rule‐based models	SAs are chemical structures that indicate or associate to toxicity.	Human‐based rules Induction‐based rules Apriori (based on breadth‐first search) Pattern growth (based on depth‐first search) such as mofa, gSpan, FFSM, and gaston.	1) It is easy to interpret and implement SAs. 2) SAs allow determining how chemicals should be altered to reduce their toxicity. 3) SAs allow identifying the structure of potential metabolites.	1) This method can indicate only the presence or absence of SAs. 2) SAs do not provide insight into the biological pathways of toxicity. 3) The list of SAs may be incomplete, which may increase false negatives.	OECD QSAR Toxtree OCES Derek Nexus HazardExpert Meteor CASE PASS cat‐SAR
Read‐across (RA)	A method of predicting unknown toxicity of a chemical using similar chemicals with known toxicity from the same chemical category	Analog approach (one‐to‐one) Category approach (many‐to‐one) Qualitative and quantitative RA Interpolation and extrapolation RA	1) RA is transparent, easy to interpret, and implement. 2) RA can model quantitative and qualitative toxicity endpoints, and uses many types of descriptors and similarity measures.	1) RA uses small datasets. 2) Accuracy depends on the number and choice of analogs, similarity metrics, strength in chemicals’ similarity, chemical properties, and category boundary. 3) RA may be inaccurate if analogs have conflicting toxicity profiles.	OECD QSAR Toxmatch ToxTree AMBIT AmbitDiscovery AIM DSSTox ChemIDplus
Dose–response (DR) and time–response (TR) models	Dose–response (or time–response) models are relationships between doses (or time) and the incidence of a defined biological effect.	Haber's law and its generalizations Bliss method (Probit model) 3D time‐dose–response models	1) Ease of interpretation and implementation 2) Consideration of dose and time of exposure 3) Interpolation of effects between different doses and exposure times	1) DR and TR models cannot extrapolate to other chemicals. 2) TR models require tested individuals to have uniform susceptibility levels. 3) DR and TR models do not consider target tissue, or chemical properties.	CEBS PubChem ToxRefDB
Pharmacokinetic (PK) and pharmacodynamic (PD) models	PK models calculate concentration at a given time. PD models calculate effect at a given concentration	One‐compartment models Two‐compartment models PBPK, PBPD and BBDR models	1) PK models determine internal doses rather than administered doses. 2) PK and PD models permit route‐to‐route, species‐to‐species, and in vitro‐to‐in vivo extrapolation. 3) BBDR is useful for extrapolating at low doses.	1) PK and PD parameters may be unavailable or inaccurate. 2) Extrapolation between species assumes that the relationship between dose and response in certain species is the same as in the other. 2) Expert knowledge is required for defining MoA, toxicity pathway and chemical interactions.	WinNonlin Kinetica ADAPT
Uncertainty factors (UFs) models	UF is a numerical value to account for variability in inter‐species, intra‐species, exposure duration, or exposed dose	Extrapolation using NOAEL, LOAEL, or BMDL. RfD and RfC models Modifying factors and safety factors	1) It easy to implement and understand UF models. 2) They provide adequate safety levels for single chemical and mixtures of chemicals. 3) They account for inter‐species and inter–individual, PK and PD differences.	Default UFs or sub‐factors are not conservative nor do they assume the worst‐case scenario.
Quantitative structure‐activity relationship (QSAR) models	QSAR is a family of models that use molecular descriptors to predict chemicals’ toxicity.	Local and global QSAR SAR, QSTR, and QSPR SAR landscapes and maps	1) QSARA models are easy to interpret if the descriptors are meaningful. 2) They can model categorical and continuous toxicity endpoints, and toxic and non‐toxic chemicals. 3) Using different types of descriptors allows for modeling complex endpoints.	1) QSARs require large datasets. 2) QSARs may require using features selection. 3) QSARs cannot be used for extrapolation between species, routes of exposure or doses, unless biological data are used. 4) QSARs do not take dose, or duration into consideration.	OECD QSAR TopKat Derek Nexus HazardExpert VEGA METEOR

SPECIAL CASES

This section focuses on special types of chemicals or toxicity endpoints that require new prediction or analysis methods.

Nanotoxicity

Nanotoxicity is the study of adverse effects caused by nanomaterials. Nanomaterials are small particles on the nanoscale (10−9 m) size range. When a particle size is decreased within the nanoscale size range, its physical and chemical properties are changed, affecting its toxicity. It was found that nanoparticles cause different or worse toxicity effects than the larger particles of the same substance.98, 99 A nanoparticle can be toxic even if the particle is not toxic at a larger size. The small size of nanomaterials facilitates cell membrane penetration and biodistribution.98, 99 Other properties that affect toxicity of nanomaterials are as follows100:

The shape of nanomaterials affects toxicological responses. For example, isolated long fiber carbon nanotubes are more inflammogenic in the outer regions of the lung than non‐fibrous nanotubes.

Large surface areas of nanomaterials increase the contact area with the biological environment and their chemical reactivity.

Surface coating material of the nanomaterials can affect biological functions. It was found that toxicity of nanomaterials that have the same metallic core could be predicted by using the properties of the coating material.

Other physicochemical properties such as electrostatic interactions between nanomaterials and biological targets can influence toxicity.

There are several mechanisms by which nanoparticles induce toxicity:

interaction and binding of the nanoparticle's surface with a biological environment (e.g., protein or cells),100

cellular entry: nanoparticles potential to enter cells,101

release of ions from the surface: ionic forms of metals can be more active,102, 103 and

generation of reactive oxygen species (ROS): overproduction of ROS can cause oxidative stress and inflammation, which disrupt normal biological functions and damage DNA and proteins.104

Quantitative structure nanotoxicity relationships (QSNR)100 (also called nano‐QSAR105) are QSAR models that use nanomaterial‐specific descriptors such as size, shape, surface area, relaxivities, solubility, zeta potential, corona composition, biodistribution, bioavailability, and surface charge in addition to structural and physicochemical properties.100, 106 However, toxicity of nanomaterials is affected by properties of their core and surface coating material, and different combinations of core and surface compositions may cause different effects. Moreover, the coating material can be modified by the biological environment.100

An example of a linear QSNR model for predicting EC50 (effective concentration for 50% enzyme activity inhibition) for organo‐coated silver nanoparticles using size and surface charge is available in ref. 107. Additionally, QSNR models that use quantum‐chemical descriptors have been developed to predict binding affinity for a set of fullerene‐C60 derivatives108 and predict cytotoxicity of metal oxide nanoparticles.109 Another model to predict metal oxide nanoparticles based on random forests classification is explained in ref. 103. Case studies on developing QSNR models are available in ref. 110.

There are some challenges in modeling nanomaterials. Nanomaterials can be composed of organic, metal, metal oxides, silica or carbon‐based nanoparticles. Therefore, it is difficult to compile datasets of congeneric nanomaterials.105 Additionally, nanomaterials are structurally diverse and act upon different MoA. However, QSNRs (similarly to QSARs) assume a common MoA for the modeled chemicals.105 The scarcity of nanomaterial experimental data and descriptors100 affects progress in this field.

However, a nano‐read‐across model (read‐across for nanomaterials) can be used because read‐across methods do not require large datasets to generate groups of sufficiently similar nanomaterials.105 An application of a nano‐read‐across model to predict cytotoxicity of metal oxide nanoparticles is provided in ref. 105.

Additionally, PBPK models may predict biodistribution, which is essential for assessing toxicity of nanomaterials.106 Customized PBPK models should be developed for nanomaterials because ADME properties, some physiological processes of nanomaterials, and transportation mechanisms are different from those of small molecules. Moreover, there are some processes that are not involved with small molecules such physical property changes, enzymatic degradation, cellular recognition, and internalization and opsonization in the blood.106 PBPK may include nanomaterial‐related descriptors such as traffic within tissues and cells, interaction with blood and tissue cells, tissue/blood partition coefficients, tissue concentration, and permeability through membranes.106 An example of a kinetic model to describe mechanisms of releasing silver ions of citrate‐coated silver nanoparticles in aqueous environments is available in ref. 102.

Mixtures

Toxicity of chemicals is affected by interactions with other chemicals. For example, mixtures may exhibit adverse effects at NOAEL doses of each chemical separately.111 Assessing toxicity of chemicals separately may underestimate or overlook the adverse effects of mixtures.111 For instance, it was found that toxicity of lead metal increases with the co‐administration of higher levels of other metals. Therefore, ‘cumulative risk assessment’ was developed to study toxicity of mixtures.112

However, there is lack of experimental datasets for toxicity of mixtures due to a large number of different combinations of chemicals,111 exposure patterns, and complex interactions. It is impossible to test all combinations of these factors. Furthermore, predictive models must address concurrent and sequential exposure to mixtures. However, a recently developed database by NoMiracle (Novel Methods for Integrated Risk Assessment of Cumulative Stressors in Europe) contains mixtures’ toxicity datasets for eco‐toxicological test species and human cell lines.112

Methods for single chemicals may not be applicable for mixtures due to difficulty in determining the combined effect.113 For example, dose–response models for mixtures vary depending on the dose ratios of chemicals in the mixture.112 Additionally, co‐administration of chemicals may alter their ADME properties,111 which should be taken into consideration when developing PBPK models for mixtures.113

Bliss generated dose–mortality curves for mixtures by changing the doses of mixtures while preserving the constituents’ ratios. He identified three categories of the joint action of mixtures.114 The first category is independent joint action. Chemicals act independently and have different modes of action. The combined effect is calculated using the effects of constituents and their interactions.114

The second category is similar joint action. Chemicals act independently and have similar MoA. The combined effect is calculated using the dose–mortality curves of constituents. This category assumes that an ingredient in the mixture can be substituted for any proportion of another ingredient without changing the combined effect.114

The third category includes synergistic and antagonistic actions. Toxicity of synergistic action is greater than that of constituents, while antagonistic action has lower toxicity than that of constituents. Synergistic effects depend upon the proportion of constituents in mixtures unlike the first two categories in which chemicals act independently, and therefore, their proportions do not alter their combined effect. Bliss developed two models (also called independent action (IA)115) to analyze synergism as explained in ref. 114.

Another additive model is the concentration addition (CA) (also called dose addition) that was developed for chemicals with a similar MoA.115 The concentration of mixtures that produces a certain effect is calculated using the proportion of a constituent and its concentration that produces the same effect.115

IA and CA models have been criticized for being ineffective for chemicals that have high potency (dose to produce a given effect) but low efficacy (maximum effect).115 Therefore, a generalized concentration addition (GCA) model was developed to address these shortcomings.116 GCA calculates the combined effect of a mixture using the potency and efficacy of the mixture's constituents.115 A comparison of CA, IA, and GCA models is available in ref. 115.

Process‐based models, however, are mechanistic models that usually use dynamic energy budgets theory such that the combined effect is calculated using the effects of constituents in addition to exposure time, toxicokinetic, and biological parameters, which allow for extrapolation between different species, chemicals, or exposure duration.117 A discussion and an example of process‐based models are available in117 and,118 respectively.

Another mechanistic model is the receptor‐oriented model, which is based on the premise that the toxicity of mixtures is caused by many chains of reactions that converge at the exposed receptor (i.e., an individual or population).112 A discussion of receptor‐oriented models is available in ref. 112. Other methods to assess toxicity of mixtures include numerical additive models such as hazard index, point of departure index, margin of exposure and cumulative risk index113; chemical interaction models such as the interaction‐based hazard index and isobole method113; and statistical models such as tree‐based clustering and weighted quartile score regression.111

RECOMMENDATIONS

The previous sections discuss advantages and limitations of different models. The flow chart in Figure 9 provides a practical guideline for choosing a method for certain types of features and toxicity endpoints, and a summary of the methods is depicted in Figure 10. In this section, we provide some recommendations for developing in silico toxicity prediction models.

Figure 9

Flowchart of in silico prediction models.

Figure 10

Summary of methods to predict toxicity of chemicals.

Datasets Quality

Large datasets from one resource must be used to achieve statistical significance and structural diversity.92, 119 While large datasets may not be readily available, it is possible to gather data from various public resources120: prepared datasets suitable for QSAR modeling such as DSSTox (www.epa.gov/comptox/dsstox/), databases that contain toxicity data such as ECOTOX (http://cfpub.epa.gov/ecotox/), global search engines that retrieve data from several databases such as ToxNet (http://toxnet.nlm.nih.gov/), and scientific literature. A discussion of advantages and limitations of these resources in addition to case studies on obtaining toxicity data are available in Ref 120. Examples and descriptions of many toxicity databases are provided in Refs 6, 8, 9, 10, 11.

However, several factors must be considered when using data from multiple resources. The experimental design must be identical or compatible.19 Additionally, the datasets should preferably be balanced (equal number of active and inactive chemicals) as some statistical methods are sensitive to unbalanced datasets.19 This limitation, however, is not absolute as there are methods that can handle imbalanced data. Moreover, experimental and computed values of the same molecular descriptor should not be mixed, and calculated molecular descriptors must be generated using the same tools.19 In addition, redundant records and outliers should be removed if necessary.119

The readers are referred to a recent review121 that discusses formal scoring methods to assess quality of datasets; chemical, biological, and endpoint‐specific factors in data variability; a practical checklist to guide the process of data assessment; data integration; and many other considerations of dataset quality.

Data Transformation

In some cases, biological data must be transformed. Continuous toxicity endpoints can be expressed in logarithmic scales to avoid statistical problems when using models such as regression.19 Additionally, continuous endpoints, which come from different resources or are generated using different experimental procedures, can be transformed to categorical values (the number of categories are problem dependent).19 Features must be normalized when using statistical similarity measures because these measures are biased toward features that have large values.47 Moreover, frequency descriptors (e.g., frequency of atoms) must be transformed to binary descriptors if the frequencies are not high enough to achieve reliable coefficients when deriving regression models.19 In addition, values must be expressed in the same units.119

Feature Selection

The number of descriptors should be as small as possible.92 Feature selection is necessary to ensure that descriptors are independent, orthogonal, uncorrelated, and non‐redundant. This will increase chances that the models generated are closer to the optimal as this should remove correlation between descriptors. The descriptors should be meaningful to simplify model interpretation. Some methods for feature selection were previously explained in the QSAR sub‐section. Additionally, it is preferable to use features that are suitable for toxicity endpoints (e.g., using environmental‐related descriptors such as acidity of the water for ecotoxicity), chemicals type (e.g., using nanomaterials‐specific features for nanomaterials), and model types (e.g., using PK descriptors for PBPK models).

Statistical Models

Simple accurate models are preferred to the complex ones. However, more sophisticated algorithms can be used to achieve better performance as long as interpretability of the model is preserved and overfitting is avoided. Interpretability of the models depends on the meaning of the descriptors and the relationship between them.90

There are two types of models12: correlative models, which use statistical algorithms that relate molecular descriptors to toxicity endpoints without providing a mechanistic insight, and mechanistic models, which explain the underlying mechanism of toxicity. Mechanistic models can be generated by statistical algorithms or designed by human experts.

However, models that are not mechanistically interpretable are still useful if they produce meaningful predictions, which is usually the case if they were correctly developed from large datasets.19 Generally, models must meet users’ expectations. If users expect models that provide mechanistic insight of toxicity, developers must try to meet these expectations.

Overfitting

Overfitting occurs when a model fits the training data extremely but cannot describe the external data well. It can be detected when models achieve high performance on training data but low performance on new data. Overfitting can be caused when modelers develop models that cover all cases in the training set or use complex models or the model includes more descriptors than necessary or wrong descriptors.12 It is customary to have at least five chemicals per descriptor. Otherwise, an overfit model may become meaningless, be overly complicated, find random correlations, or have a limited applicability domain ref. 97.

There are two main statistical parameters to assess goodness of fit of models (e.g., QSARs generated using regression): coefficient of multiple determination (R 2) and standard error of estimate (s).92, 122 R 2 estimates how successful the regression line is in explaining variation in the biological data (e.g., measured response), while s estimates deviation of the measured response from the regression line.92 A well‐fit model should have large values of R 2 (close to 1) but low values of s (close to 0).92, 122 However, an overfitted model can be generated if its statistical fit is greater than the experimental variation (caused by errors in measured response or calculated descriptors).122 It is possible to increase R 2 by adding more (possibly irrelevant) molecular descriptors. This can be avoided by using adjusted R 2 (R 2 (adj)). Unlike R, 2 R 2 (adj) does not increase when adding irrelevant descriptors. Mathematical definitions of these and other statistical parameters for fitting are explained in,92 and a discussion of sources of experimental variability is available in ref. 122.

Applicability Domain

The applicability domain (AD) is ‘a theoretical region in physicochemical space (the response and chemical structure space) for which a QSAR model should make predictions with a given reliability.’92 The AD determines types of molecules and toxicity endpoints to which the model can be applied. For example, global and local QSARs have large and small ADs, respectively.6 Interpolation within the domain is more reliable than extrapolation outside the domain.19 AD determines uncertainty in predicted toxicity by measuring the difference between the new molecules and the training set.97 The model may not be reliable or perform well if the new molecule is sufficiently different than the training set.6, 90

A stepwise approach to defining ADs was developed by Dimitrov et al.123 and further expanded by Hewitt and Ellison.124 Five categories of AD definition approaches were explained and compared in124 descriptor, structural fragments, structural similarity, mechanistic, and metabolism‐based approaches.

The descriptor‐based category includes four types of methods: range, geometric, distance, and probability density‐based methods, which were compared in ref. 124. Range‐based methods require determining the range (i.e., minimum and maximum) of each descriptor. However, geometric methods such as convex hull define the smallest convex area containing the whole training set. Probability density distribution‐based methods identify probability density of the dataset and then determine the highest density region. Distance‐based methods use distance measures (e.g., Euclidean, Hamming, or Tanimoto). Then, the average distance to k‐nearest neighbors or the distance to the centroid is taken to determine how similar/dissimilar the new chemical is from the training set. A threshold can be used to determine if the new chemical is within the AD. A comparison of different distance measures to define AD is available in ref. 97. Another application using cluster analysis is illustrated in ref. 125.

Structural fragment‐based ADs require that all structural fragments in the new chemical be present in the training set.124 However, the structural similarity‐based ADs determine structural similarity of the whole chemical rather than fragments.124 Mechanistic ADs have been discussed in refs. 123, 124, 125, 126. Another AD method using simplified molecular‐input line entry system (SMILES) attributes is illustrated in ref. 127. Moreover, examples of tools for defining ADs (e.g., The OECD QSAR Toolbox and ToxTree, etc.) are explained in ref 124.

Model Evaluation

Cross‐validation on training sets is insufficient for a realistic evaluation because models are usually applied to new chemicals. Models must be evaluated using testing sets that are large enough to achieve statistical significance, diverse, within the AD, and different than the training sets. Evaluating models’ accuracy is discussed in ref. 92.

Model Application

Model development processes and AD must be transparent to users to ensure appropriate application of the models. For example, molecular descriptors of new chemicals must be calculated and expressed in the same units as in the training set, and the new data must be transformed as the training set. Additionally, the output of the models must be clearly explained. Moreover, models preferably should provide a rationale (e.g., a list of the rules that were used to form the prediction). Also, models should provide a confidence score of the prediction. This will help users decide whether to accept the predictions of the models. Consensus can be used to determine reliability of the prediction as explained in ref. 15. For example, if a chemical was predicted to be carcinogenic by several tools, there might be higher confidence in the prediction than if it was predicted by a single tool. However, if the tools provide conflicting predictions, then it is necessary to resolve such conflicting situations, which is usually done by an additional machine learning model.

TWENTY‐FIRST CENTURY TOXICOLOGY

The phrase ‘21st century toxicology’ (Tox‐21c128) refers to ‘the transformation underway in the tools and approaches used to evaluate chemical substances for possible effects on human health.129’ Tox‐21c focuses on toxicity pathways,130 mechanisms/modes of action, and adverse outcome pathways (AOP)131 in humans. Another related concept is the 3Rs (replace, reduce, and refine), which was proposed in 1958.128 Tox‐21c and 3Rs overlap, but a discussion of their differences is presented in ref. 128.

Several strategies have been proposed to implement Tox‐21c. In 2004, the National Toxicology Program (NTP) published its report A National Toxicology Program for the 21st century, which aims ‘to support the evolution of toxicology from a predominantly observational science at the level of disease‐specific models to a predominantly predictive science focused upon a broad inclusion of target‐specific, mechanism‐based, biological observations.132’ In 2007, the National Research Council (NRC) published another report Toxicity Testing in the 21st Century: a Vision and a Strategy, which proposed using computational methods to decrease the number of tested animals, make toxicity testing more relevant to humans by using human cells, and make toxicity testing cheaper and faster.133

Several approaches have been developed to implement these visions. Integrated Testing Strategy (ITS) aims to combine information from testing and non‐testing methods.134 Another approach is Pathways of Toxicity (PoT), which aims to produce a comprehensive list of all human PoT.135 However, Integrated Approaches to Testing and Assessment (IATA) focuses on using hypotheses to prioritize chemicals for testing.131

Several projects implement the above approaches. For example, ReProTect is an integrated testing project specializing in reproductive toxicity.136 The ToxCast (Toxicity Forecaster) project uses quantitative high throughput screening to test chemical toxicity on many biological pathways,137 and it is part of the Tox21 project that follows NTP and NRC guidelines.138 Additionally, the Human Toxicology Project Consortium implements the NRC guidelines and focuses on implementing AOP .139

Conclusion

The field of in silico toxicology has been in a continuous development through the introduction of new methods, improvement of the existing ones, or discarding of some of them. Unfortunately, a method that is suited for certain types of toxicity endpoints or chemicals may not work properly (or not work at all) for others. If used correctly, in silico tools can be very effective in assessing chemicals’ toxicity. Therefore, to ensure accurate and effective application of in silico models, it is necessary to (1) understand the methods’ strengths, limitations, scope of application, and interpretation; (2) choose the most effective method for the problem at hand; and (3) customize these methods for each problem if necessary. Users of toxicity prediction models can follow these three steps only if the data and processes to develop the model are transparent, applicability domains are well defined, the outputs of the models are clearly explained, and models are simplified.

The Tox‐21c stresses replacing animal testing with human‐relevant testing methods, either in vitro or in silico. With the increasing number and variety of alternative testing methods, it is necessary to apply strategies (e.g., ITS) to intelligently combine and use this information for toxicity assessment and decision making. Clearly, in silico toxicology is a useful component of the toxicity assessment process. Looking to the future, computational methods are likely to expand to include models for special and new types of toxicity endpoints and chemicals, provide insight into toxicological pathways, combine and compare results from different models, customize models to meet users’ expectations, and refine models as new data become available.

Table S1. Different types of PK and PD models.

Table S2. Different method to determine uncertainty factors.

Table S3. Different types of models and molecular descriptors in the QSAR family.

Click here for additional data file.