Title: ETH_Sem6 / DC / W12_OptimizationProblems

Q: OFAT optimization
A: One Factor At a Time:
Bad because it ignores potential interaction between features, but is super easy to use

Q: DoE optimization
A: Design of Experiments:
Identify significant factors by smart screening and optimize them. Then test the robustness upon small changes.
Example for radiopharmaceuticals:

Q: Gaussian (normal) distribution
A: Defined by mean \(\mu\) and variance \(\sigma^2\) with standard deviation \(\sigma\):
\[\mathcal N(x|\mu, \sigma) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}\]

Q: Bayes theorem and example application in GPR
A: \[p(\vec w|\mathcal D) = \frac{p(\mathcal D | \vec w) p(\vec w)}{p(\mathcal D)}\]
For example for fitting a linear regression model \(y = f(x) = w_0 + w_1 x\)

Q: Bayesian optimization in a nutshell
A: 1. Train Gaussian Process Regression on currently known data samples
2. Use posterior predictive distribution to select new samples
3. Measure these samples, go back to 1. if not satisfied

Q: Aquisition functions: Difference Exploitation and Exploration
A: Exploitation: Sample to improve on currently optimal values
Exploration: Sample to discover new minima

Q: Chemistry-oriented packages for Bayesian optimization
A: GAUCHE: https://leojklarner.github.io/gauche/index.html
BoFire: https://github.com/experimental-design/bofire
BayBE: https://emdgroup.github.io/baybe/
ATLAS: https://github.com/aspuru-guzik-group/atlas

Q: How to generate SELFIES in python
A: from rdkit.Chem import MolFromSmiles
from selfies import encoder, decoder
smile = 'C=O'
mol = MolFromSmiles(smile)
selfie = encoder(mol)

Q: Python function to mutate a SELFIE
A: def mutate_selfie(selfie, max_molecules_len, write_fail_cases=False):
    '''Return a mutated selfie string (only one mutation on selfie is performed)
    Mutations are done until a valid molecule is obtained
    Rules of mutation: With a 33.3% probability, either:
        1. Add a random SELFIE character in the string
        2. Replace a random SELFIE character with another
        3. Delete a random character
    Parameters:
    selfie            (string)  : SELFIE string to be mutated
    max_molecules_len (int)     : Mutations of SELFIE string are allowed up to this length
    write_fail_cases  (bool)    : If true, failed mutations are recorded in "selfie_failure_cases.txt"
    Returns:
    selfie_mutated    (string)  : Mutated SELFIE string
    smiles_canon      (string)  : canonical smile of mutated SELFIE string
    '''
    valid=False
    fail_counter = 0
    chars_selfie = get_selfie_chars(selfie)
    while not valid:
        fail_counter += 1
        alphabet = list(selfies.get_semantic_robust_alphabet()) # 34 SELFIE characters
        choice_ls = [1, 2, 3] # 1=Insert; 2=Replace; 3=Delete
        random_choice = np.random.choice(choice_ls, 1)[0]
        # Insert a character in a Random Location
        if random_choice == 1:
            random_index = np.random.randint(len(chars_selfie)+1)
            random_character = np.random.choice(alphabet, size=1)[0]
            selfie_mutated_chars = chars_selfie[:random_index] + [random_character] + chars_selfie[random_index:]
        # Replace a random character
        elif random_choice == 2:
            random_index = np.random.randint(len(chars_selfie))
            random_character = np.random.choice(alphabet, size=1)[0]
            if random_index == 0:
                selfie_mutated_chars = [random_character] + chars_selfie[random_index+1:]
            else:
                selfie_mutated_chars = chars_selfie[:random_index] + [random_character] + chars_selfie[random_index+1:]
        # Delete a random character
        elif random_choice == 3:
            random_index = np.random.randint(len(chars_selfie))
            if random_index == 0:
                selfie_mutated_chars = chars_selfie[random_index+1:]
            else:
                selfie_mutated_chars = chars_selfie[:random_index] + chars_selfie[random_index+1:]
        else:
            raise Exception('Invalid Operation trying to be performed')
        selfie_mutated = "".join(x for x in selfie_mutated_chars)
        sf = "".join(x for x in chars_selfie)
        try:
            smiles = decoder(selfie_mutated)
            mol, smiles_canon, done = sanitize_smiles(smiles)
            if len(selfie_mutated_chars) > max_molecules_len or smiles_canon=="":
                done = False
            if done:
                valid = True
            else:
                valid = False
        except:
            valid=False
            if fail_counter > 1 and write_fail_cases == True:
                f = open("selfie_failure_cases.txt", "a+")
                f.write('Tried to mutate SELFIE: '+str(sf)+' To Obtain: '+str(selfie_mutated) + '\n')
                f.close()
    return (selfie_mutated, smiles_canon)

