Acetylecholinesterase (AChE) is an important enzyme that breaks down acytelcholine (a neurotransmitter). It is commonly found throughout the body but most commonly at neuromuscular junctions. AChE is a target of various drugs; inhibitors of AChE are used to treat disease such as myasthenia gravis and Alzheimer’s disease.
In ChEMBL, we identify 8,832 molecules/drugs with IC50 data targeting AChE. Some of these drugs are more potent than others; we consider IC50 standard values less than 1,000 as active, more than 10,000 as inactive, and in between as intermediate. In drug discovery, we want to know how features of these molecules interact to determine the potency. Framed as a causal inquiry, if we knew the causal relationships of the molecular features, we would be better informed on how and where to manipulate the features to improve drug efficacy (potency).
The features we use come from Lipinski’s Rule of Five with 881 PaDEL fingerprints. The Lipinski features are continuous while the PaDEL fingerprints are binary. Not all PaDEL fingerprints are used since some of them have low variance (ultimately, 144 of 881 PaDEL fingerprints were used). The target variable (IC50) is continuous but can also be discretized as described above (into active, intermediate and inactive).
Three Phase Dependency Analysis (TPDA) is the structure learning algorithm applied to final data consiting of 149 features (144 PaDEL fingerprints, 4 Lipinski descriptors, 1 IC50 feature) across 6,157 molecules (records with the same SMILES are removed keeping only 1). In particular, we used pybbn-spark to learn the full causal model. The tree structure is shown below with the
- IC50 node as blue,
- molecular weight as orange,
- octonal-water partition coefficient (solubility) as green,
- number of hydrogen bond donors as purple,
- number of hydrogen bond acceptors as gray, and
- all PaDEL nodes as red.

As can be seen, with the exception of the number of hydrogen bond donors node, the Lipinksi and IC50 features are very close to one another graphically. The full directed acyclic graph (DAG) is shown below. It turns out that the octonal-water partition coefficient and number of hydrogen bond donors nodes are connected.

Let’s direct our attention to the Markov blankets of the Lipinski features and IC50. Notice that the molecular weight (mol_wt, orange) is a confounder of the octonal-water partition coefficient( mol_log_p, green) and number of hydrogen acceptors (num_h_acceptors, gray).

To estimate the true causal effect of “num_h_acceptors” (number of hydrogen acceptors) on “mol_log_p” (octonal-water partition coefficient), we can apply the do-operator.
Denote “mol_log_p” as and “num_h_acceptor” as
, then we estimate the causal effects of
on
. For this estimation, we use py-bbn.
| do-operation | Causal Effect |
| 13% | |
| 20% | |
| 21% |
For drug discovery and potency, this model suggest that the molecular solubility has a direct causal relationship with IC50. The probability of the IC50 being potent (or active) is highest when the molecular solubility is high. One way to get molecular solubility to be high is by decreasing the number of hydrogen acceptors (demonstrated by the do-operator above). Further analysis can be conducted to determine which other features can improve IC50 by looking at other features in the Markov blanket of the Lipinski nodes.
Contact us if you would like to learn more about using causal learning and inference for drug discovery. We would love to hear from you and learn about the problems you are trying to solve.

Leave a Reply