Skin Sensitization: Modeling Based on Skin Metabolism Simulation and Formation of Protein Conjugates

Skin Sensitization Data

The model was derived from a data set compiled from chemicals tested in the LLNA (185 chemicals), GPMT (307 chemicals), as well as from the BgVV list (248 chemicals).

The LLNA was initially developed as an alternative approach to hazard identification (Kimber and Basketter 1992; Kimber et al. 1994, 2002; Dearman, Basketter, and Kimber 1999; Gerberick et al. 2000; Basketter et al. 2002). In the local lymph node assay, the response of the immune system to the topical application of chemicals is measured by the lymph node cell proliferation. Typically, groups of mice are exposed to various concentrations of a test chemical or to the vehicle alone for 3 consecutive days on the dorsum of both ears. The practice is to select three consecutive suitable concentrations of the chemical in ideally a standard vehicle, looking for the highest possible concentration, whilst avoiding the unacceptable dermal trauma or systemic toxicity. Five days after the initiation of the exposure, mice are injected intravenously with [³H]thymidine. Five hours later, the animals are sacrificed and for each dose group, a stimulation index (SI) is determined by measuring the increase in [³H]thymidine incorporation in the lymph nodes relative to the vehicle-treated control group. Currently, the preference is to estimate the concentration of a chemical required to generate a three-fold greater increase in radioisotope incorporation as compared to the control, known as the EC3 value. Although the current European legislation requires dichotomous classification of skin sensitizers, based on EC3 values, different potency categories of contact allergens can be differentiated. Table 1 illustrates the LLNA classification scheme based on 10-fold variation of EC3 values (ECETOC 2003).

The guinea pig maximization test was specifically designed for hazard identification and is not well suited to potency estimation (ECETOC 2000). First, a pilot study is performed to determine the maximum nonirritating dose. Test and control groups contain 15 to 20 animals. Testing lasts approximately 5 weeks and includes an induction phase, rest period, and 24-h topical challenge. Red and swollen sites on test animals indicate sensitization. Sensitization potential is characterized based on the percentage of animals exhibiting reactions. The classification of sensitization potential for GPMT is presented in Table 1 (Barratt et al. 1994; ECETOC 2003). Classification schemes for the LLNA and both guinea pig test methods have recently been proposed by ECETOC (2003).

At the German Federal Institute for Health Protection of Consumers and Veterinary Medicine (BgVV), a group of experts including dermatologists from universities and representatives of chemical industry and from regulatory authorities was established in 1985 (Schlede et al. 2003). The aim of the project was to bring together a broad range of expertise to collect and evaluate data from the literature on substances with documented contact allergic properties in humans and animal experiments. The evaluation results for selected chemicals were published as a loose-leaf book (Kayser and Schlede 2001). In this publication, chemicals were listed as category A (significant contact allergen), category B (solid-based indication for a contact allergenic potential), or category C (contact allergen with insignificant or questionable contact allergenic potential). The definition of these categories is given in Table 1. The authors realize that the BgVV classification scheme differs from the pure hazard identification schemes applied in the LLNA or guinea pig testing because classification according the BgVV can also take into account prevalence in the population or clinical relevance.

Ideally modeling would have been carried out using one experimental protocol, preferentially the LLNA. Because the goal was to derive a QSAR as widely applicable as possible, the different training set test protocols were combined to expand the chemical diversity and the sensitization data reassessed using a unified potency scale. This did mean that some of the Category C BgVV data needed to be re-reviewed by the authors and reclassified from weak sensitizers to nonsensitizers. This unified classification scheme (Table 2) mirrors that for the LLNA and the EC3 as presented in Table 1. The first three categories of the EC3 scheme (extreme, strong, and moderate sensitizers) were merged into a single category (S = significant sensitizer) including chemicals with EC3 less than 10%. The other two categories coincide with the EC3 scheme: weak (W) for 10% ≤ EC3 ≤ 100% and nonsensitizers if stimulation index did not exceed threefold effect compared to the controls. The present unified classification scheme could easily be expanded for any subsequent refinements to the QSAR model, e.g., to discriminate between extreme, strong, and moderate sensitizers. Table 2 illustrates the relationship of the outlined three- category scheme with those from the other two data sets (GPMT and BgVV).

