Using QML to aid in Antibiotic drug discovery

19 min readMar 28, 2023

Quantum is a word that seems quite daunting, in actuality it is tricky but any reader of this article can learn about it. However, explaining both quantum machine learning and quantum chemistry in one article is quite daunting. Yes, in the future maybe I will be writing explanations of the two, but today I want to focus on using these techniques to create new antibiotics.

So, If you would like to learn a bit about quantum, you can find info here: https://qiskit.org/textbook/ch-states/introduction.html

And without further ado, let’s get into the article!

Analyzing the properties of a molecular structure can be useful for creating new antibiotics, and using quantum computing can be beneficial for this because it allows us to identify novel drug targets and validate them at a very high speed.

But why do we need new antibiotics?

Play along with me… You are experiencing cold and flu symptoms — headache, sore throat, fever, runny nose, etc — so naturally you go to your doctor. He says that you should rest up, drink fluids (regular stuff) and take paracetamol for the headache. You then follow up with “could you prescribe me antibiotics” — this is where the mistake happens

Bacteria have actually evolved into “superbugs”, this basically means that they are resistant to antibiotics — so if your lying in a hospital bed with no medicine to help you with a superbug, blame people who take antibiotics for no real reason.

If nothing is done to halt the rise of superbugs, global deaths from drug-resistant infections could rise from the currently estimated 700,000 every year to 10 million people every year by 2050.

What do we do about it?

So, it is clear that we need to do something about it, but what?

Well, the obvious answer would be to make new drugs — but that is quite a complex task.

Nanotechnology researchers are working on using nanobots to deliver medicine, this is quite a good way to get past the barriers bacterial cells have in place to prevent a medicine from attaching to its ribosomes — hence killing it.

Another way is using AI to aid drug discovery!

This can be done in more things than antibiotics, but AI can be a very useful tool to aid drug discovery.

Predicting drug response: Predicting drug response is an important aspect of drug discovery because it helps researchers identify which drugs are likely to be effective in treating a particular disease, and which patients are likely to respond well to a given drug.
GANs: We can use GANs to simulate molecular structure, Simulating molecular structures helps in drug discovery by allowing researchers to understand how potential drug molecules interact with their target proteins or enzymes at a molecular level. This information can help in designing new drugs with greater specificity, efficacy, and fewer side effects.
QML: QML allows for faster and more accurate predictions of molecular properties and drug interactions. QML can also be used to predict molecular properties such as solubility, stability, and toxicity, which are important factors in drug development. By accurately predicting these properties, researchers can identify potential drug candidates that are safe, effective, and have a higher chance of success in clinical trials.

What are the benefits of QML?

If you are in the AI world, you must be super excited about QML

It can have so many benefits and is probably the most exciting thing in AI currently — if you ask me at least.

Faster computing: Quantum computers can solve certain problems much faster than classical computers, which could have a significant impact on fields such as cryptography, optimization, and drug discovery.
Improved data security: Quantum computers can potentially break certain encryption schemes used to secure sensitive data, but they can also be used to create more secure encryption protocols that are impossible to break using classical computers.
Simulation of complex systems: Quantum computers can simulate the behavior of complex systems that are difficult to model with classical computers, such as large molecules or materials.
Advances in machine learning: Quantum computers can be used to improve machine learning algorithms and enable more accurate and efficient data analysis.
Scientific breakthroughs: Quantum computers could help scientists make significant advances in fields such as chemistry, physics, and materials science by allowing them to simulate and analyze complex systems with greater accuracy and detail.

What are we doing today — an overview

Today we will use quantum computing and machine learning techniques to simulate and analyze the properties of a molecular structure, as said in the title, but what are the steps to doing that?

Quantum circuit design: The first step is to design a quantum circuit that can efficiently encode the molecular structure of a drug. This requires expertise in quantum chemistry and quantum computing. The goal is to encode the relevant information about the molecular structure in a way that can be processed by a quantum computer.
Quantum simulation: Once the quantum circuit is designed, the next step is to simulate the molecular structure of the drug using a quantum computer. This involves using quantum algorithms to calculate the electronic structure of the drug and its interactions with other molecules.
Feature extraction: The output of the quantum simulation is a complex quantum state that contains information about the molecular structure of the drug. The next step is to extract relevant features from this quantum state that can be used as input to a classical machine learning model.
Classical machine learning: The extracted features are then used to train a classical machine learning model to predict the effectiveness of the drug. The model can be trained using a variety of techniques, including supervised learning, unsupervised learning, or reinforcement learning.
Optimization: The final step is to use the machine learning model to optimize the molecular structure of the drug. This involves searching for a molecular structure that maximizes the desired properties of the drug, such as its effectiveness or safety.

