Chemical Informatics Fundamentals

The evolution of artificial intelligence in chemistry involves several stages. Early efforts focused on digitizing chemical knowledge into databases using various representation methods. Mid-stage approaches relied on manual feature engineering to encode existing data, followed by traditional machine learning algorithms for prediction tasks. The current era is defined by the widespread adoption of deep learning networks, which has spurred the development of novel molecular fingerprinting techniques.

Molecular Representation: SMILES

SMILES (Simplified Molecular Input Line Entry System) is a standard notation used in chemoinformatics to represent chemical structures as strings of ASCII characters. Proposed by Weininger et al., this system encodes atoms, bonds, charges, and structural features like rings or side chains according to specific syntax rules:

• Atoms are represented by their chemical symbols.
• Bonds are denoted by specific characters (= for double, # for triple).
• Special groups or ions use square brackets, e.g., [Cu+2] represents a copper ion with a +2 charge.

Molecular Fingerprints

A molecular fingerprint is a fixed-length bit vector that acts as a unique identifier for a molecule's structural characteristics. Each bit posiiton corresponds to a specific chemical substructure. For instance, a 2-bit vector might indicate the presence of a methyl group at index 0 and a benzene ring at index 1. A vector [0, 1] implies the molecule contains a benzene ring but lacks a methyl group.

Tooling: RDKit

RDKit is an open-source cheminformatics library supporting Windows, macOS, and Linux. It facilitates integration with Python, Java, and C, serving as a primary tool for processing molecular data in computational chemistry workflows.

Modeling Methodologies

Machine learning tasks generally fall into classification (discrete outputs) or regression (continuous outputs). In earlier experiments involving reaction data, Random Forest models were employed, combining multiple decision trees trained on random subsets of data.

Deep Learning, a subset of machine learning, automates feature extraction through neural networks, eliminating the need for laborious hand-crafted features. Since SMILES sequences resemble language, Natural Language Processing (NLP) techniques are applicable here.

Recurrent Neural Networks (RNN)

RNNs are designed for sequential data. They process inputs step-by-step, maintaining a hidden state that carries information from previous steps. The mathematical update rule typically follows:

h_n = σ(W_hh * h_{n-1} + W_hx * x_n + b)
y_n = Softmax(V * h_n + c)

The hidden state allows the network to "read" past sequence context to inform future predictions. However, standard RNNs face limitations regarding long-term dependency due to vanishing gradients and poor parallelization capabilities during training, as computations must occur sequentially.

Implementation Details

Imports and Dependencies

import re
import time
import pandas as pd
from typing import List, Tuple
import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import Dataset, DataLoader, Subset

Model Architecture

class SequenceYieldNetwork(nn.Module):
    def __init__(self, vocab_size, feature_dim, hidden_units, out_units, layers, dropout_rate, device):
        super(SequenceYieldNetwork, self).__init__()
        self.embedding_layer = nn.Embedding(vocab_size, feature_dim)
        self.recurrent_cell = nn.RNN(
            input_size=feature_dim,
            hidden_size=hidden_units,
            num_layers=layers,
            batch_first=True,
            dropout=dropout_rate,
            bidirectional=True
        )
        self.output_projection = nn.Sequential(
            nn.Linear(2 * layers * hidden_units, out_units),
            nn.Sigmoid(),
            nn.Linear(out_units, 1),
            nn.Sigmoid()
        )

    def forward(self, x):
        embedded = self.embedding_layer(x)
        _, final_states = self.recurrent_cell(embedded)
        states_transposed = final_states.transpose(0, 1)
        flattened = states_transposed.reshape(states_transposed.size(0), -1)
        prediction = self.output_projection(flattened).squeeze(-1)
        return prediction

Tokenization Logic

class MoleculeSeqTokenizer():
    def __init__(self, pad_symbol, pattern, vocab_path, max_seq_len):
        self.pad_token = pad_symbol
        self.regex_pattern = pattern
        self.vocab_path = vocab_path
        self.max_length = max_seq_len

        with open(self.vocab_path, 'r') as f:
            lines = [line.strip('\n') for line in f.readlines()]
        
        self.token_map = {token: idx for idx, token in enumerate(lines)}

    def _apply_regex(self, smiles_list):
        compiled_re = re.compile(f"({self.regex_pattern}|.)")
        processed_tokens = []
        for smi in smiles_list:
            match = compiled_re.findall(smi)
            if len(match) > self.max_length:
                match = match[:self.max_length]
            processed_tokens.append(match)
        return processed_tokens

    def tokenize(self, smiles_list):
        raw_tokens = self._apply_regex(smiles_list)
        tagged_tokens = [["<CLS>"] + t + ["<SEP>"] for t in raw_tokens]
        padded_tokens = self._pad_sequences(tagged_tokens, self.pad_token)
        indexed_tokens = self._map_to_indices(padded_tokens)
        return indexed_tokens

    def _pad_sequences(self, seqs, pad_val):
        max_len = max(len(seq) for seq in seqs)
        return [seq + [pad_val] * (max_len - len(seq)) for seq in seqs]

    def _map_to_indices(self, tokens):
        indices_list = []
        for token_group in tokens:
            current_idx = []
            for token in token_group:
                if token in self.token_map:
                    current_idx.append(self.token_map[token])
                else:
                    new_id = max(self.token_map.values()) + 1
                    self.token_map[token] = new_id
                    current_idx.append(new_id)
            indices_list.append(current_idx)
        return indices_list

Data Loading Pipeline

def load_reaction_dataset(file_path, is_training=True):
    df = pd.read_csv(file_path)
    
    r1 = df["Reactant1"].tolist()
    r2 = df["Reactant2"].tolist()
    prod = df["Product"].tolist()
    addi = df["Additive"].tolist()
    solv = df["Solvent"].tolist()
    
    yields = df["Yield"].tolist() if is_training else [0] * len(r1)
    
    combined_input = []
    for react_a, react_b, product, additive, solvent in zip(r1, r2, prod, addi, solv):
        connection_string = ".".join([react_a, react_b])
        full_string = ">".join([connection_string, product])
        combined_input.append(full_string)
    
    paired_data = list(zip(combined_input, yields))
    return paired_data

class ReactionDataSubset(Dataset):
    def __init__(self, dataset_items: List[Tuple[List[str], float]]):
        self.records = dataset_items

    def __len__(self):
        return len(self.records)

    def __getitem__(self, idx):
        return self.records[idx]

def collate_fn(batch):
    regex_pattern = r"[[^\\]]+]|Br?|Cl?|N|O|S|P|F|I|b|c|n|o|s|p|\\(|\\)|\\.|=|#|-|\\+|\\\\\\|/|:|~|@|\\?|>|\\*|\\$|[0-9]{2}|[0-9]"
    tokenizer = MoleculeSeqTokenizer("<PAD>", regex_pattern, "../vocab_full.txt", max_length=300)
    
    smi_batch = []
    yield_batch = []
    
    for item in batch:
        smi_batch.append(item[0])
        yield_batch.append(item[1])
    
    token_tensors = torch.tensor(tokenizer.tokenize(smi_batch)[1])
    target_tensor = torch.tensor(yield_batch)
    
    return token_tensors, target_tensor

Training Procedure

def execute_training_loop():
    subset_size = 10
    vocab_dimension = 294
    input_dim = 300
    hidden_dim = 512
    output_dim = 512
    layers = 10
    drop_prob = 0.2
    clip_value = 1
    epochs = 100
    lr = 0.0001
    
    start_time = time.time()
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    
    train_data = load_reaction_dataset("../dataset/round1_train_data.csv")
    main_dataset = ReactionDataSubset(train_data)
    train_loader = DataLoader(main_dataset, batch_size=128, shuffle=True, collate_fn=collate_fn)
    
    model = SequenceYieldNetwork(vocab_dimension, input_dim, hidden_dim, output_dim, layers, drop_prob, device).to(device)
    model.train()
    
    optimizer = optim.Adam(model.parameters(), lr=lr)
    loss_fn = nn.L1Loss()
    best_error = 10
    
    for epoch in range(epochs):
        total_loss = 0
        for src, targets in train_loader:
            src, targets = src.to(device), targets.to(device)
            optimizer.zero_grad()
            preds = model(src)
            loss = loss_fn(preds, targets)
            loss.backward()
            torch.nn.utils.clip_grad_norm_(model.parameters(), clip_value)
            optimizer.step()
            total_loss += loss.item()
        
        avg_loss = total_loss / len(train_loader)
        print(f'Epoch: {epoch+1:02} | Train Loss: {avg_loss:.3f}')
        
        if avg_loss < best_error:
            torch.save(model.state_dict(), '../model/RNN.pth')
    
    elapsed_mins = (time.time() - start_time) / 60
    print(f"Total execution time: {elapsed_mins:.2f} minutes")

if __name__ == '__main__':
    execute_training_loop()

Inference and Submission

import torch
from torch.utils.data import DataLoader
from dataset import ReactionDataSubset, collate_fn
from model import SequenceYieldNetwork
from data_utils import load_reaction_dataset

def generate_submission_results(model_checkpoint, output_dest):
    params = {
        'embed_size': 294, 'input_size': 300, 'hidden_size': 512,
        'output_size': 512, 'layers': 10, 'dropout': 0.2
    }
    device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
    
    test_records = load_reaction_dataset("../dataset/round1_test_data.csv", train=False)
    test_dataset = ReactionDataSubset(test_records)
    test_loader = DataLoader(test_dataset, batch_size=64, shuffle=False, collate_fn=collate_fn)

    network = SequenceYieldNetwork(**params).to(device)
    network.load_state_dict(torch.load(model_checkpoint))
    network.eval()
    
    results = []
    with torch.no_grad():
        for src, _ in test_loader:
            src = src.to(device)
            pred = network(src)
            results.extend(pred.detach().tolist())
    
    submission_lines = ['rxnid,Yield']
    for i, val in enumerate(results):
        submission_lines.append(f'test{i+1},{val:.4f}')
    
    with open(output_dest, 'w') as f:
        f.write('\n'.join(submission_lines))

    print("Submission generation complete.")

generate_submission_results("../model/RNN.pth", "../output/RNN_submit.txt")

Performance Evaluation

Empirical testing was conducted on local server hardware. The initial run required approximately two hours of processing time. Comparative analysis against baseline tree-based ensemble models indicated that the RNN architecture achieved lower predictive accuracy in this specific context. This suggests that while RNNs are theoretically capable of handling sequential molecular representations, they may not be the optimal choice for this particular task without further architectural adjustments or hyperparameter optimization.