Q: SAScore
A: Estimates the synthesizability: Score between 1 (easily synthesizable) and 10 (hard to synthesize). The score is calculated as a sum of fragment scores and complexity penalty. The fragment score is based on statistics of the frequency of Extended Connectivity Fingerprints of diameter 4 (ECFP4) fragments from Pipeline Pilot on almost one million molecules obtained from the PubChem database. In python:
import sascorer
from rdkit.Chem import MolFromSmiles
smiles = 'C=O'
mol = MolFromSmiles(smiles)
score = sascorer.calculateScore(mol)
Sources:
Review: https://doi.org/10.1186/s13321-023-00678-z
ECFP4: https://doi.org/10.1021/ci100050t
Pipeline Pilot: https://link.springer.com/doi/10.1007/s11030-006-9041-5

Q: Example python function to identify forbidden substructures (that are theoretically possible but chemically unfeasible)
A: def forbidden_structure(mol):
    """
    Check for forbidden fragments
    Return True: contains a substructure violation
    Return False: No substructure violation
    """
    violation = False
    forbidden_fragments = [
        '[A;R]=[*;R2]', # Bredt's rule
        #'a~[*;R2]~a', # Atom in aromatic ring connected to two rings
        '*1=**=*1', # Cyclobutadiene of all elements
        '*1***1', # Cyclobutene of all elements
        '*=*1*=***=*1', # Unstable atom double bond arrangment
        '[#6-]', # Negatively charged carbon
        '[#16-]', # Negatively charged sulphur
        '[#8-]', # Negatively charged oxygen
        '[#7-]', # Negatively charged nitrogen
        '[*+]', # Positively charged any atom, this makes the ones above redundant
        '[*-]', # Negatively charged any atom, this makes the ones above redundant
        '[#15;H]', # Phosphorous with hydrogen atom
        '[#15]', # Phosphorous atoms
        '[#7&X5]', # Pentavalent nitrogen
        '*=[#16;!R]', # Double bond to non-ring S
        '[#16&X3]', # Trivalent sulphur
        '[#16&X4]', # Tetravalent sulphur
        '[#16&X5]', # Pentavalent sulphur
        '[#16&X6]', # Hexavalent sulphur
        '[#5,#7,#8,#16]~[F,Cl,Br,I]', # Halogens connected to B, N, O or S
        '*=*=*', # Allene
        #'*#*', # Triple bond
        '[#8,#16]~[#8,#16]', # Bonds between O-O or S-S
        '[#7,#8,#16]~[#7,#8,#16]~[#7,#8,#16]', # N-O-S bonds including permutations
        '[#7,#8,#16]~[#7,#8,#16]~[#6]=,:[#7,#8,#16;!R]', # N-N-C=N bonds where last N is not in ring (N can be O or S as well)
        '*=N-[*;!R]', # *=N-* bond where last * is not in a ring
        '*~[#7,#8,#16]-[#7,#8,#16;!R]', # *-N-N bonds where last N is not in a ring
        'c#ccccc', # benzene triple bond (originally separate "bad" filter)
        "[#5]", # Boron
        "[S;X2]", # non-aromatic sulfide
        "[c;$(c1nc(Br)cc([OH])c1),$(c1nc(Br)ccc([OH])1),$(c1nc(Br)c([OH])cc1),$(c1([OH])nc(Br)ccc1)]", # OH on pyridine ring
        "[N;H1]([OH])", # hydroxylamine group
        "[OX2H][#6X3]=[#6]", # enol
        "[C;X3]=[N;X2]", # imine group
        #"[CX3;$([C]([#6])[#6]),$([CH][#6])]=[NX2][#6]", # substituted imine
        #"[$([CX3]([#6])[#6]),$([CX3H][#6])]=[$([NX2][#6]),$([NX2H])]", # substituted or unsubstituted imine
        "[NX3][CX3]=[CX3]", # enamine
        "C1C=C1", # cyclopropene
        "C1CC1", # cyclopropyl
        #"C=C-[N;H]" # cyclic enamine already captured in more general emanine
        "C1=CCC=C1" # cyclopentadiene
        ]
    for ni in range(len(forbidden_fragments)):
        if mol.HasSubstructMatch(MolFromSmarts(forbidden_fragments[ni])):
            violation = True
            # for indices in mol.GetSubstructMatches(MolFromSmarts(forbidden_fragments[ni])):
            #     display(Draw.MolToImage(mol, highlightAtoms = indices))
            # print(forbidden_fragments[ni])
            break
        else:
            continue
    return violation

Q: Consider the plot of bayesian optimizations using different representations of the molecules for a categorical model:
What can you tell about the difference in performance?
A: The BO search performs clearly better than the random search. The BO campaigns perform similarly with all the different representations used, most likely due to the categorical nature of our search space (it can be noted that Gaussian processes - the most commonly used underlying model for BO - are designed for continuous functions, not categorical searches).