Automatic design of functional
molecules and materials

Koji Tsuda
University of Tokyo / NIMS / RIKEN AIP

津田宏治 （つだ こうじ）
• 経歴
– 1972 京都生まれ– 1998 京大で博士号取得、旧電子技術総合研究所入所– 2000 ドイツGMD FIRSTで在外研究– 2001 産総研CBRCに配属– 2003-2004, 2006-2008 ドイツ・マックスプランク研究所– 2014 東京大学大学院教授• スキル– 機械学習、データマイニング、統計学、離散アルゴリズ
ム
2

第４のパラダイム：
データ中心科学• 経験科学・理論科学・計算科学（シミュレーション）に代わる21世紀の科学• 従来： 直感・経験に基づく仮説生成 ⇒ データに基づく検証• 第４： データからの知識発見による仮説生成⇒
データに基づく検証3

様々な分野からの要望が、データ科学に集中： Computational Thinking
化学

薬学
物理学

生物学

データ科学
環境・エネルギー

マクロ経済学
スマートグリッド4

Materials Genome InitiativeStrategic Report by JST (2013)

Automatic Materials Design

Experimental DesignMachine
Learning

Simulation(DFT etc)
DataExperiments

Materials Design Examples• Design of Si-Ge nanostructures (Phys Rev X2017)• De novo design of organic compounds withdesired absorption wavelength• Machine-learning-guided mutagenesis
platform for accelerated discovery of novel
functional proteins8Screening by first principles
calculations alone
Mat.
１

Mat.
２

Mat.
３

Mat.
４

Mat.
５

Mat.
６

Mat.
７

Mat.
８

Mat.
９Mat.
１０First Principles Calc.
Score Score Score Score Score Score Score Score Score Score123456789109

Bayesian Optimization(Jones et al., 1998)• Find best data points with minimum numberof observations• Choose next point to observe to discover thebest ones as early as possible

Bayesian Optimization (1)
Mat.
１Mat.
２

Mat.
３Mat.
４

First Principles Calc.
Score Score Score
123Mat.
５Mat.
６Mat.
７Mat.
８Mat.
９Mat.
１０

Bayesian Optimization (2)
Mat.
１Mat.
２

Mat.
３

Mat.
４Mat.
５

Mat.
６Mat.
７Mat.
８Mat.
９Mat.
１０

First Principles Calc.
Predicted Scores
Score Score Score Pred. Pred. Pred. Pred. Pred. Pred. Pred.
123Score Score Score Score Score Score Score
45678910Var.
4Var.
5Var.
6Var.
7Var.
8Predicted Variances

Var.
9Var.
10

Bayesian Optimization (3)
Mat.
１Mat.
２Mat.
３Mat.
８Mat.
４Mat.
５Mat.
６Mat.
７Mat.
９Mat.
１０First Principles Calc.
Score Score Score Score
123813

Bayesian Optimization (4)
Mat.
１Mat.
２Mat.
３Mat.
８Mat.
４Mat.
５Mat.
６Mat.
７Mat.
９Mat.
１０First Principles Calc.
Score Score Score Score Pred. Pred. Pred. Pred. Pred. Pred.
1238Score Score Score Score Score Score
456
7
910Var.
4Var.
5Var.
6Var.
7Var.
9Var.
1014

Where to observe next?
Current
Maximum

Measured Value
Explanatory Variable

15

Gaussian Process
Current
Maximum

Measured Value
Explanatory Variable16

Maximum probability of improvement
Current
Maximum

Measured Value
Explanatory Variable17

Alloy Structure Optimization (Phys Rev X, 2017)
Question: How to organize 16 alloy atoms (Si: 8, Ge: 8) to obtain the largest andsmallest interfacial thermal conductance?

311 151 5 9 134 8 12 162 6 10 14yz Lead
Descriptors: C168 = 12,870
Case 1 2 3 411 1 1 121 1 1 131 1 1 1… … … … …

7Si/Ge alloy region
5111…

6111…

7111…

8100…

90
10…

Lead