It should be noted that 17 chemicals exist in all three data sets; 34 chemicals were present in both the LLNA and BgVV data - sets; 44 chemicals were common to both GPMT and BgVV data sets, and 41 chemicals in both the LLNA and GPMT - data sets. The strength of the association between the different data sets was estimated using the adjusted Pearson’s contingency coefficient, C* (Agresti 1996). The values showed - a well-defined correlation between observed sensitization potency of chemicals in different databases—from moderate (C ^* _GPMT–BgVV = 0.48) to significant (C ^* _LLNA–BgVV = 0.67 and C ^* _LLNA–GMPT = 0.61). For chemicals whose potency was assessed by more than one method, the LLNA value was accepted as representative for the structure; where a LLNA test was not available, the GPMT value was taken. Excluding the coincident chemicals reduced the training set to 634 chemicals.

Modeling Approach

The single biodegradation pathway probabilistic scheme used in CATABOL (microbial metabolism simulation model) (Jaworska et al. 2001) was modified in TIMES (Mekenyan et al. 2004) to enable the modelling of multiple biotransformation pathways in the skin. The multipathway scheme was predicated on the condition that a chemical could be metabolically transformed across both the most probable pathway as well as less significant pathways. The abiotic molecular transformations and enzyme-mediated reactions such as oxidative, redox, reductive, hydrolytic, and conjugative reactions, and reactions with skin proteins, were mimicked by a hierarchically ordered list of principal molecular transformations.

Each molecular transformation comprises a parent submolecular fragment, transformation products, inhibiting fragment (masks), and a reactivity SAR/QSAR model. The masks act as reaction inhibitors; thus, a fragment assigned as a mask attached to a target subfragment could prevent the transformation on the parent chemical from occurring. The presence of groups that can promote or inhibit enzymatic reactions significantly increases the number of principal transformations. A probability of occurrence is assigned to each transformation to determine its hierarchy in the transformation list.

A parent molecule with the source fragment is matched against a list of hierarchically ordered transformations, starting with those having the highest probability. This effectively produces a set of first level metabolites. Each of these derived metabolites is then submitted to the same list of hierarchically ordered transformations, to produce a second level of metabolites. The procedure is continued until a constraint for metabolism propagation is satisfied [e.g., low probability of obtaining a metabolite; application of Phase II reaction (glucuronidation, sulfate conjugation, etc.); protein conjugation or formation of a hydrophilic metabolite]. The metabolic profile map that is generated is then used to calculate both the quantities of metabolites and the probabilities of the predicted metabolites. TIMES is completely integrated with both a QSAR library and a modelling engine to provide an estimate of physicochemical properties and the toxicity of each metabolite. A detailed explanation of the mathematical formalism for the TIMES approach for metabolite simulation and a discussion of the practical applications for forecasting plausible activated metabolites have been published elsewhere (see Mekenyan et al. 2004).

The initial simulation of metabolism by TIMES is based mainly on two-dimensional (2D) chemical structures. Whilst it is a widely accepted hypothesis that skin sensitization involves the covalent binding of the sensitizing chemical to protein in skin, models based only on 2D molecular structure may be insufficient to account for reactivity and, hence, sensitization potential. For this reason, the COmmon REactivity PAttern (COREPA) method (Mekenyan et al. 1997, 1999; Bradbury et al. 2000) was incorporated into the model to analyze certain structural alerts that are distributed between sensitizing and nonsensitizing chemicals in the training set. COREPA is a probabilistic technique for identifying common stereoelectronic (reactivity) patterns of structurally diverse chemicals, which may exert similar or differential biological effects. All energetically reasonable conformers are used to establish conformer distributions across the global and local stereoelectronic descriptors associated with the activity under study. The COREPA model is derived in the form of a decision tree. Its logic boxes consist of decision rules based on the reactivity patterns described by a combination of global descriptors of molecular steric and electronic structure and local reactivity parameters associated with specific alerting groups (Mekenyan et al. 1997, 1999; Bradbury et al. 2000).

