Tutorial 2: Human Health Impact Prediction Using ANN Model

# Cell 1: Import Required Libraries / Packages

import numpy as np                    # Stores numbers in arrays and runs fast calculations on them
import pandas as pd                   # Opens the data file as a table and finds rows by their SMILES string
import matplotlib.pyplot as plt       # Plot Graphs
import ipywidgets as widgets          # Creates UI elements such as the slider below
from IPython.display import display   # Displays UI elements

import tensorflow as tf               # Loads the trained neural network (ANN) and runs it to make predictions

import os                             # Allows interactions with the file system
import requests                       # Makes requests to internet files

import joblib                         # Parallel Computing

# Cell 2: Define Function to Download Weight (.h5 file) of ANN Model

def download_if_needed(url, local_path):
    """Download a file from GitHub if it doesn't already exist."""
    if os.path.exists(local_path):
        print(f"✓ Using existing {local_path}")
        return

    print(f"Downloading {local_path}...")
    response = requests.get(url)
    response.raise_for_status()

    with open(local_path, "wb") as f:
        f.write(response.content)

    print("Download complete.")

# Cell 3: Create Path to Dowloaded Files and Database

# URL / path to download (.h5 file)
HH_MODEL_URL = ("https://raw.githubusercontent.com/cacherowan/CACHE-Rowan/main/Reference_Files/Chemical_Property_Database/model_v5_HH_okayvaltest.h5")

# Name of weight file
HH_MODEL_PATH = "model_v5_HH_okayvaltest.h5"

# URL / path to download (.joblib file)
MIN_MAX_SCALAR_URL = ("https://raw.githubusercontent.com/cacherowan/CACHE-Rowan/main/Reference_Files/Chemical_Property_Database/minmax_scaler_0_5-1.joblib")

# Name of min_max_scalar file
MIN_MAX_SCALAR_PATH = "minmax_scaler_0_5-1.joblib"

# URL / path to database (the spreadsheet holding the molecular properties)
DATABASE_PATH = "https://raw.githubusercontent.com/cacherowan/CACHE-Rowan/main/Reference_Files/Chemical_Property_Database/Processed_Solvent_DF_v6_TEST.xlsx"

# Call to Function to Download Weight (.h5 file)
download_if_needed(HH_MODEL_URL, HH_MODEL_PATH)

# Call to Function to Download Min Max Scalar
download_if_needed(MIN_MAX_SCALAR_URL, MIN_MAX_SCALAR_PATH)

# The trained model file (.h5) that predicts Human Health Impact
HH_MODEL = HH_MODEL_PATH

# Min Max Scalar
Min_Max_Scalar = MIN_MAX_SCALAR_PATH

# Cell 4: Load the trained model
# Load the .h5 file, which contains the trained network: its structure and its learned weights.  
# This rebuilds the model as "model_CC," ready to make predictions.  
# compile=false in the code below helps skip the training setup, since we are only using the model to predict (not train it).  

import warnings
warnings.filterwarnings("ignore")

model_HH = tf.keras.models.load_model(
    HH_MODEL_PATH,
    custom_objects={"LeakyReLU": tf.keras.layers.LeakyReLU},
    compile=False
)

# Cell 5: Load the molecule database
# Opens the spreadshet into a table, and uses each molecule's SMILES code as its row name so molecules are easy to look up

db = pd.read_excel(DATABASE_PATH)
db = db.set_index('SMILES')

# Cell 6: Select the thermodynamic properties and the molecular descriptors used as inputs for the Human Health model
# The Human Health model was trained on a specific set of 10 features.
# We must feed the model the exact same features, in the exact same order, every time we ask it to make a prediction.

# Thermodynamic properties:
thermo_feat_HH = [
    'Heat of Vaporization(J/mol)',    # energy needed to turn the liquid into vapour
    'Heat Capacity (kJ/kgC)',         # heat needed to raise the molecule's temperature
    'XLogP',                          # a measure of how "fat-loving" vs "water-loving" it is
    'Pitzer’s Acentric Factor [-]',   # a measure of how non-spherical the molecule is
    'Critical Temperature [K]'        # temperature above which it can't be liquefied
]
# Molecular descriptors:
mol_desc_feat_HH = [
    'Chi0n',                # a connectivity index describing molecular structure
    'HallKierAlpha',        # a shape-related descriptor
    'SMR_VSA7',             # surface-area descriptor related to molar refractivity
    'VSA_EState6',          # surface-area descriptor related to electronic state
    'NumValenceElectrons'   # total number of valence electrons in the molecule
]

# Cell 7: Load the scaler
scaler = joblib.load(Min_Max_Scalar)

feature_names = list(scaler.feature_names_in_)
data_min = scaler.data_min_
data_max = scaler.data_max_
lo, hi = scaler.feature_range

