Applications of Artificial Intelligence to Computational Chemistry
Nicholas J. Browning
Supervisor: Prof. Dr. Ursula Roethlisberger
- Introduction
- Motivation
- EVOLVE
- GAML
- ML-MTS
Overview
1
Science, 2018, 361, 6400, pp. 360-365
Chemical Design
2
Chemical Space
Space defined by all possible molecules given a set of construction rules and boundary conditions
Chemical spaces are vast
J. Am. Chem. Soc., 2009, 131, 25, pp 8732–8733
ACS Chem Neurosci., 2012, 3, 9, pp 649–657.
3
> 10^{60}
Chemical Accuracy
Full CI
HF
DFT
MP2
Computational Cost
Accuracy
CC
\text{Experiment}
\infin
4
Motivating Remarks
Designing new compounds with interesting properties requires searching incomprehensibly large chemical spaces
Requiring chemical accuracy for predicted properties engenders the use of computationally expensive electronic structure methods
5
1. Searching Chemical Space
2. Improving Chemical Accuracy, for Less
Evolutionary Algorithms
Optimisation algorithms inspired by biological processes or species behaviour
6
Particle Swarm Optimisation
Genetic Algorithms
Machine Learning
Replace expensive electronic structure calculation with cheaper, data-driven regression model
\left < \Psi| \hat H | \Psi \right >\rightarrow f(\{Z, R\})
7
Motivating Remarks
8
1. Searching Chemical Space
2. Improving Chemical Accuracy, for Less
EVOLVE: An Evolutionary Toolbox for the Design of Peptides and Proteins
Nicholas J. Browning, Marta A. S. Perez, Elizabeth Brunk, Ursula Roethlisberger, in submission, JCTC, 2019.
- Introduction
- Motivation
- EVOLVE
- GAML
- ML-MTS
Haemoglobin (PDB 1A3N)
- DNA Signalling
- Acid-base catalysis
- atom/group transfer
- Electron transfer
- Redox catalysis
Design of proteins is attractive:
- Tune or change activity
- Improve solvent tolerance
- Increase thermal stability
Proteins fulfill a wide scope of biological functions ranging from regulation to catalytic activity:
Protein Design
Issues:
- Search space is large
- E.g 50-residue protein:
20^{50} \approx 1 \times 10^{65}
possible sequences
9
Short 20-mer polyalanine
Central 8 residues open to mutation (purple spheres)
EVOLVE Algorithm and Prototype
Side chain dihedrals discretised into set of 177 low-energy conformers [1]
1. Proteins, 2000, 40, pp 389-408
10
\vec{\chi} = \{\chi_i\}
f = \Delta E^\text{mutant}_{\text{system}} - \Delta E^{\text{ref}}_\text{system} = E^\text{mutant}_\text{system} - E^\text{ref}_\text{system} + \sum_i^N (E^{\text{ref}}_\text{aa}(i) - E^{\text{mutant}}_\text{aa}(i))
E_\text{system} = E^{FF99SB}_\text{system} + \Delta G^{GB}_{sol}
Amber FF99SB force field
Solvent effects modelled by implicit solvent model
Backbone stability measure
\epsilon = 80
\epsilon = 5
0
\epsilon = 10
Buried Protein Environment
Solvent Exposed
Intermediary
EVOLVE Fitness Function
11
94 Rotamers
177 Rotamers
177 Rotamers
-40.6 kcal/mol
-28.4 kcal/mol
-42.3 kcal/mol
EVOLVE Optimisation Curves
12
Fittest Individuals - Hydrophobic Environment
94 Rotamers
f = -28.4 kcal/mol
Buried Core Protein
13
Solvent Exposed Environment
Fittest Individuals - Solvent Exposed Environment
14
177 Rotamers
f = -40.6 kcal/mol
Homohelix Study
-40.6 kcal/mol
-28.4 kcal/mol
-42.3 kcal/mol
-28.1 kcal/mol
-28.0 kcal/mol
-27.9 kcal/mol
W04
15
Selection probabilities show consistent convergence towards a subset of chemical space
Selection Probabilities
16
Can use this subset to focus GA exploration on more interesting regions
55 Rotamer subset of residues: TRP, GLU, ASP, HIP, ARG
New fitness: -42.4 kcal/mol
Subspace Search
17
Outlook
J. Am. Chem. Soc., 2018, 140, 13, pp 4517–4521
Mutant protein solvent and pH tolerant
Thermostability greatly increased compared to wild-type
B1 Domain Protein G
18
1. Searching Chemical Space
2. Improving Chemical Accuracy, for Less
Genetic Optimisation of Training Sets for Improved Machine Learning Models of Molecular Properties
Nicholas J. Browning, Raghu Ramakrishnan, O. Anatole von Lilienfeld, Ursula Roethlisberger, JPCL, 2017
- Introduction
- Motivation
- EVOLVE
- GAML
- ML-MTS
-
Can the error be improved by intelligent sampling - is less more?
-
Are there systematic trends in “optimal” training sets?
Chemical Interpolation
19
Chemical Database
8 heavy-atom subset of the QM-9 database [2]
10 thermochemical and electronic properties computed with DFT/B3LYP
2. Sci. Data, 2014, 1.
Out-of-sample test error used as fitness metric
Complex optimisation problem:
Database QM-8
Kernel Ridge Regression
\approx 10^{1449}
C(10900, 1000)
20
Significant reduction in out-of-sample errors
Fewer training molecules required to reach given accuracy
random sampling
GA optimised
Learning Curves - Enthalpy of Atomisation
21
Direct Learning
Delta Learning
Systematic Trends - Moments of Inertia
22
Certain molecules are more representative
Systematic shifts in distance, property, and shape distributions
Systematic Trends - Property Distributions
23
Certain molecules are more representative
Systematic shifts in distance, property, and shape distributions
Motivating Remarks
1. Searching Chemical Space
2. Improving Chemical Accuracy, for Less
24
Machine Learning Enhanced Multiple Time Step Ab Initio Molecular Dynamics
Nicholas J. Browning, Pablo Baudin, Elisa Liberatore, Ursula Roethlisberger
- Introduction
- Motivation
- EVOLVE
- GAML
- ML-MTS
Molecular Dynamics
Provides a 'window' into the microscopic motion of atoms in a system
F_A = -\nabla_A E_0(\{R, r\})
Forces acting on (classical) nuclei calculated from a potential energy surface
25
Key Issues
- the accuracy of the potential energy surface : ML
- timestep used to integrate atomic forces : MTS
Multiple Time Step Integrators
Address the time step issue by observing that the forces acting on a system have different timescales
Separate "fast" and "slow" components of the total force, which can be integrated using smaller and larger timesteps
Slow forces typically expensive to compute therefore MTS yields a computational gain
26
MTS Algorithm
\gamma = (q_1, q_2, \dots, q_N, p_1, p_2, \dots, p_N)
iL = iL_\text{fast} + iL_\text{slow}
a(\gamma_t) = e^{iLt}a(\gamma_0)
e^{iL\Delta t} \approx e^{iL_{slow} \frac{\Delta t}{2}} \times e^{iL_\text{fast}\Delta t} (\Delta t) \times e^{iL_{slow} \frac{\Delta t}{2}}
e^{iL\Delta t} \approx e^{iL^\text{slow} (\Delta t/2)} \Big[ e^{iL^\text{fast} (\Delta t/n)} \Big]^n e^{iL^\text{slow} (\Delta t/2)}
Propagatation of phase space vector
Separate slow and fast components
phase space vector
Discrete time propagator obtained from Trotter factorisation
# initialise x, p
dt = 0.36
Dt = dt*n
Fslow = compute_slow_forces(x)
Ffast = compute_fast_forces(x)
for t in range(0, total_time, Dt):
p = p + 0.5*Dt*Fslow # Integrate slow components
for i in range(1, n): # Integrate fast components
p = p + 0.5*dt*Ffast
x = x + dt * p/m
Ffast = compute_fast_forces(x)
p = p + 0.5*dt*Ffast
Fslow = compute_slow_forces(x)
p = p + 0.5*Dt*Fslow # Integrate slow components
27
\delta t
Ab Initio MTS
28
Machine Learning Potentials
Replace expensive electronic structure method with data-driven kernel model
Delta-learning model naturally results in an MTS scheme
F^\text{fast} = F^\text{low}
F^\text{slow} = \Delta F^\text{ML}
F^\text{fast} = F^\text{low} + \Delta F^\text{ML}
F^\text{slow} = F^\text{high} - F^\text{fast}
Scheme I ML-MTS
Scheme II ML-MTS
\Delta F^\text{ML} = F^\text{high} - F^\text{low}
Exact slow forces
Cost of high level, evaluated infrequently
Approximate slow forces
Cost of low level
29
\Delta E^\text{ML}(M_I) = E^\text{high}(M_I) - E^\text{low}(M_I) = \sum_I\sum_J k(M_I, M_J) \alpha_J
\Delta F_L^\text{ML}(M_I) = -\frac{\partial \Delta E^\text{ML}}{\partial R_L}(M_I) = -\sum_I\sum_J \frac{\partial k(M_I, M_J)}{\partial M_I} \frac{\partial M_I}{\partial R_L} \alpha_J
Use a computationally efficient kernel method [3]:
- Generalisable
- No 2nd derivatives required, only kernel and descriptor gradients
Representation:
-
Atomic environments: Atomic Spectrum of London-Axilrod-Teller-Muto representation (aSLATM) [4]
MPI-ML code implemented into CPMD
Machine Learning Potentials
3. Chem. Phys., 2019, 150, 064105
4. ArXiv, 2017, https://arxiv.org/abs/1707.04146
30
Training Data
Trajectories of small isolated N-water molecule clusters
Fast components = LDA
Slow components = PBE0 - LDA
Predictions made on liquid water
31
Results - Errors
32
Results - Errors
Significant improvement upon LDA - forces almost identical to PBE0
33
Results - Time Step and Speedup
Significant increase in outer timestep wrt. standard MTS
Resonance effects limit maximum outer timestep
Encouraging increase in overall speedup
Speedup dictated by number of SCF cycles to converge wfn
\epsilon = \frac{1}{N}\sum_i \frac{E_i - \bar{E}}{\bar E}
34
Results - Properties
MTS (n=4dt) identified as yielding properties with excellent agreement wrt. VV-PBE0 in previous work [5]
Both scheme I and II ML-MTS result in significantly improved radial distributions
5. J. Chem. Theory Comput., 2018, 14, 6, pp 2834–2842
35
1. Searching Chemical Space
Conclusions
- EVOLVE is capable of searching through inordinately large chemical spaces
- Meaningful chemical space insights can be gained from the EVOLVE procedure
- Potential for experimental reach
36
2. Improving Chemical Accuracy, for Less
-
Within a specified chemical space, the accuracy of ML models can be significantly improved with intelligent sampling
-
kernel methods can provide a significant increase in timestep used in MTS integrators
- Maintain good accuracy
- Significantly reduce computational cost
Acknowledgements and Thanks
Dr. Raghu Ramakrishnan
Dr. Esra Bozkurt
Dr. Marta A. S. Perez
Dr. Pablo Baudin
Dr. Elisa Liberatore
Prof. Ursula Roethlisberger
37
Prof. O. Anatole von Lilienfeld
Prof. Clemence Corminboeuf
Dr. Albert Bartók-Pártay
Dr. Ruud Hovius
Acknowledgements and Thanks
38
Thermostability Measure
f = \Delta E^\text{mutant}_{\text{system}} - \Delta E^{\text{ref}}_\text{system} = E^\text{mutant}_\text{system} - E^\text{ref}_\text{system} + \sum_i^N (E^{\text{ref}}_\text{aa}(i) - E^{\text{mutant}}_\text{aa}(i))
E_\text{system} = E^{FF99SB}_\text{system} + \Delta G^{GB}_{OBC}
Backbone stability measure
M. A. S. Perez, E. Bozkurt, N. J. Browning U. Roethlisberger, in preparation
39
Thermostability Measure
40
Manual Mutations
\epsilon = 80
D03-E05 -> R*-E05... = -30.1 kcal/mol
-40.6 kcal/mol
-28.4 kcal/mol
-42.3 kcal/mol
D03-E05 -> E*-E05... = -39.2 kcal/mol
\epsilon = 10
W04-W04-W04-W04-F04-Q05-W04 -> {8 x W*} = -28.2 kcal/mol
41
EVOLVE: MOGA
42
Amino Acids
43
GA-ML Small Optimal Training Sets
44
GA-ML Distance Distributions
45
\boldsymbol{M}_I^{(2)}(r) = \frac{1}{2} \sum_{J\ne I} Z_IZ_J \frac{1}{\sigma_r \sqrt{2\pi} r^6} e^{-\frac{\left \lVert r - R_{IJ} \right \rVert ^2}{2\sigma_r^2}}
\boldsymbol{M}_I^{(3)}(\theta)= \frac{1}{3}\sum_{I \ne J \ne K} Z_I Z_J Z_K \frac{1}{\sigma_\theta \sqrt{2\pi}} e^{-\frac{\left ( \theta - \theta_{IJK} \right )^2}{2\sigma_\theta^2}} \frac{1 + \cos\theta \cos\theta_{JKI}\cos\theta_{KIJ}}{\left (R_{IJ} R_{IK} R_{KJ}\right)^3}
\boldsymbol{M}_I = \left \{ \boldsymbol{M}^{(1)}_I, \boldsymbol{M}^{(2)}_I, \boldsymbol{M}^{(3)}_I\right \}
SLATM Parameters
46
R_c = 4.8 A
\theta_\text{min} = -0.35 \ \text{rad}
\theta_\text{max} = 3.5 \ \text{rad}
\text{grid resolution}= 0.1A
\text{angular grid resolution}= 0.1
\ \text{rad}
Two-body term:
Three-body term:
\sigma_r = 0.05
\sigma_{\theta} = 0.05
Thermostability
47
Applications of Artificial Intelligence to Computational Chemistry
By Nick Browning