Model Construction

There is a long-established connection between the ability of chemicals to react with proteins to form covalently linked conjugates and their skin sensitization potential (Landsteiner and Jacobs 1936; Dupuis and Benezra 1982; Basketter et al. 1995; Leppoittevin et al. 1998). For skin sensitization to occur, once a chemical has penetrated, it must be able to partition into relevant cellular compartments of the epidermis, in order to be sufficiently bioavailable. More importantly, in terms of mechanism, the chemical, or its metabolite, must be sufficiently electrophilic to react covalently with nucleophilic groups on skin protein to produce the complete antigens capable of invoking an immune response. The skin sensitization model was built as a composite of the following submodels:

Skin metabolism simulator: This mimics the metabolic fate of parent chemical controlled by skin enzymes and thus the potential formation of protein adducts with reactive agents. 2D structural information of parent chemicals is used to model metabolism. Metabolic pathways are generated based on a set of hierarchically ordered principal transformations including spontaneous reactions, enzyme-catalyzed phase I and phase II drug metabolism reactions, and reactions with protein nucleophiles. The formation of macromolecular immunogens was used to identify probable structural alerts in parent chemicals or their metabolites.

Skin Metabolism Model

The core of the skin metabolism model is a set of 203 principal transformations. These were defined on the basis of empirical and theoretical knowledge and hence peer-reviewed by human experts. A representative subset of these transformations, as well as their reliability, probability, and inhibiting masks are outlined in Table 3. Due to the limited quantitative skin metabolism data reported for chemicals, many transformation reliabilities were assigned on the basis of available information or expert knowledge in one of five confidence levels. The highest confidence level of reliability (R _tr = 1.00) was assigned when a reference source existed which supported enzymatic transformation or when this was an abiotic transformation commonly known. A reliability weight of 0.75 was assigned when no reference source existed which supported transformation; but the enzymatic reaction was well known. A reliability weight of 0.50 was assigned when the transformation was expected but the enzymatic system causing it was unknown. A reliability level of 0.25 was assigned to a transformation when expert judgement supported its feasibility. Finally, a confidence level of reliability of 0 was assigned to transformations where no knowledge of feasibility. The reliability of metabolites resulting from the simulated skin metabolism maps was evaluated by the product of the transformation reliabilities across the pathway reaching these metabolites (Mekenyan et al. 2004).

The principal transformations were separated into two major classes: non–rate-determining and rate-determining reactions. Non–rate-determining included abiotic and biologically mediated transformations, which occur at very high rates. They also included nine non–rate-determining reactions predominantly hydrolysis of salts. Eighty-three other transformations simulating the reactions of highly reactive groups and intermediates such as quinones, α,β-diketones, oxiranes, acyl halides, thiocarboxylic acids, hydroperoxides, nitrenes, geminal diols, etc., with cell proteins were also included. All non–rate-determining reactions and covalent interaction between allergens and cutaneous proteins were assigned with the highest possible probability of 1.00. Interactions with proteins were grouped into three types: reactions causing a significant sensitization effect [extreme, strong, moderate: (S)], reactions causing a weak (W) sensitization effect, and interactions with proteins requiring a QSAR(s) model to be utilized in order to quantify the likely potency of a sensitization-alerting group.

Phase I reactions were modelled through 102 transformations. These transformations included amongst others, C-hydroxylation, ester hydrolysis, dealkylation, oxidation, deamination, and oxirane hydration. Decomposition of geminal derivatives and keto-enol tautomerization were also modeled in this group. Phase II reactions such as glucuronidation and sulfate conjugation were modelled using nine principal transformations. The probability of phase I and phase II reactions were determined by expert knowledge.