So now we know what we are doing — lets do it!

The dataset:

To do this today, I used the PubChem dataset.

So, PubChem is a public database that contains information about millions of chemical compounds. It’s like a giant library, but instead of books, it has information about chemicals!

Scientists and researchers use PubChem to study chemical compounds and their properties. You can search for a specific chemical or browse through the database to find chemicals that match certain criteria.

The dataset contains a lot of information about each chemical compound, such as its molecular formula, weight, structure, and other properties. There are also links to research articles and other resources related to the compound.

The PubChem dataset is a large collection of chemical compounds and their associated properties, such as molecular weight and boiling point. These properties can be used as input features for machine learning models, including those based on quantum machine learning (QML).

QML algorithms are designed to take advantage of the quantum properties of matter to perform certain types of computations more efficiently than classical algorithms. In the case of chemistry, QML can be used to simulate the behavior of molecules and predict their properties. By using the PubChem dataset as training data for QML models, it is possible to develop highly accurate models for predicting the properties of new molecules, which can have applications in drug discovery and materials science.

Step 1: Import our libraries

As usual, the first thing we do is import our libraries, these are gonna be all the libraries needed for all the steps.

As it is QML, we need loads of libraries, So today I will not be explaining each one of them, but I do want to talk briefly about the qiskit library.

Qiskit is a software development kit (SDK) for building quantum computing applications. It’s like a set of tools that allows developers to create and run quantum circuits on real or simulated quantum devices.

Qiskit is a popular open-source platform that’s widely used in the quantum computing community. It has a lot of built-in functionality for things like circuit design, quantum algorithm development, and optimization.

One of the great things about Qiskit is that it’s designed to be accessible to people with different levels of quantum computing knowledge. You can start using Qiskit with just a basic understanding of quantum mechanics, and then gradually build up your skills as you go along.

Qiskit also has a lot of documentation, tutorials, and other resources available online, which makes it easier to learn and get started with. So, if you’re interested in quantum computing, Qiskit is definitely a platform worth exploring.

#import libraries
import numpy as np
from qiskit import Aer, execute
from qiskit.aqua import aqua_globals
from qiskit.aqua.components.optimizers import COBYLA
from qiskit.aqua.components.variational_forms import RYRZ
from qiskit.aqua.operators import Z2Symmetries, WeightedPauliOperator
from qiskit.aqua.algorithms import VQE
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from rdkit import Chem
from rdkit.Chem import AllChem
from qiskit.aqua.utils import split_dataset_to_data_and_labels, map_label_to_class_name
from qiskit.chemistry import FermionicOperator
from qiskit.chemistry.core import Hamiltonian, QubitMappingType
from qiskit.aqua.algorithms import VQE
from qiskit.aqua.components.optimizers import L_BFGS_B
from qiskit.aqua import aqua_globals
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from qiskit.result import Result
from qiskit.quantum_info import partial_trace
from qiskit.circuit.library import PauliFeatureMap
from qiskit.opflow import Z, X, I
from qiskit.aqua.algorithms import VQC
from sklearn.linear_model import LogisticRegression

Step 2: Load and preprocess the data

Now let’s get our dataset in the mix, because without a dataset what even is AI?

# Load the bioassay data from a CSV file
bioassay_df = pd.read_csv('PubChem_BioAssay.csv')

Now we need to make sure we are looking at the part of the data about antibiotics

# Filter the data to include only compounds with relevant bioactivity data
filtered_df = bioassay_df[bioassay_df['Antibiotic Activity'] != '']

Now we standardize the data

# Standardize the molecular structures using RDKit
filtered_df['Molecule'] = filtered_df['SMILES'].apply(lambda x: Chem.MolFromSmiles(x))
filtered_df['Molecule'] = filtered_df['Molecule'].apply(lambda x: Chem.AddHs(x))
filtered_df['Molecule'] = filtered_df['Molecule'].apply(lambda x: Chem.MolToSmiles(x))
filtered_df['Molecule'] = filtered_df['Molecule'].apply(lambda x: Chem.MolFromSmiles(x))