10 11 12 13 14 15 160 0 0 0 0 0 00 0 0 0 0 0 0
1 0 0 0 0 0 0… … … … … … …

Calculator: Atomistic Green’s Function (AGF): Phonon transmission
Evaluator: Interfacial Thermal Conductance (ITC)
Optimization method: Thompson Sampling (Bayesian Optimization)
Department of Mechanical Engineering, Thermal Energy Engineering Lab

18

Alloy Structure Optimization
ITCSi-SiSi-GeMax
MinOptimal structures were obtained by calculating only 3.4% of all candidates.
Department of Mechanical Engineering, Thermal Energy Engineering Lab

19

Superlattices Structure Optimization
Topic: Arrange 10-layer superlattices structure (5 layers Si + 5 layers of Ge) between
Si and Si to obtain minimal thermal conductance (1 layer thickness = 5.43 A)

Si/Ge superlattice region

LeadDescriptors:

C105 = 252
341
1
111
1
…
…

Case12111211311…
…
…
Best Structure: (1101010001)

5100…

6010…

7001…

8000…

Department of Mechanical Engineering, Thermal Energy Engineering Lab

Lead
9000…

10000…

20

COMBO: COMmon Bayesian
Optimization Python Libraryhttps://github.com/tsudalab/combo• Fast learning by random feature maps
• Automatic hyperparameter initialization & update

Ueno et al.,Materials Discovery,
2016.

Computational Time of COMBO
2×10 4#Features: 2000
#Features: 5000
scikit-learn

Time (sec)

1.5

10.500500

1000

1500

2000

Number of Experiments

2500

3000

Monte Carlo Tree SearchAtom assignment problem
Question: How to organize 16 alloy atoms (Si: 8, Ge: 8) to obtain the largest andsmallest interfacial thermal conductance?

3yz Lead

11 15
1 5 9 134 8 12 16
2 6 10 14
7

Si/Ge alloy region

Lead• Discrete black-box optimization:

Maximize f (x), x 2 {0, 1, k 1}NPNs.t.
`=1 I(x` = j) = nj , j = 0, . . . , k124

Search space is too large forBayesian Optimization(Jones et al., 1998)• Need to compute scores for all candidates10 3

Design time in seconds

10 2
10 1
10MDTS
Bayesian optimization

010 -1
10 -2
10 -3
10 -4
16

18

2022

24

Number of atoms in the structure26

Monte Carlo Tree Search• Candidates at leafs ofsearch tree
• Merit at leafs: Thermal
Conductivity• Merit at intermediate
node by Rollout– Random traverse to aleaf (engineered)– UCB score: average
merit + PremiumVariables at each node i• fi : Immediate merit by rollout• wi : Cumulative merit (downwardsum of fi )• vi : Visit count• UCB scorewiui =+Cvir2 ln vparentvi

MCTS steps• Initially, only root node exists• j-th level corresponds to xj• Tree grows by iterating the following 4 steps
Selection

Expansion

Simulation

Backpropagation

Root

Root

Root

Root
000

100

1

10

10

1000

10

1010110110001011100

10101110111010101001010

101

1. Selection• Reacha leaf byExpansion
following
Rootthe childwith largest UCB score

Selection
Root0011• In case of tie, choose
01
randomly

1001010102. Expansion
Expansion

Simulation

Root0011010• Generate k children
Root
under the selected
node01• If atom j is exhausted,01
nochild
is made• Initialize
01the variables
010110110001011100

10111010101001010

3. Simulation
Simulation
Root
001

10100110110001011100

10111010101001010

Backpropagation• Traverse randomly to a
Root
leaf
01
• Respectcomposition
constraints
01• Call first principle
01
calculation and obtainthermalconductivity
01

4. Backpropagation
Backpropagation
Root00110101• Visit count vi of allancestor nodes isincremented• Cumulative merit wi isalso updated

Experiments: Design time
• Total time = Design time + Simulation time# calculations vs success rate
• Success rate = Fraction of discovering the best in 100 runs

# calculations vs success rate
• Success rate = Fraction of discovering the best in 100 runs