The set of hierarchically ordered principal transformations within TIMES effectively simulate metabolism by generating plausible and likely metabolic pathways. The established hierarchy between the transformations and hydrophobicity effectively control the propagation of the generated metabolic pathways. The process of metabolic tree generation continues until user-defined thresholds for probability of obtaining of metabolites and their hydrophobicity (log K _ow ) were reached. In this study the thresholds for probability of obtaining of metabolites were set at 10⁻⁵. The thresholds for hydrophobicity were not taken into account. A maximum of 10 rival transformations at each level was applied. The propagation of the metabolic tree was interrupted if a phase II or protein conjugation reaction occurred. Figure 1 illustrates the predicted skin metabolism of isoeugenol (CAS no. 97-54-1).

Isoeugenol is shown to undergo demethylation to give rise to the corresponding catechol, followed by oxidation to the o-quinone (Barratt and Basketter 1992); or formation of semiquinone free radicals or to undergo sulfate or glucuronide conjugation (phase II reactions). The reactive intermediates are quinone and semiquinone structures, which are capable of reacting covalently with skin proteins.

Modeling skin metabolism led to the following percentages of correctly classified S sensitizers: LLNA 78%, GPMT 81%, and BgVV 82%. It is worth mentioning that these percentages were lower than 50% when transformations modeling metabolic activation were not considered. Good predictions were also obtained for nonsensitizers (N): 68%, 70%, and 57%, respectively, for the LLNA, GPMT, and BgVV datasets. The percentage of correctly classified weak (W) sensitizers was 26%, 22%, and 39%, respectively, in the three data sets. The poor classifications for weak sensitizing chemicals may be attributed primarily to the use of 2D structural information and some weakness in the underlying biological data set. The majority of chemicals with structural alerts were classified as significant (S) sensitizers.

QSAR Models of Reactivity

Pattern recognition models were derived for ketones; sulfuric, sulfonic, or phosphonic acid esters; aldehydes; and conjugated double bonds containing O, S, or N atoms in order to improve the classification of chemicals (containing these functional groups) on the basis of skin metabolism simulation. Models were derived at a probability threshold of 0.6 (Mekenyan et al. 1997, 1999; Bradbury et al. 2000). These four specific alerting groups and the number of experimentally observed S, W, and N sensitizers are presented in Table 4.

Due to limited data, chemicals having only one of the above mentioned four alerts from the three datasets were combined to form the corresponding training sets. For some of alerts, such as ketones, acid esters and aldehydes, QSAR models were only trained to discriminate between S and N (or between S and W) sensitizers as there were insufficient data for weak (or nonsensitizing) chemicals available.

For example, the COREPA decision model for classifying chemicals with conjugated double bonds containing either O, S, or N atoms as strong sensitizers uses as discriminating parameters in the logic boxes of the decision tree the lowest unoccupied molecular orbital (E _LUMO ), HOMO-LUMO energy gap (E _GAP ), and molecular weight (MW). The second decision tree discriminates weak from non-sensitizers using E _GAP , E _LUMO , and log K _OW as key parameters. The COREPA models for ketones, acid esters, and aldehydes contain only one logic box. The key parameters for ketones were log K _OW and the net atomic charge on the ketone C-site, whilst for sulfuric, sulfonic, and phosphonic acid esters the discriminating parameters were electronegativity (EN) and log K _OW , and for aldehydes the key descriptors were MW, log K _OW , and the net atomic charge on the aldehyde C-site. The sensitivities of the COREPA models varied from 60% to 100% across different sensitization categories.

Integration of Skin Metabolism and COREPA Models

QSAR models derived for identifying skin sensitization potency for some alerting groups were used to reduce misclassifications based on skin metabolism only. During metabolism simulation, chemicals having select alerts (such as ketones; sulfuric, sulfonic, or phosphonic acid esters; aldehydes; and conjugated double bonds containing O, S, or N atoms) were screened by the COREPA models to determine the level of activation of the alerting groups and subsequently their ability to interact with protein nucleophiles. The results of the combined application of skin metabolism and partial COREPA models in predicting the skin sensitization potential (as S, W, or N) of the chemicals in the LLNA, PGMT, and BgVV datasets are summarized in Figure 2.