First, we take the ‘SMILES’ column of the dataset, which contains the SMILES notation for each molecule, and use the RDKit function ‘Chem.MolFromSmiles’ to convert each SMILES notation into a molecule object.

SMILES stands for Simplified Molecular Input Line Entry System. It’s basically a way to represent the structure of a molecule using a short string of characters.

For example, instead of drawing out the structure of caffeine on a piece of paper, you could use the SMILES notation to represent it as “CN1C=NC2=C1C(=O)N(C(=O)N2C)C”.

This notation is really useful for computer programs because it’s easy to read and write, and can be used to search large databases of molecules.

Then, we add hydrogen atoms to the molecule objects using ‘Chem.AddHs’. Hydrogen atoms are not usually included in SMILES notation but are important for determining the geometry and properties of the molecule.

The next line of code converts the modified molecule objects back into SMILES notation using ‘Chem.MolToSmiles’, and then again back into molecule objects using ‘Chem.MolFromSmiles’. This step is important for standardizing the molecule structures and ensuring that they are consistent and can be properly processed by the subsequent steps of the analysis.

# Calculate the relevant molecular descriptors using RDKit
filtered_df['Descriptors'] = filtered_df['Molecule'].apply(lambda x: AllChem.GetMorganFingerprintAsBitVect(x, 2))

In this part of the code, we are calculating the relevant molecular descriptors for our filtered dataset using the RDKit library in Python. These descriptors are important because they provide numerical representations of the molecules that we can use to train our machine learning model.

We start by taking the ‘Molecule’ column from our filtered dataset, which contains the molecular structures standardized using the SMILES notation, and apply the RDKit function ‘AllChem.GetMorganFingerprintAsBitVect’ to each molecule. This function calculates the Morgan fingerprint of the molecule, which is a binary vector that encodes the presence or absence of molecular substructures (called “features”) at different radii around each atom in the molecule.

The function takes two arguments: the first argument is the molecule itself, which we pass as an argument to the lambda function. The second argument is the radius of the fingerprint, which we set to 2 in this case. This means that the function will consider all features that are within two bonds of each atom.

Finally, we store the calculated molecular descriptors in a new column called ‘Descriptors’ in our filtered dataset. This column contains a binary vector for each molecule, which we can use as input features for our machine learning model.

# Calculate the relevant molecular descriptors using RDKit
filtered_df['Descriptors'] = filtered_df['Molecule'].apply(lambda x: AllChem.GetMorganFingerprintAsBitVect(x, 2))

Now, We are labeling the compounds in the filtered dataframe based on their antibiotic activity. If the antibiotic activity is ‘Active’, we label it as 1, indicating that it is active against bacteria. Otherwise, we label it as 0, indicating that it is not active. This will help us in training our machine learning model to predict the antibiotic activity of compounds based on their molecular descriptors.

# Label the compounds according to their antibiotic activity
filtered_df['Antibiotic Label'] = filtered_df['Antibiotic Activity'].apply(lambda x: 1 if x == 'Active' else 0)

And we can consider this data preprocessed!

Step 3: Define the quantum circuit to encode the molecular structure

As said before, now we design a quantum circuit that can efficiently encode the molecular structure of a drug. This requires expertise in quantum chemistry and quantum computing. We want to encode the relevant information about the molecular structure in a way that can be processed by a quantum computer.

# Define the molecule to encode
molecule = 'O=C(NC1CC1)c1ccc(Cl)cc1'

This line defines the molecule that we want to encode. In this case, it's a molecule called "O=C(NC1CC1)c1ccc(Cl)cc1".

# Define the feature map to encode the molecule
num_qubits = 4 # You may need to adjust this based on the size of your molecule
paulis = []
for i in range(num_qubits):
    for j in range(i):
        paulis.append(Z ^ (i & j) * X ^ ((i + j) % num_qubits))
feature_map = PauliFeatureMap(feature_dimension=num_qubits, reps=1, paulis=paulis)