Total time vs achieved TC (N=22):
fast design = more calculations

MDTS: Materials Design with Tree Searchhttps://github.com/tsudalab/MDTS• MCTS Implementation for atom assignment
problems with composition constraints• Parameter freede novo Molecular Generation
Machine learning (prediction)

Properties:
Binding Score
Bio Activity
ADMET…

ONONHOO

Generation of moleculeswith desired properties

SMILESONON

HOO

Cc3ccc(c2nc(CCCCO/N=C(CCC(O)=O)c1ccccc1)c(C)o2)cc3

Atom: {C, c, o, O, N, F, [C@@H], n, -, S,Cl, [O-],[C@H], [NH+],[C@], s,Br, [nH], [NH3+], [NH2+], [C@@], [N+], [nH+], [S@], [N-],[n+],[S@@], [S-], I, [n-], P, [OH+],[NH-], [P@@H], [P@@], [PH2], [P@],[P+], [S+],[o+], [CH2-], [CH-], [SH+], [O+], [s+], [PH+], [PH], [S@@+] }Bonds: {/,=, \# }Ring: {1,2,3,4,5,6,7,8,9}
Branch: {(, )}De novo Molecular Generation• Most existing methods make molecules bycombining predetermined fragments• De novo generation by deep neural networks– Variational autoencoder (Kusner et al., ICML 2017)– Recurrent neural network (Segler et al., Arxiv, 2017)• ChemTS (https://github.com/tsudalab/ChemTS)– MCTS + Recurrent Neural Network

Generation with VAE• Variational Auto Encoder (VAE).
Input:
SMILES

DNN(encoder)

Latent
DNNSpace (decoder)

Output:SMILES

Optimization inlatent space• Problem– Low chance of creating valid SMILES (0.4%〜5%)– SlowGómez-Bombarelli et al., “Automatic chemical design using a datadriven continuous representation of molecules,” Arxiv, 2016.
Kusner et al., "Grammar Variational Autoencoder," ICML, 2017.

ChemTS: MCTS and RNNrootCNO…

n[… Cn…

…

• Search tree of SMILES• Grow tree by MCTS
• Intermediate node =prefix of SMILES• Rollout: String
completion by RNN

RNNCompletion

4High chance of valid SMILES creation (>40%)
Fast generation

RNN：trained by 250,000 compounds• Input：String s1,…,sT• Output：Distribution of string P(y1),…,P(yT)• Output string is right-shift of input stringP(yt)Softmax activationht

Gated Recurrent
Unit x 2

ht+1

One-hot coding

stKeras & Tensorflow

MCTS(a) Selection(b) Expansion

1. Selection
R!"!"!$!#!%

2. ExpansionR!"!$!#!%!"!#!%!$!&!'

3. Simulation(c) Simulation(d) Backpropagation

RR!"!#!$!#!%!"!%!$!&

4. Backpropagation!"!$!#!%!"!%!$!&!'!#!$!#!%!'RNNUCB Score
2 ln 𝑉,-./01
𝑤"𝑢" =+𝐶
𝑣"𝑣"
AveragerewardPremium((Evaluate drug-likeliness score

𝐽 𝑚 = 𝑙𝑜𝑔𝑃 𝑚 − 𝑆𝐴 𝑚 − 𝑅𝑖𝑛𝑔𝑃𝑒𝑛𝑎𝑙𝑡𝑦(𝑚)Generating molecules with desired
absorption wavelength• Question: Can ChemTS generate synthesizable,
functional compounds within realistic
computational time?

Sumita et al., in preparation

47• ChemTS + TD-DFT (Gaussian)• 2 days computation on a 12-core server foreach of five target wavelengths• 6 molecules synthesized, 5 experimentally
validated
Target wavelength(nm)
Generated
DFT wavelength error< 20nmSynthesized
Confirmed inexperiments

200

300

400

500

600

646

757

629607

638

34

26

13

12

122

1

10

1211048

• Six compounds with known synthetic routes
and high oscillation strengths are synthesized
200nm

200nm

I300nm

IIIII

300nm

400nm

IV500nm

VVI49

UV spectroscopy• Five out of six reproduced TD-DFT wavelength

Tautomer50

51

Electronic origin of molecule I• n-π* excitation, not
HOMO-LUMOexcitation• Typical n-π*
excitations appear
at lower level (300400 nm)

52ChemTS Conclusion• ChemTS worked well in a proof-of-concept
study• Five functional compounds in 10 days with an
inexpensive server !
• Future work
– Automated detection of quantum simulation flaws– Synthetic route planning– High-throughput experiments (robots?)53

理研AIP分⼦情報科学チーム (PI: 津⽥宏治)• 人工知能を用いて、分子・物質を設計する手法を確立• 試行錯誤による設計よりも時間的・金銭的コストダウン• 医療・材料・環境に関連する物質設計への応用
医療材料・デバイス環境・がん免疫療法・ケモバイオ抗体

ナノデバイス

バイオマス分解

活性化
リンパ球

サイトカインなど放出

細胞傷害

がん細胞

54

タンパク質工学(protein engineering)• 抗体医薬（オブジーボ)、酵素など、タンパク質は広範な産業応用を持つ• タンパク質は、２０種のアミノ酸の鎖• 変異を導入することで機能強化・改変• K箇所に変異導入： ２０のK乗の可能性

55

通常の手法既存タンパク質

ランダム変異の入ったタンパク質

特性データ取得
優れたタンパク質を選択
56

AI支援タンパク質工学プラットフォーム
既存タンパク質

機械学習モデルを学習

ランダム変異の入ったタンパク質

機械学習によって設計された
変異の入ったタンパク質

特性データ取得

特性データ取得
優れたタンパク質を選択
57

機械学習による実験計画
物質
１物質
２物質
３物質
４物質
５物質
６物質
７物質
８物質
９物質
１０実験
特性値1

特性値2

特性値358

機械学習による実験計画
物質
１物質
２物質
３物質
４物質
５物質
６物質
７物質
８物質
９物質
１０

特性値3

予測値4

予測値5

予測値6

予測値7

予測値8

予測値9

予測値 10分散
4分散
5分散
6分散
7分散
8分散
9分散
10

実験
特性値1

特性値259

機械学習による実験計画
物質
１物質
２物質
３物質
４物質
５物質
６物質
７物質
８物質
９物質
１０

特性値3

予測値4

予測値5

予測値6

予測値7

予測値8

予測値9

予測値 10分散
4分散
5分散
6分散
7分散
8分散
9分散
10

実験
特性値1

特性値2

推薦順位

２位１位３位60

どれを観測すべきか？ 推薦順位をつける
現在の最大値

特性値

説明変数

61

ガウス過程による機械学習
現在の最大値

特性値

62

「改善確率」に基づく推薦
現在の最大値

特性値

３位

１位

２位

63

緑色蛍光タンパク質 (GFP)• 標的遺伝子に結合して、標的遺伝子の生体組織における発現を標識するために用いられる• 下村脩により発見 (2008年ノーベル化学賞)

64実証実験：AI支援プラットフォームを用いてGFPの機能を改変する• GFPの４箇所を改変 (65,72,77,203)• 候補数は16万 ： 網羅的実験は不可能• より波長が長く、蛍光が強いほど、大きくなるよ
うに統合特性値を定義• アミノ酸の説明変数：T-scale (Tian et al., 2007)
65

蛍光強度

ランダム変異蛍光変化率66

機械学習による推薦リスト(100個）

67

機械学習で設計した変異を導入したタンパク質

蛍光強度

12個の新規黄色蛍光タンパク質を発見！

蛍光変化率

68

今回発見した黄色蛍光タンパク質

GFP
69

まとめ：人工知能技術により物質設計は加速できる
⼤腸菌培養
蛍光観察
AIによる配列設計…MVSKGEELFT…
…MTSKVEGLFT…

定量評価

共同研究：
梅津光央 (東北大)斎藤裕 (産総研)
亀田倫史（産総研）