Model predictions for the group of W sensitizers were most significantly improved with the percentage of chemicals correctly classified increased to 42%. As can be seen from Figure 2, COREPA models do not significantly improve the identification of S and N sensitizers. The present QSAR model has some limitations with chemicals such as acrylates, fused acyclic hydrocarbons, and multifunctionalized chemicals containing phenolic, ether, ester, lactone, azolactone, or isothiocyanate moieties in the same molecule. In addition, metals are limited primarily to sodium (predominately as sulfonate salts or acetate salts); lead as lead acetate; mercury as methylmercury ion; and zinc as its thiocarbamate salt.

Applicability Domain of the Skin Sensitization Model

The number and breadth of chemicals underpinning the model are limited and not universally applicable. Consequently, it is critical to establish whether the model is applicable or appropriate for a given new query chemical submitted for screening. The basic principles of the theory of chemical structure introduce the concept of the mutual influence amongst the atoms in a molecule and this can be used to help define the applicability domain of the model. It is common knowledge that besides the first neighbors, all other atoms in a molecule affect each other. For example, a chlorine atom attached to a carbon atom differs substantially depending on whether that carbon atom belongs to an alkane, alkene, or aromatic-ring fragment.

Based on the concept of mutual influence between the atoms in a molecule, the structural domain of an atom participating in a number of molecules, ions, or radicals can be defined by the set of all atom-centered fragments containing this atom and its first, second, etc., neighbors. The influencing effect of neighbors decreases with the increase of their distance to a specified atom; it also depends on the atomic type of the neighbors. For example, conjugated systems and aromatic-ring systems may transmit the electron withdrawing or donating effects of an atom much further. The following rules are proposed to reflect the effect of different neighbors on the specified atom:

Hydrogen atoms are treated as an inherent part of the atom to which they are bound; so effectively, they are ignored and are considered to be a characteristic of the atoms to which they are bound.

These rules were used to determine the structural domain of the training set and applicability domain of the QSAR model. In the first case, the set of substructures is extracted from all chemicals that belong to the training set. In the second case, the fragments are extracted among those chemicals from the training set for which the QSAR model correctly predicts the modeled end point. The extracted atom-centered fragments can be used to determine whether a certain chemical belongs to the training set domain and applicability domain of the model or not. The latter is a part of the training set domain where one could expect a good performance of the model.

In the present work, the chemicals from the training set that were correctly classified as significant, weak, and nonsensitizers were used to extract the domain of the model. The first or second neighbors of the atom-centered fragments were used to define structural similarity between chemicals for which the model worked correctly. The above approach was applied to assess the structural domain of the training set and QSAR model for skin sensitization. Table 5 summarizes the statistics of the obtained model domain. The categorization of the chemicals from the training set identified all correctly predicted chemicals as belonging to the model domain. Additional chemicals were incorrectly classified as a part of the model domain. This misclassification may potentially be explained by such factors as experimental error, model inadequacy, or the empirical approach used to determine the structural domain. Regardless of the empirical nature of the outlined approach, the goodness of the determined model domain was very high if the atom-centered fragments included the second neighbors. As can be expected, determination of the model domain accounting for first neighbors only resulted in a significant increase of chemicals that were incorrectly classified as elements of the model domain.

External Model Validation

From the various validation approaches available, external validation was applied. This is the most demanding way of testing the validity of a model. The approach consists of making predictions for an independent set of data not used in the model training (Eriksson, Jonsson, and Berglind 1993).

A set of 96 chemicals tested for skin sensitization by LLNA (47 chemicals) and GPMT (49 chemicals) tests were collected randomly from unpublished and published data sets (Cronin and Basketter 1994). These chemicals were not used in the model building process. The chemicals identified as in the model domain and their predictions for their skin sensitization potentials are presented in Table 6. As can be seen from the table, the model predicts the external data compounds reasonably well if the model domain is determined by accounting for the second neighbors of centered atoms. In this case, the correctness of predictions is 87% as compared with 71% if the model domain is determined by accounting for first neighbors only. It should be mentioned that ignoring the model domain reduces the predictability of the model to only 52%.