This section defines the feature map that we will use to encode the molecule into a quantum circuit. The feature map is essentially a set of instructions that we will use to translate the properties of the molecule into a format that can be represented by a quantum computer.

First, we define the number of qubits that we will use to represent the molecule. In this case, we have set the number of qubits to 4, but this number may need to be adjusted depending on the size of the molecule.

Next, we create an empty list called "Paulis". We will use this list to store the set of instructions that make up our feature map.

The next few lines of code create the set of instructions for our feature map. Specifically, we are creating a set of Pauli gates that we will use to represent the molecule. Pauli gates are a type of quantum gate that are commonly used in quantum computing.

Pauli gates are a set of quantum logic gates in quantum computing that are used to manipulate the state of qubits. There are three types of Pauli gates: the X gate, the Y gate, and the Z gate.

We used the X and Z gates in the context of quantum computing, specifically in the context of encoding a molecule using a quantum circuit. In this context, we used the X and Z gates as part of a feature map to encode the molecular structure of a given molecule onto a quantum circuit.

The X gate is often referred to as the “bit-flip” gate because it maps the computational basis state 0 to 1 and vice versa. It is a one-qubit gate that flips the state of a qubit from |0⟩ to |1⟩ or vice versa. The Z gate, on the other hand, is often referred to as the “phase-flip” gate because it maps the computational basis state |0⟩ to itself, while mapping the state |1⟩ to -|1⟩. It is also a one-qubit gate that applies a phase shift to a qubit.

In the feature map we used, we constructed Pauli strings that consisted of tensor products of X and Z gates acting on different qubits. These Pauli strings acted as a way of encoding the molecular structure of the molecule onto the quantum circuit, with each Pauli string representing a different aspect of the molecule’s structure.

By using X and Z gates in this way, we were able to effectively map the molecular structure of the given molecule onto the quantum circuit, which is important for further analysis and simulations.

# Initialize the circuit to encode the molecule
circuit = feature_map.construct_circuit(molecule)

Finally, we use our feature map to initialize the quantum circuit that will encode the molecule. We do this by calling the "construct_circuit" method of our feature map and passing in the molecule that we want to encode. The resulting circuit is a representation of the molecule that can be used as input to a quantum algorithm.

Step 4: Simulate the molecular structure using a quantum computer

Of course I do not have a quantum computer laying around at home, but this is what I would do if I did have access to one.

So, in this code, we’re basically trying to simulate the behavior of a molecule called ‘O=C(NC1CC1)c1ccc(Cl)cc1’ using a quantum computer.

# Define the molecule
molecule = 'O=C(NC1CC1)c1ccc(Cl)cc1'
basis = 'sto-3g'
charge = 0
spin = 0
driver = PySCFDriver(molecule=molecule, basis=basis, charge=charge, spin=spin)
molecule_data = driver.run()

We're defining a molecule using its chemical formula O=C(NC1CC1)c1ccc(Cl)cc1. We're also setting other parameters such as the basis set to sto-3g, the charge to 0 and the spin to 0. Then we're creating a PySCF driver object to compute molecular data for the given molecule and its parameters. We store this data in a variable called molecule_data.

# Define the Hamiltonian
hamiltonian = Hamiltonian(qubit_mapping=QubitMappingType.PARITY, two_qubit_reduction=True, freeze_core=True, orbital_reduction=[-3, -2])

Here we're defining a Hamiltonian which is an operator that describes the energy of a quantum system. The Hamiltonian takes parameters such as the qubit mapping, which is set to QubitMappingType.PARITY, and the two-qubit reduction flag, which is set to True to allow for more efficient computations. We also set the freeze_core flag to True to remove the core electrons and orbital_reduction to remove some of the virtual orbitals.

In quantum computing, the Hamiltonian is used to represent the problem to be solved as a set of quantum states, which can then be manipulated and measured using quantum algorithms. This allows for the simulation of complex quantum systems, such as chemical reactions or materials properties, which are difficult to study using classical computers.

By using the Hamiltonian to represent a system, we can take advantage of the properties of quantum mechanics, such as superposition and entanglement, to perform computations in a way that is exponentially faster than classical algorithms.

