Modeling drug-induced liver injury

This case study is about training machine learning models for drug-induced liver injury (DILI) prediction. DILI is a significant concern in drug development, as it can lead to severe health issues and drug withdrawals and thus to financial losses.

By conducting this case study, we aim to explore the application of machine learning techniques in predicting DILI based on molecular structures, mostly represented as SMILES strings. We separate our work into the study of different molecular pathways that can lead to DILI, such as BSEP inhibition and mitochondrial membrane potential disruption (more on these pathways in the Background section).

We make the code for training these models opensource and clearly indicate the license of the data for each model. With this, we hope that scientists will be able to use our suite of models to derisk drug development. The different pathway models are individually packaged and versioned, we further report exhaustive and transparent evaluation metrics for each version.

Our motivation is fourfold. First, we aim to demonstrate our expertise at producing reproducible, packaged and versioned artifacts of research. Second, this research task allowed us to deep dive into concrete and important biological topics on which we have learned a lot. Third, our goal is to create a fully opensource alternative library to existing commercial ones. Finally, we enjoy solving problems, and we are committed to producing public and reproducible artifacts of the best quality.

In cellular biology, one often refers to the concept of "pathways", which are sequences of actions among molecules in a cell that lead to a certain product or change in the cell (metabolic pathways, signaling pathways). DILI is associated with several key pathways, including those involved in drug metabolism (mostly via CYP enzymes), drug transport, DNA damage, mitochondrial membrane disruption, etc.

At the time of writing, we have packaged three of them: the BSEP inhibition pathway (which we will call BSEP), the mitochondrial membrane potential disruption pathway (OXPHOS), and the NRF2 oxidative stress pathway (NRF2).

Since the Swiss physician Paracelsus, it is known that "the dose makes the poison". In this case study however, we will not consider the dose of the drugs. Our virtual pathways are binary classifiers, and what a positive means depends on the pathway: each outputs a probability of being associated with DILI as the authors of its dataset defined it. In most cases the labels are binary in the datasets themselves. Where a dataset provides continuous values instead — IC50 values for BSEP inhibition, for example — we regress those values and provide them as a first output, then threshold them consistently with the literature to give the probability alongside.

Mitochondrial membrane potential disruption

As a well-known sentence says, "mitochondria is the powerhouse of the cell". Mitochondria are organelles that produce energy for the cell in the form of ATP. In order to produce ATP, mitochondria need to manipulate oxygen and electrons, which can lead to the production of reactive oxygen species. Because of their reactivity, reactive oxygen species can damage the cell and lead to cell death if they happen to leak out of the mitochondria. Hence mitochondrial membrane disruption is a key pathway that can lead to cell death and DILI. It often happens specifically in the liver because it is the organ that metabolizes drugs and is therefore exposed to high concentrations of potentially reactive molecules.

BSEP inhibition

Bile salt export pump (BSEP) is a protein that is responsible for transporting bile salts from the liver to the bile ducts. Bile salts are important for digestion and absorption of fats. If the pump is inhibited, bile salts can accumulate in the liver and cause cholestatic liver injury. BSEP inhibition is a well-known pathway that can lead to DILI, and it is often tested in drug development.

NRF2 oxidative stress

This pathway acts as a detection mechanism of reactive oxygen species in the cell. We find the mechanism of NRF2 activation to be particularly interesting. An enzyme called KEAP1 (Kelch-like ECH-associated protein 1) is responsible for regulating the activity of NRF2. Under normal conditions, KEAP1 binds to NRF2 and promotes its degradation, preventing it from activating antioxidant genes. However, when reactive oxygen species are present, they modify KEAP1 and prevent it from binding to NRF2. This allows NRF2 to accumulate and translocate to the nucleus, where it can activate the expression of antioxidant genes.

The virtual pathways (models for each pathway) rely on the structure of the molecules, which is represented as SMILES strings (but we also process these strings using various cheminformatics techniques to provide more specific features to our models). Therefore, models will be mostly blind to the pharmacological properties of the drugs, and will not be able to predict DILI for drugs that for example inhibit a very specific enzyme involved in a DILI pathway. Because of this, we made sure to select pairs of drugs that are chemically similar and have similar pharmacological properties, and of which the DILI mechanism depends on the chemical information.