def scale_subset(values):
    """Apply the fitted MinMax scaler to a subset of columns, by name."""
    scaled_values = []

    for feature_name in values.index:
        if feature_name not in feature_names:
            print(f"Warning: '{feature_name}' was not seen by the scaler at fit time.")
            continue

        # Find where this feature is in the scaler's arrays
        position = feature_names.index(feature_name)

        x = values[feature_name]
        x_min = data_min[position]
        x_max = data_max[position]

        # MinMax formula: rescale x from [x_min, x_max] to [lo, hi]
        x_scaled = (x - x_min) / (x_max - x_min) * (hi - lo) + lo

        scaled_values.append(x_scaled)

    return pd.Series(scaled_values, index=values.index)

# Cell 8: Predict Human Health Impact for one molecule
# This function takes a molecule, looks up its properties in the database, feeds them to the trained model, and returns the predicted Human Health Impact.
# Run this cell once.  It won't show any output on its own, it just sets yp the function so the cells below can use it.

def predict_human_health_impact(molecule_name, smiles):
    """
    Predict the Human Health Impact of one molecule.

    molecule_name : name for printing only (e.g. "Methanol") — doesn't affect the result.
    smiles        : SMILES string; must exactly match a row in the database.
    Returns the predicted Human Health Impact in DALY per kg of chemical.
    """
    descriptors = db.loc[smiles]
    descriptors_hh = descriptors.loc[thermo_feat_HH + mol_desc_feat_HH]
    descriptors_hh_scaled = scale_subset(descriptors_hh)

    print(f"\n{molecule_name}  ({smiles})")
    print("-" * 40)
    for feature_name, value in descriptors_hh.items():
        print(f"  {feature_name:<32} {value:>8.4f}")
    print("-" * 40)

    model_input = np.array([list(descriptors_hh)])
    prediction = model_HH.predict(model_input, verbose=0)
    impact_value = prediction.item()

    print(f"  Human Health Impact: {round(impact_value, 4)} DALY/kg\n")
    return impact_value

# Cell 9: Calculate the Human Health Impact of a chemical
# Give the function below two things: the chemical's name (for display) and its SMILES code, which the model uses to look it up.  Run the cell to see its features and its predicted human health impact impact.
# To calculate the impact for another chemical, copy the line below into a new cell and change the name and SMILES.

predict_human_health_impact("Methanol", "CO");

# Cell 10: Example with Ethanol

predict_human_health_impact("Ethanol", "CCO");

# Cell 11: Example with Benzene

predict_human_health_impact("Benzene", "c1ccccc1");

# Cell 12: Example with Toluene

predict_human_health_impact("Toluene", "Cc1ccccc1");

# Cell 13: Example using Phenol (Which will be displayed later on a slider and graph)

# 1. Example: Phenol
MOLECULE = "PHENOL"
SMILES = "Oc1ccccc1"

impact_value_slider = predict_human_health_impact(MOLECULE, SMILES)

# Cell 14: Phenol Continued, Creating Slider to Visualize Impact

kg_slider = widgets.FloatSlider(value=(round(impact_value_slider, 4)), min=(round(impact_value_slider, 4)) * 0.5, max=(round(impact_value_slider, 4)) * 1.5, step=0.1, description='kg amount:')
output_slider = widgets.Output()

def update_slider(change):
    with output_slider:
        output_slider.clear_output()
        kg_slider_kg = kg_slider.value
        total_hh_slider = (round(impact_value_slider, 4)) * kg_slider_kg
        print(f"Health impact:  {total_hh_slider:.6f} Disability-Adjusted Life Years (DALY)")

kg_slider.observe(update_slider, names='value')
display(kg_slider, output_slider)
update_slider(None)

# Cell 15: Phenol Continued, Creating Slider and Graph to Visualize Impact

kg_slider_graph = widgets.FloatSlider(value=(round(impact_value_slider, 4)), min=(round(impact_value_slider, 4)) * 0.5, max=(round(impact_value_slider, 4)) * 1.5, step=0.1, description='kg amount:')
output_graph = widgets.Output()

hh_per_kg = (round(impact_value_slider, 4))
kg_range = np.linspace((round(impact_value_slider, 4)) * 0.5, (round(impact_value_slider, 4)) * 1.5, 20)

def update_graph(change):
    with output_graph:
        output_graph.clear_output(wait=True)
        kg_graph_kg = kg_slider_graph.value
        total_hh_graph = hh_per_kg * kg_graph_kg

        fig, ax = plt.subplots(figsize=(5, 3.5))

        ax.plot(kg_range, hh_per_kg * kg_range, color='#e0685c')
        ax.scatter([kg_graph_kg], [total_hh_graph], color='#e0685c', zorder=5, s=80)
        ax.set_title(f"Health: {total_hh_graph:.6f} Disability-Adjusted Life Years (DALY)")
        ax.set_xlabel(f"kg {MOLECULE}")

        plt.tight_layout()
        plt.show()

kg_slider_graph.observe(update_graph, names='value')
display(kg_slider_graph, output_graph)
update_graph(None)