# Map the fermionic operator to the qubit operator
ferm_op = FermionicOperator(h1=molecule_data.one_body_integrals, h2=molecule_data.two_body_integrals)
qubit_op = ferm_op.mapping(map_type=hamiltonian.qubit_mapping, threshold=hamiltonian.mapping_threshold)
qubit_op = qubit_op.two_qubit_reduced_operator(hamiltonian.two_qubit_reduction)

We're mapping the fermionic operator, which describes the behavior of electrons in the molecule, to a qubit operator which is used in quantum computing. We're using the molecular data computed earlier to create the fermionic operator. Then, we're using the Hamiltonian object to map the fermionic operator to a qubit operator. Finally, we're reducing the number of qubits in the qubit operator using the two_qubit_reduced_operator function.

# Define the VQE algorithm
optimizer = L_BFGS_B(maxfun=1000)
vqe = VQE(qubit_op, ansatz=None, optimizer=optimizer)

We're defining a Variational Quantum Eigensolver (VQE) algorithm, which is used to find the minimum energy of a molecule. We're using the qubit operator we mapped earlier and an optimizer called L_BFGS_B. We also set the maximum number of function evaluations to 1000.

# Run the VQE algorithm
backend = Aer.get_backend('statevector_simulator')
result = vqe.run(backend)

Here we're running the VQE algorithm using the statevector_simulator backend provided by IBM's Qiskit. The statevector_simulator backend simulates the behavior of a quantum system using a classical computer. We store the result of the computation in a variable called result.

# Print the results
print(f"Optimal energy: {result.optimal_value}")
print(f"Optimal parameters: {result.optimal_point}")

Finally, we're printing the optimal energy and parameters found by the VQE algorithm

Step 5: define relevant features

The result of a quantum simulation aimed at understanding the molecular structure of a drug is a complex quantum state. This state contains a wealth of information that can provide insights into the drug’s properties and behavior. However, to make use of this information in practical applications, it is necessary to extract relevant features from the quantum state. This involves analyzing and interpreting the data to identify the most significant details, which can then be used as input to a classical machine learning model.

# Get the result of the quantum simulation
result = vqe.run(backend)

In this block, the code is running a quantum simulation using a Variational Quantum Eigensolver (VQE) algorithm on the specified backend. The result of this simulation is stored in the result variable.

# Extract the statevector from the result
statevector = result.get_statevector()

Here, the statevector is extracted from the result of the quantum simulation and stored in the statevector variable. The statevector represents the quantum state of the system at the end of the simulation.

# Compute the density matrix from the statevector
density_matrix = np.outer(statevector, np.conj(statevector))

This code block computes the density matrix from the statevector. The density matrix is a mathematical object that describes the statistical properties of a quantum system. The np.outer() function calculates the outer product of the statevector with its complex conjugate, which results in a matrix representation of the density matrix.

# Trace out the qubits corresponding to the electrons
electrons = [0, 1, 2, 3]  # Example: trace out the first 4 qubits
rdm = partial_trace(density_matrix, electrons)

This block of code performs a partial trace on the density matrix to obtain the reduced density matrix (RDM) for a subset of qubits. The electrons variable specifies the subset of qubits that correspond to the electrons, and the partial_trace() function calculates the trace over these qubits to obtain the RDM.

# Flatten the RDM to a 1D array
features = np.ndarray.flatten(rdm)

Finally, the RDM is flattened into a 1D array using the np.ndarray.flatten() function and stored in the features variable. This is done to obtain a feature vector that can be used for further analysis or processing.

Step 6: Classic machine learning model

After extracting features from the data, a classical machine learning model can be trained to predict the efficacy of the drug using different techniques such as supervised learning, unsupervised learning, or reinforcement learning.

# Split the dataset into training and test sets
X_train, X_test, y_train, y_test = train_test_split(features, labels, test_size=0.2, random_state=42)

This line of code uses the train_test_split function from the scikit-learn library to randomly split the dataset into two parts: X_train and y_train, which will be used to train the model, and X_test and y_test, which will be used to evaluate the model's performance. features and labels are the input data and output data, respectively. test_size=0.2 specifies that 20% of the data should be used for testing, while the remaining 80% will be used for training. The random_state=42 parameter ensures that the random split is reproducible.

# Train a logistic regression model
model = LogisticRegression(max_iter=1000)
model.fit(X_train, y_train)