Several of the compounds below are old, well-studied drugs that public screening libraries are full of. Hence, before analyzing our predictions themselves, we measured how distant each compound was from the training dataset of the models. More specifically, we looked at four indicators. First we looked whether the compound was effectively in the table or not. Then we count the number of training compounds sharing its Murcko scaffold (note that our split between train and val is based on the scaffold). We also define the nearest neighbour as the largest Tanimoto similarity to a training compound over Morgan fingerprints. Finally, the descriptor distance is the distance to the nearest training compound in the standardised descriptor space that is given as an input to the model. We call a compound held out when it is absent from the table and has no near duplicate in it.

Amodiaquine / chloroquine

Amodiaquine (LiverTox) and chloroquine (LiverTox) are both 4-aminoquinoline antimalarials that share a nearly identical core scaffold and mode of antiparasitic action. Yet, only amodiaquine is associated with DILI, which has restricted its clinical use largely to short prophylactic or combination-therapy regimens. Chloroquine, in contrast, is only a rare and idiosyncratic cause of clinically apparent liver injury despite decades of widespread use. The divergence is thought to trace to amodiaquine's extensive hepatic first-pass metabolism (by CYP2C8) into the reactive metabolite desethylamodiaquine, which can then form protein adducts and drive oxidative stress. On the other hand, chloroquine undergoes only minor hepatic metabolism and is mostly excreted unchanged. This pair is therefore well suited to probing the NRF2 pathway.

Distances to the training set

Amodiaquine can be considered held out of the training set of the NRF2 predictor: first, it is effectively absent from the training set. The nearest neighbour is at 0.50, against a median of 0.57 for the table's own compounds, and only 2 training compounds share its scaffold. Chloroquine however is in the table, and the assay labels it inactive.

Stanozolol / nandrolone

Stanozolol (LiverTox) and nandrolone (LiverTox) are both anabolic-androgenic steroids used therapeutically for performance enhancement, but they differ critically in chemical structure at the C17 position: stanozolol is 17-alpha-alkylated, while nandrolone is not. This structural difference explains a divergence in liver safety: 17-alpha-alkylated steroids are consistently linked to a distinctive "bland cholestasis". The alkyl substitution slows hepatic clearance and prolongs hepatic exposure, and 17-alpha-alkylated androgens have also been shown to disturb bile acid transport and antioxidant defenses in the liver. Hence, we studied this pair under the lens of the BSEP inhibition pathway.

Distances to the training set

Nitrofurantoin / trimethoprim

Nitrofurantoin (LiverTox) and trimethoprim (LiverTox) are both commonly prescribed for urinary tract infections and are frequently co-administrated or considered interchangeably in clinical practice. Nitrofurantoin undergoes nitro-reductive metabolism that generates reactive oxygen species through redox cycling, and it has been shown to directly inhibit complex I of the mitochondrial electron transport chain in hepatocytes, depleting glutathione and impairing mitochondrial respiration. Trimethoprim causes hepatotoxicity through a distinct, largely immune-mediated idiosyncratic mechanism (which we do not model for complexity limitations), and is rare enough that some large cohort comparisons found similar background hepatic event rates between trimethoprim alone and the combination. This pair therefore probes the mitochondrial membrane potential disruption (OXPHOS) pathway specifically.

Distances to the training set

Both compounds are effectively present in the training set.

Three pathways are studied. Each of which is an independently versioned Python package with its own license, dataset and model card. All three models take SMILES strings as an input and return the probabilities of a binary classification as an output.

We compare the AUC of the classifiers obtained with a mean of five seeds of a scaffold split. We further report the protocol we used for the split, which is clearly explained in our code. Basically, scaffold-shuffle assigns scaffold groups at random between train and validation sets whereas scaffold-balanced packs each group into the fold that it overfills least (required when the folds are smaller than the size of the groups). We invite the reader to refer to the code for further understanding.

Packagesgithub.com/Sagace-Labs/virtual-pathways-hubOpensource repository of our published pathways.

The pairs are scored on the corresponding pathways. We report the classification. We also report a counter-screen which is a cell-viability read taken from the same library as the primary assay. It indicates whether the compound was killing the cell. We report it when available (BSEP has none, hence the dashes). cLogP is a computed estimate of lipophilicity, which is a confound to DILI caused by BSEP inhibition.

The pairs are correctly classified by our virtual pathways models. The gaps are +0.505 on NRF2, +0.393 on BSEP, +0.906 on OXPHOS.

These numbers are not calibrated probabilities of liver injury in a patient. They represent the probability that a compound is active in a particular cell assay, as that assay's authors defined activity.

One issue with the virtual pathways is that the harm may come from a metabolite of a tested compound. Assays might capture this as a residual signal, when tested on a cell model that correctly expresses the CYP metabolisation pathway. It is not clear however that our models will be able to predict the DILI activity of compounds that do harm through their metabolites, directly via their SMILES string. We tested this directly on five pairs from the medicinal chemistry literature.

A second boundary is immune-mediated injury. It is not directly predictable by a chemistry model and highly varies from individual to individual. Ximelagatran, withdrawn after transaminase elevations linked to HLA-DRB1*07:01, is a standard example. Our virtual pathways score it backwards (OXPHOS: 0.027 against 0.164 for dabigatran).

We intend to develop virtual pathways that predict the metabolites of a given drug and use these predictions as inputs to these more classical virtual pathways. That is a substantially harder problem and it will be the subject of its own case study.

The three pairs above are textbook cases, chosen for simplicity. A possibly harder test is to rank three compounds. We provide the example of the endothelin receptor antagonists. Sitaxentan, bosentan and ambrisentan represent in a sense a gradient of regulatory-flagged DILI. Sitaxentan was withdrawn worldwide in 2010 after fatal hepatotoxicity. Bosentan remains on the market with dose-dependent transaminase elevations and a monitoring requirement. Ambrisentan is in the same drug class, but its liver monitoring requirement was removed in 2011.

Immune-mediated reactions depend on a patient's HLA type, that is not predictable from the chemistry of the drug. Same for an injury caused through an enzyme that we have not modelled. A low score on the panel is therefore not a safety claim in itself.

We model small molecules only. The screens used to train the models were run on small-molecule libraries and we did not find any comparable public resource for biologics. Though we are confident in the possibility of creating virtual pathways for biologics with sufficient datapoints.

Most outputs are binary. Most of the time, we did not design the threshold ourselves. Where the data presented continuous values, we regressed them and reported the estimate alongside (that is the case for BSEP potencies and the endothelin table).

We cherry-picked the pairs to showcase. We swept a much larger set and are first showing the most successful pairs. We present our predictions for other pairs below.

Each virtual pathway is a published package with an immutable version directory holding a manifest, the weights, and a model card. A version also records data provenance. We also declare permissive licenses for the code and make the licenses about the data explicit.

Packagesgithub.com/Sagace-Labs/virtual-pathways-hubOpensource repository of our published pathways.

The predictions on this page were produced by installing the three packages and calling each pathway's predict on the standardised SMILES. Compound structures were resolved from PubChem and were then processed by the same standardiser as the training tables.

1. Li X. Q., Björkman A., Andersson T. B., Ridderström M., Masimirembwa C. M. (2002). Amodiaquine Clearance and Its Metabolism to N-Desethylamodiaquine Is Mediated by CYP2C8: A New High Affinity and Turnover Enzyme-Specific Probe Substrate. The Journal of Pharmacology and Experimental Therapeutics 300 399-407.doi.org/10.1124/jpet.300.2.399Amodiaquine is bioactivated to a reactive quinone imine that forms protein adducts.
2. Rendic S., Guengerich F. P. (2020). Metabolism and Interactions of Chloroquine and Hydroxychloroquine with Human Cytochrome P450 Enzymes and Drug Transporters. Current Drug Metabolism 21 1127-1135.doi.org/10.2174/1389200221999201208211537Chloroquine undergoes only limited hepatic metabolism.
3. Ducharme J., Farinotti R. (1996). Clinical Pharmacokinetics and Metabolism of Chloroquine. Clinical Pharmacokinetics 31 257-274.doi.org/10.2165/00003088-199631040-00003Clinical pharmacokinetics of chloroquine.
4. Nunes V., Schinoni M. I., Bessone F., Lucena M. I., Medina-Cáliz I., Hernandez N. et al. (2025). Stanozolol-induced Liver Injury: A Distinctive Cholestatic Clinical and Biochemical Phenotype at Presentation. Journal of Clinical and Experimental Hepatology 15 102506.doi.org/10.1016/j.jceh.2025.10250617-alpha-alkylated androgens and the bland cholestasis they cause.
5. Lim L. O., Bortell R., Neims A. H. (1986). Nitrofurantoin inhibition of mouse liver mitochondrial respiration involving NAD-linked substrates. Toxicology and Applied Pharmacology 84 493-499.doi.org/10.1016/0041-008x(86)90254-1Nitrofurantoin impairs mitochondrial respiration in isolated hepatocytes.