The LogisticRegression function from scikit-learn is used to create a logistic regression model. max_iter=1000 specifies the maximum number of iterations the solver should run for convergence. model.fit(X_train, y_train) fits the model to the training data X_train and y_train.

Logistic regression is a good choice for binary classification problems where the goal is to predict one of two possible outcomes. In this case, the aim is to predict whether a molecule is stable or not based on its features, which is a binary classification problem. Logistic regression models the probability of a binary response variable as a function of the predictor variables, which in this case are the features extracted from the reduced density matrix.

Logistic regression is also relatively simple and fast to train compared to other machine learning models, making it a popular choice for binary classification problems with small to medium-sized datasets. It is also less prone to overfitting than more complex models, which can be an advantage when working with limited amounts of data. Finally, the output of logistic regression is easily interpretable, as it provides a probabilistic estimate of the likelihood of a molecule being stable based on its features.

# Evaluate the model on the test set
accuracy = model.score(X_test, y_test)
print("Accuracy:", accuracy)

The score method of the logistic regression model is used to evaluate the accuracy of the model on the test set. The resulting accuracy score is stored in the accuracy variable. Finally, the accuracy score is printed to the console using print("Accuracy:", accuracy).

Step 7: Optimise

Finally, we use the machine learning model to optimize the molecular structure of the drug.

# Define the molecule to optimize
molecule = 'O=C(NC1CC1)c1ccc(Cl)cc1'

This line of code sets the SMILES string of the molecule to be optimized.

# Set up the PySCF driver to compute the molecular integrals
driver = PySCFDriver(atom=molecule, unit=UnitsType.ANGSTROM, charge=0, spin=1)

This uses the PySCFDriver class from the Qiskit Chemistry module to set up the driver to compute the molecular integrals. atom=molecule specifies the molecular structure to compute the integrals for, unit=UnitsType.ANGSTROM sets the unit of length to angstroms, and charge=0 and spin=1 specify the charge and spin of the molecule.

# Use the driver to compute the one- and two-electron integrals
molecule_data = driver.run()
one_body_integrals = molecule_data.one_body_integrals
two_body_integrals = molecule_data.two_body_integrals

Here we use the run method of the driver object to compute the one- and two-electron integrals of the molecule. The results are stored in molecule_data. one_body_integrals and two_body_integrals are then extracted from molecule_data for further processing.

# Convert the integrals to a fermionic operator
fermionic_op = FermionicOperator(h1=one_body_integrals, h2=two_body_integrals)

This code creates a FermionicOperator object from the one- and two-electron integrals. A FermionicOperator is a class that represents a molecular Hamiltonian in the second quantization formalism.

# Map the fermionic operator to a qubit operator using the parity mapping
num_particles = molecule_data.num_alpha + molecule_data.num_beta
num_spin_orbitals = molecule_data.num_orbitals * 2
qubit_op = fermionic_op.mapping(map_type='parity', threshold=0.00000001).two_qubit_reduced_operator(num_particles)

This maps the fermionic_op object to a qubit_op object using the parity mapping. The map_type='parity' parameter specifies that the parity mapping should be used, and threshold=0.00000001 specifies the threshold for eliminating small terms in the mapping. The two_qubit_reduced_operator method reduces the number of qubits required to represent the operator. The num_particles and num_spin_orbitals variables are used to determine the number of particles and spin orbitals in the system.

# Convert the qubit operator to a weighted Pauli operator
qubit_op = WeightedPauliOperator(qubit_op)

This code converts the qubit_op object to a WeightedPauliOperator object. A WeightedPauliOperator is a class that represents a weighted sum of Pauli operators.

# Define the optimizer and variational form for the VQE algorithm
optimizer = COBYLA(maxiter=1000)
var_form = RYRZ(qubit_op.num_qubits, depth=3)

Finally, this defines the optimizer and variational form for the Variational Quantum Eigensolver (VQE) algorithm. COBYLA is the optimizer used, and

And we are done!!!

Conclusion

Quantum Machine Learning (QML) can be a powerful tool in drug discovery due to its ability to potentially handle and analyze vast amounts of complex data with high accuracy and efficiency. QML algorithms have the potential to analyze large datasets in parallel, which can be useful for tasks such as identifying correlations between molecular structure and biological activity.

Key takeaways: