Proceedings of the

4th Biosupercomputing Symposium- International Symposium for Next-Generation Integrated Simulation of Living Matter（ISLiM）-

December 3-5, 2012
Tokyo International Forum ( Hall D7 )"Next-Generation Integrated Simulation of Living Matter (ISLiM)" programcommissioned by Ministry of Education, Culture, Sports, Science and Technology ( MEXT )Proceedings of the

4th Biosupercomputing Symposium-International Symposium for Next-Generation Integrated Simulation of Living Matter ( ISLiM )-December 3-5, 2012
Tokyo International Forum（Hall D7）"Next-Generation Integrated Simulation of Living Matter（ISLiM）" programcommissioned by Ministry of Education, Culture,
Sports, Science and Technology（MEXT）

4th Biosupercomputing Symposium (December 3 - 5, 2012)

AgendaDecember 3 (Mon)
10:00-10:20

Opening Talk

Koji KAYA (RIKEN)(Session Chair: Shu TAKAGI, RIKEN/University of Tokyo)

10:20-11:10

ISLiM Keynote: Current status of software development in ISLiM and its future

11:10-12:00

What can Multiscale Models do for Precision Medicine?

12:00-13:00

Ryutaro HIMENO (RIKEN)

Grace Peng (National Institute of Health)Lunch break(Session Chair: Ryutaro HIMENO, RIKEN)

13:00-13:50

Development of Multiscale Thrombosis Simulator

13:50-14:40

Multiscale simulations of tumor angiogenesis with therapeutic applications

14:40-15:10
15:10-16:00

Shu TAKAGI (RIKEN/University of Tokyo)

Aleksander Popel (Johns Hopkins University)
Coffee break

Multi-scale, multi-physics heart simulator "UT-Heart" for heart research
Seiryo SUGIURA (University of Tokyo)(Session Chair: Hideo YOKOTA, RIKEN)

16:00-16:50

Is eukaryotic cell biology ready for supercomputing?
A perspective from the Virtual Cell projectIon I. Moraru (University of Connecticut Health Center)

December 4 (Tue)(Session Chair: Makoto TAIJI, RIKEN)

10:00-10:50

K computer Keynote: Real Peta-scale computing, stairway to Exa

10:50-11:40

Optimization of Life-Science applications on the K computer

11:40-12:30

Large Scale Biomolecular Modeling with IBM Blue Gene

12:30-13:30

Motoi OKUDA (Fujitsu)

Yousuke OHNO (RIKEN)
Ruhong Zhou (IBM)
Lunch break(Session Chair: Ryutaro HIMENO, RIKEN)

13:30-13:55

Research activity of Cell Scale Simulation Team

13:55-14:20

Dissection of regulatory mechanisms for metabolic systems by quantitative
imaging mass spectrometry

Hideo YOKOTA (RIKEN)

Makoto SUEMATSU (Keio University)

4th Biosupercomputing Symposium (December 3 - 5, 2012)(Session Chair: Shin ISHII, Kyoto University)

14:20-15:10

Closing the loop: simulation of the whole sensory-motor neural network in action
Kenji DOYA (Okinawa Institute of Science and Technology〈OIST〉)

15:10-15:40

Coffee break

15:40-16:30

Multi-scale modelling of the motor system underlying
goal-directed behaviour in a vertebrate model organism

16:30-17:20

Brain-scale neuronal network simulations on K

(17:30-19:00)

Social (optional)

Sten Grillner (Karolinska Institutet)

Markus Diesmann (Forschungszentrum Jülich)

December 5 (Wed)(Session Chair: Akinori KIDERA, RIKEN/Yokohama City University)

10:00-11:00

Molecular machines and nuclear processes studied by coarse-grained
molecular simulations

11:00-12:00

Biomolecular simulations under cellular environment

12:00-13:10

Shoji TAKADA (Kyoto University)
Yuji SUGITA (RIKEN)
Lunch break(Session Chair: Mitsunori IKEGUCHI, Yokohama City University)

13:10-14:10

Large-scale simulations of biomolecular machines

14:10-15:10

Use of Modeling and Simulation in Drug Discovery and Development:
It is all about the Questions

15:10-15:40

Coffee breakKarissa Sanbonmatsu (Los Alamos National Laboratory)

Sandra R. B. Allerheiligen (Merck)(Session Chair: Satoru MIYANO, University of Tokyo)

15:40-16:05

Supercomputing for Next-Generation Cancer Research

16:05-16:30

Large-scale protein-protein interaction network prediction
by an exhaustive rigid docking system MEGADOCK

Satoru MIYANO (University of Tokyo)

Yutaka AKIYAMA (Tokyo Institute of Technology)(Session Chair: Yutaka AKIYAMA, Tokyo Institute of Technology)

16:30-16:55

Supercomputing accelerates genomic medicine

16:55-17:20

Towards efficient improvement of transcriptional circuit models by Life Science
Data Assimilation System (LiSDAS)

17:20-17:30

Closing Talk

Tatsuhiko TSUNODA (RIKEN)

Tomoyuki HIGUCHI (Institute of Statistical Mathematics)
Ryutaro HIMENO (RIKEN)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Table of Contents
Title

Speaker

Preface

Koji KAYA(RIKEN)

Page
1

Keynotes
ISLiM Keynote:
Current status of software development in ISLiM and its future

Ryutaro HIMENO

K computer Keynote:
Real Peta-scale computing, stairwayto Exa

Motoi OKUDA(RIKEN)(Fujitsu)

5

13Invited Talks
What can Multiscale Models do for Precision Medicine?

Grace Peng
(National Institute of Health)Multiscale simulations of tumor angiogenesis
with therapeutic applications

Aleksander PopelIs eukaryotic cell biology ready for supercomputing?
A perspective from the Virtual Cell projectIon I. Moraru

Large Scale Biomolecular Modeling with IBM Blue Gene

Multi-scale modelling of the motor system underlying
goal-directed behaviour in a vertebrate model organism

Large-scale simulations of biomolecular machines
Use of Modeling and Simulation in Drug Discovery andDevelopment: It is all about the Questions(Johns Hopkins University)(University of Connecticut
Health Center)

Ruhong Zhou
(IBM)

Sten Grillner(Karolinska Institutet)

Karissa Sanbonmatsu(Los Alamos National
Laboratory)

Sandra R. B. Allerheiligen(Merck)

21

25

29

33

37

39

43

4th Biosupercomputing Symposium (December 3 - 5, 2012)

ISLiM Talks
Development of Multiscale Thrombosis Simulator
Multi-scale, multi-physics heart simulator "UT-Heart" forheart research
Optimization of Life-Science applications on the K computer
Research activity of Cell Scale Simulation Team

Shu TAKAGI(RIKEN/University of Tokyo)

Seiryo SUGIURA(University of Tokyo)

Yousuke OHNO(RIKEN)

Hideo YOKOTA(RIKEN)

Dissection of regulatory mechanisms for metabolic systems
by quantitative imaging mass spectrometry

Makoto SUEMATSU

Closing the loop: simulation of the whole sensory-motor
neural network in action

Kenji DOYA

Brain-scale neuronal network simulations on K
Molecular machines and nuclear processes
studied by coarse-grained molecular simulations
Biomolecular simulations under cellular environment
Supercomputing for Next-Generation Cancer Research
Large-scale protein-protein interaction network prediction
by an exhaustive rigid docking system MEGADOCKSupercomputing accelerates genomic medicine
Towards efficient improvement of transcriptional circuit
models by Life Science Data Assimilation System (LiSDAS)(Keio University)(Okinawa Institute of Scienceand Technology〈OIST〉)

Markus Diesmann(Forschungszentrum Jülich)

Shoji TAKADA(Kyoto University)

Yuji SUGITA(RIKEN)

Satoru MIYANO(University of Tokyo)

Yutaka AKIYAMA(Tokyo Institute of Technology)

Tatsuhiko TSUNODA(RIKEN)

Tomoyuki HIGUCHI(Institute of Statistical
Mathematics)

47
53
59
65
71
77
83
91
97
103
111
117
125

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Preface

Koji Kaya
Program Director
Computational Science Research Program, RIKENProfile:1966, Doctor of Science, Tokyo University
1966, Research Staff, Theoretical Organic Chemistry Lab., RIKEN
1970, Associate Professor, Department of Chemistry, Tohoku University
1973, Member of Technical Staff, Bell Telephone Laboratories, U. S.
1981, Professor, Department of Chemistry, Keio University
1999, Director in General, Institute for Molecular Science
2004, Directors of Discovery Research Institute, and Wako institute, RIKEN
2006 - Program Director, Next Generation Computational Science Research Program1

4th Biosupercomputing Symposium (December 3 - 5, 2012)

In 2006, The Next Generation Super Computer Project started at RIKEN under the leadership of President
Ryoji Noyori sponsored by MEXT Japan. Since then, this project experienced and overcome many hurdles bothin technical and political points. In October 2011, we reached the final goal, i.e., establishment of 10 peta-flops
computation.
In parallel with the development of the super computer (hard ware), the grand challenge program forthe software development for peta-flops computation in life science started in October 2006 which is named
Integrated Simulation of Living Matter (abbreviated as ISLiM). ISLiM program consists of 6 teams including 4
computer simulation teams, a data analysis team and software sophistication technology team. As a result of hard
works by the young reseachers, most of the developed software programs finally reached to the final stage of thedevelopment worthwhile to challenge to peta-flops computation.
In this symposium, we focus our attention on the reports and discussions of the results of our softwaredevelopment and possible applications of the software programs to life sciences including drug discovery, medical
applications in addition to the sophisticated lectures by the invited speakers from abroad.
Finally, we wish to express our sincere gratitude to the 7 invited speakers from abroad who kindly accepted
our invitation to ISLiM 4-th symposium for their great contribution to make this symposium fruitful for all theparticipants.
Koji Kaya
Program Director
Computational Science Research Program
RIKEN

2

Keynotes

4th Biosupercomputing Symposium (December 3 - 5, 2012)

ISLiM Keynote:

Current status of software development
in ISLiM and its future
Ryutaro Himeno
Group Director, Research and Development ofIntegrated Simulation of Living Matter, RIKENProfile:Ryutaro HIMENO received the B.E. and the M.E. from Kyoto University, Kyoto, Japan in 1977 and 1979,
respectively. He received Doctor of Engineering degree from the University of Tokyo in 1988. In 1979, he joined
Nissan Motor Co., Ltd., Yokosuka, Japan, where he has been engaged in the research of applying Computational
Fluid Dynamics analysis to the car aerodynamic development. From 1984 to 1986, he served as a researcher
at the Institute of Space and Astronautical Science, Tokyo, Japan. In 1999, he joined RIKEN. He is now thedirector of Advanced Center for Computing and Communication, Deputy Program Director of Computational
Science Research Program and Group Director of Research and Development of Integrated Simulation of Living
Matter at RIKEN. He is a visiting Professor at the Kobe University, Hokkaido University and Tokyo Denki
University. He is also well known researcher in Baseball Physics.
He was a winner of 2006 Gordon Bell Prize (Honorable Mention, Peak Performance, with Dr. Tetsu Narumi,
et al) and Computer Visualization Contest in 2000 by Nikkei Science. He received both JSME Computational
Mechanics Award and Computational Award by Japan Association of Computational Mechanics in 2011 as well
as JSME Computational Mechanics Achievement Award in 1997 and JSME Youth Engineer Award in 1988. He
was also awarded Paper Award by NICOGRAPH in 1993, Giga FLOPS Award by CRAY Research Inc. in 1990
and other awards.

5

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Our project ISLiM started in 2006 and six years havepast. The project will end in the end of March, 2013.
This international symposium was planed to showwhat we have developed and achieved in past six
years and to announce what we are going to do fromnow on.

2) Innovative simulated results shall be shown using
full system of K computer.
3) All codes shall be available on K computer whenusers want to use
3. Current status
Three codes: ZZ-EFSI, cppmd and UT-Heart achieved
effective performance larger than 1 PetaFLOPs on K
computer. 13 codes showed good liner scaling up to
more than 10,000 nodes on K computer.
Computation using K computer is currently going on
and the results will be shown from each team in this
symposium. Major codes we have developed can be
down-loadable from our web site and the number ofavailable codes are increasing. We will complete thelibrary at the end of March, 2013.
We will keep the library for users who want to use
and take care of those users even after we finish our
project.

2. What we plane in 2006
In 2006, we started a grand challenge project called
ISLiM for K computer to demonstrate its performance.
The ISLiM stands for Integrated Simulation of
Living Matter to reproduce life phenomena on asupercomputer for understanding them and developing
new medicine or new medical treatments. We have 6
research teams: Molecular scale team, cell scale team,organ and body scale team, data analysis fusion team,brain and neural system team and HPC team. We
started a high performance software package for life
science for K computer which contains 31 application
software.

References
1. http://www.csrp.riken.jp/index_e.html

The target was as follows:
1) At least 2 application codes achieved effective one
PetaFLOPs on K computer.

Figure1 project target and structure of research teams

6

7

Applications

Buildings

Researchbuilding

Computerbuilding

Next-Generation
Integrated
Life Simulation

Next-Generation
Integrated
Nanoscience
Simulation

System

FY2008

Prototype,
Prototype,evaluation
evaluation

FY2009Construction
Construction

Construction
Construction

4th Biosupercomputing Symposium

DesignDesignDesignDesignDevelopment,
Development, production,
production, andand evaluation
evaluation

FY2011

1

3

TuningTuning andandimprovement
improvement

FY2012

Verification
Verification

Verification
Verification

Production,
Production, installation,
installation,andand adjustment
adjustment

FY2010Life Science Grand Challenge:
ISLiM

Development,
Development, production,
production, andand evaluation
evaluation

Detailed
Detailed design
design

FY2007

Conceptual
design

FY2006

Next-generation Supercomputer R&D project and Life
Science Grand Challenge

4th Biosupercomputing Symposium

RIKEN
himeno@riken.jp

Group Director, Research and Development ofIntegrated Simulation of Living Matter, RIKEN

Ryutaro Himeno

Current status of software development in
ISLiM and its future10-3~-2
Macro

100

Organ/Wholebody scale

Cell
scale

Basic principal(Equations)

Analytical
approach

Molecule
scale

New

Experimental data

Fusion of Data based methodology
approach and Anaｌ
ytical approach

Data-Basedapproach toanalysis2

4•Understand the life
phenomena:
reproducing life
program on computer•Develop new medicine
and medical equipment

Peta-scale Life
Science Software
Package

Approachcombination of analytical approach and Data-driven approach on Peta-scalecomputerGrand Challenge in Life Science at 1st stage

Meso

10-5

Goal
Understand the life phenomena: reproducing life program on computer
Develop new medicine and medical equipment

Micro

10-8

Movement of many particles will becalculated on Supercomputer insome day.
Life phenomena will become predictable while current
biology is descriptive science
We decided to develop our own application codes
optimized for the newly developed K computer.

Life Phenemena are so complicated and still
unknown.
But in molecular scale, they are governed bytheories in Physics

Motivation

4th Biosupercomputing Symposium (December 3 - 5, 2012)8

Abstraction by modeling→out

Experimental data

Diffusion coefficient

Mitochondria

Aquaporin

Ion channel

Cell membrane

Organelle velocity vector

Metabolic pathway

Support of integrated simulation

Construction of a differentiated cell model
Extension of each program
Coupled analysis of various simulation

Live cell model

Realization of integrated simulationin the intracellular coordinate systemin→y=f(x)+g(x)exp(az)･･･

Fundamental research forcell model construction

Differentiated cell model

Velocity of matter

Fusion of Data based
approach and Anaｌ
ytical approach

Brain & Neural System

BrainsystemNerve cell
local circuit

TGCTCTTATCTTCCTCCCACAGCTCCT
GGCCCATCATGGCTGGTGTGGCTAATGCCCGCTATTAAACTTAATAAAAAACATTTATTTTCATTGCAATATTTTACTGATATTATACAGGGGGCATGGTTTGACT……

Data-Based approach toanalysis

Local simulationHeart

Pancreatic cell

Pancreas

Cardiac muscle cell

Erythrocyte

Hepatocyte

Blood disc

Lever

Basic principal(Equations)

Organ/Wholebody scale

Cell
scale

Ion
channelModeling beyond scale

Molecule scale

Analytical
approach

Added Brain & Neural System Research

Grand Challenge in Life Science at 2nd stage

75(1) Bonding

AMP(2) Structure change(3) Chemical reaction

Molecular Dynamics

Amino acid as fundamental degrees of freedom
Empirical potential obtained from database6

Larger size, longer time,and less accuracy

Coarse Grained Model.

Whole body voxel model

Surface modelMulti physics model

velocity

Muscle force

Stress-strain

Blood tube
blood
muscleWhole bodycoronal

Leg model

Eyeball model

8

Cerebral artery model

Local precise model

Organ/Whole body simulation

Calc. of molecular force field

Smaller size, shorter time,and greater accuracy

Medical images

Ab-initio calculation

Quantum Chemistry

ADP x 2

Example: Phosphoryl transfer reaction of adenylic asid kinase

ATP

Dynamics of molecules

Method for the Molecule-Scale Research: Combination of the following three models

Simulating protein function in life
phenomena to reveal important
subjects in molecular and cellbiology

Objectives

Molecule-Scale Research

4th Biosupercomputing Symposium (December 3 - 5, 2012)

9

axon

Groth corn

Symposium

Grand challenge targets in longterm
4th Biosupercomputing

Nerve cell

Nerve process

High precision neuron simulation
with development & plasticity
mechanism （No. of molecules: 10,102 compartments）

Neuron simulation

Investigate Neurodegenerative disorder and contribute to11
its medical treatment

Establish a total model integrating input, decision andcommand of brain.

Whole brain simulation
Challenge image processing of
105pixel picture byamphiblestrode modelwith106cells & visual cortex withcolumns & 105 cells

Local circuit simulation

High precision simulation of
information processing &plasticity in brain cortex column(10105cells,107synapse)

Challenge brain and neural system simulation combining scales: neuron, local
circuit and brain. This team includes Prof. Kanzaki, Dr. Usui, Dr. Diesman andProf. Doya.

4th Biosupercomputing SymposiumEx) Ultrasound propagation
simulation for HIFU

and more….

9

In-silico analysis ofdrug-target gene

PPARαModeling of functions

Peta scale computing

Peta-scale
data-analysis

10

Organ /
whole body

DataAnalysis

others

4th Biosupercomputing
Symposium

Academic and industrial usage

Feedback Software family optimized for the supercomputer
after
applying
Improve by
Improve software for commercial use
real
ourselves
phenomenaEasy to use from user’s pointof view

Common libraries + Core software

Support to derive higher performance

Cell scale

GUI / Workflow

Molecular
scale

12

Softwarehouse

Support toutilizedeveloped
codes
Grand Challenge
Goal

Commonsoftwareplatform andoptimization

Each softwaredevelopment

・Fusion of simulation models andexperimental data based on data
assimilation- Automation of large-scale modeling- Information technology for fusion ofsimulation and data analysis

HPC Team for development of common software platform and optimizationon the new supercomputer

Estimation of gene network

1049 genes- Order-made medicine and therapy based on
optimized peta-scale data analysis- Peta-scale analysis of gene network for drug
target discovery and validation

Expected outcome- Higher dimensions: from gene to exon
- Heterogeneity: array, structure, dynamics,…
- Invisibility: deficit of information

Exponential increase of experimental data

Data-driven large-scale
simulation technology

Brain & Neural System Team lead by Prof. Ishii (started
in Oct. ’08)

By Prof. Hisada andProf. Sugiura

Planned simulation based on precise Human model

Overview of data analysis research
Fusion of large-scale heterogeneous data
analysis and life-science simulation

4th Biosupercomputing Symposium (December 3 - 5, 2012)10

Whole genome association study, data assimilation,
prediction of protein-protein interaction, genome
sequencer data processing, Haplotype wholegenome association study

Parallel middleware, parallel visualization software,
high-speed software core library

Data based
analysis

HPC

80000

90000

Neural sim., whole brain sim., cortical micro circuit
sim.,

Brain & neuralsystem

0

10000

20000

30000

40000

50000

60000

70000

PFLOPsNumber of Nodes

4th Biosupercomputing Symposium

Test operation of K computer

4th Biosupercomputing Symposium

Heart sim., Lung sim., HIFU sim., Fluid-structure sim. 6
for blood flow,

Organ/wholebody scale

4

9

5

1

Voxel-based multi-compartment transport-diffusion
simulator

Cell scale

Available No. of Nodes

9

MD, Quantum Chemistry, Coarse Grain MD

0

2

4

6

8

15

10

12

13

No. ofapplication
SoftwareMolecule scale

Applications
Platypus-MM/CG
Platypus-REIN
MARBLE
CafeMol
ProteinDF
Platypus-QM/MM-FE
Platypus-QM
Platypus-QM/MM
Platypus-CGM/MM
RICS
ZZ-EFSI
ZZ-DOSEZZ-HIFUZZ-LUNG
UTHeart
ParaHaplo
NGS analyzer
ExRAT
SiGN-BN
SiGN-L1
SiGN-SSM
SBiP
LiSDASMEGADOCK
NEST
CMDN
VSM
NeuroMorphoKit
IOSSIM
cppmd
VLSVL

No. of Cores

Life Science Software Package

0

100

1000

10000

1000001000000

1

2

3

4

Phase II < 8192

5

2011年12月

2012年2月

2011年11月

2012年1月

2011年10月
2011年9月

2011年9月

2011年8月

2011年7月

2011年6月

2011年5月

2011年5月

4th Biosupercomputing Symposium

UTHeart

ZZ-HIFU

ZZ-EFSI

MARBLE

PlatypusMM/CG
Platypus-REIN

cppmd

IOSSIM

NEST

SiGN-BNMEGADOCKLiSDAS

14

16

Cppmd ~
4Peta
7
Goal: > 1 Pata FLOPS

6

24/27 running on K

UT-Heart~ 3Peta

ZZ-EFSI :4.5Peta

882944nds

Hybrid parallelization
with MPI + OpenMP

12288nds

Phase III on K

1024nds

Code development history: scalbility

Phase I < 256

Current status of development

4th Biosupercomputing Symposium (December 3 - 5, 2012)2012年6月

2012年5月

2012年4月

2012年3月

2012年2月

2011年4月

2011年3月

2011年2月

2010年12月

2010年12月

2011年1月

2010年11月

2010年10月
2010年9月

PFLOPS

11

4th Biosupercomputing Symposium

0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

ZZ-EFSI

42.7%

4th Biosupercomputing Symposium

UT-Heart

27.7%

82,944 nodes

cppmd

38.0%

Measured performance on K computer using fullsystem

1

10

100

1,000

10,000PFlopsNodes

100,000

Max Node No. of Linear Speedup

19

17

End of July,2012
End of April, 2012

Parallelization ratio p (%)

4th Biosupercomputing Symposium

Weak Scaling on K: ZZ-EFSI•Using FFT to solve Poisson equation on PC•Multi-color SOR on PC Cluster•On K computer which has 80,000 CPUs
Implicit solver  explicit solver using quasicompressible formulation
less communication

Changing solvers in ZZ-EFSI

98.6

98.8

99

99.2

99.4

99.6

99.8

100

Parallel efficiency

20

18

4th Biosupercomputing Symposium (December 3 - 5, 2012)

12

4th Biosupercomputing Symposium• We have developed wide variation of high performance
application codes in Life Science for K computer, frommolecular scale to organ/whole body scale as well as data
analysis and brain/ neural system.
• 12 codes linearly scale more than 10,000 nodes of K
• ZZ-EFSI achieved 4.5 PetaFLOPs and 2 other codesachieved more than 3 PetaFLOPs.
• Currently, we are analyzing computed results right now.
• The codes are available from our web:
http://www.csrp.riken.jp/index_e.html

Summary

Organization 13(RIKEN, U. of Tokyo, Osaka U., Kyoto U. Tohoku U.,
Tokai U., JAIST, Chiba U., Keio U. Yokohama CityU., Kobe U., TiTch, ISM)
About 200 person, including 60 postDoc

23

21

For your references

Workshop in 2009

Our homepage:http://www.csrp.riken.jp/index_e.html
Youtube:Approach to the Life Sciences Grand Challenge-Next
Generation Integrated Simulation of Living

Workshop in 2008

Contributors

24

22

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

K Computer Keynote:

Real Peta-scale computing, stairway to Exa

Motoi Okuda
Executive Architect, Technical Computing Solutions Unit, Fujitsu Ltd

Profile:Motoi Okuda Ph.D.
Executive Architect, Technical Computing Solutions Unit, Fujitsu Ltd.
Work Experience· Planning, marketing and support of High Performance Computer· Supervising Fj activities of national projects such as National Grid Project and Next Generation Supercomputing
Projects· Developing computational science and engineering applications, such as fluid dynamics nuclear code,crashworthiness, and molecular dynamics· Project leader of Parallel Computing Center in Fujitsu Lab.
Education· Ph.D., Information Science from Japan Advanced Institute of Science and Technology, Ishikawa.
· Master of Nuclear Engineering from Nagoya University.
Membership in Professional Societies includes· Executive board member of The Japan Society for Computational Engineering and Science (JSCES)

13

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
FUJITSU and RIKEN have been working together todevelop K-computer.
After six- year project, the K-computer entered into
service at the end of September and many users are
now using the huge amount of computing power We
are now expecting innovative outcomes of the mostadvanced K computer use.
My talk aims to cover the outline of the K computer
project, its design concept, how to implement
the system, impacts of adopted technologies onapplication execution, and some examples ofapplication performance achievement. K computer’sCPU, SPARC64TM VIIIfx, and Tofu interconnect are

the key technologies for achieving the project targets
such as performance, reliability, operability and power
consumption. I will talk about these features and an
example of performance on application program.
The talk also refers to our new challenges to thenext step, Exa-scale computing, based on lessons
we learned from K computer project and technology
trends. Successful project management is crucial topromoting such a huge project. So I would like tomention it briefly for the future project planning.
Japanese government has just begun their new
challenge to the future supercomputer system. I will
also refer to the outline of Japanese activities.

14

4th Biosupercomputing Symposium (December 3 - 5, 2012)

15

16

Communication

Communication

4th Biosupercomputing Symposium

7

Courtesy of RIKEN
Copyright 2012 FUJITSU LIMITED* : SC ‘11 Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis, Article No.
1 ,First-principles calculations of electron states of a silicon nanowire with 100,000 atoms on the K computer New high speed allreduce algorithm : Communication throughput 3.2GB/s.
 Application performance : 3.08Pflops (43.6% efficiency)(107,292 atoms, 55,296 nodes)

Computation

19%

Modified for K computer

Tuning for computation
Tuning for communication(8,000 atoms, 256 nodes)

Computation

32%

Original

receiving 2011 Gordon Bell Award* New allreduce algorithm for Tofu interconnect made a great contribution to

Tofu Interconnect - New allreduce algorithm -4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)17

Invited Talks

4th Biosupercomputing Symposium (December 3 - 5, 2012)

What can Multiscale Models do for Precision
Medicine?
Grace C.Y. Peng
Program Director, National Institute of Biomedical Imaging and Bioengineering(NIBIB), National Institutes of Health (NIH)

Profile:Grace C.Y. Peng received the B.S. degree in electrical engineering from the University of Illinois at Urbana, theM.S. and Ph.D. degrees in biomedical engineering from Northwestern University. She performed postdoctoraland faculty research in the department of Neurology at the Johns Hopkins University. In 2000 she became theClare Boothe Luce professor of biomedical engineering at the Catholic University of America. Since 2002, Dr.
Peng has been a Program Director in the National Institute of Biomedical Imaging and Bioengineering (NIBIB),
at the National Institutes of Health. Her program areas at the NIBIB include mathematical modeling, simulation
and analysis methods, and next generation engineering systems for rehabilitation engineering, neuroengineeringand surgical systems. In 2003, Dr. Peng lead the creation of the Interagency Modeling and Analysis Group (IMAG),
which now consists of program officers from ten federal agencies of the U.S. government and Canada (www.
imagwiki.org/mediawiki). IMAG has continuously supported funding specifically for multiscale modeling (ofbiological systems) since 2004. IMAG facilitates the activities of the Multiscale Modeling (MSM) Consortiumof investigators (started in 2006). Dr. Peng is interested in promoting the development of intelligent tools andreusable models, and integrating these approaches in engineering systems and multiscale physiological problems.
21

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction – What is Precision
Medicine?
In 2011, the National Academy of Sciences (NAS)published a report entitled, “Toward Precision
Medicine: Building a Knowledge Network for
Biomedical Research and a New Taxonomy of
Disease.”1This report has spurred many talks anddiscussions in the scientific arena in the UnitedStates to examine how we can manage the explosionof disease-relevant data, and better integrate theknowledge gained from basic biomedical researchwith medical histories and health outcomes ofindividual patients. The NAS Committee concluded
it is time now to modernize the human diseasetaxonomy that informs healthcare decisions, by more
precisely defining and classifying diseases. A New
Taxonomy would precisely define diseases basedon their intrinsic biology, in addition to traditional“signs and symptoms”; and incorporate a deeper
understanding of disease mechanisms, pathogenesis,
and treatments in a dynamic, iterative fashion –continuously incorporating newly emerging disease
information. The New Taxonomy would requirean Information Commons in which data on large
populations of patients become broadly available forresearch use, and a Disease Knowledge Network that
adds value to these data by highlighting their interconnectedness and integrating them with evolving
knowledge of fundamental biological processes. The
result would be “Precision Medicine”.

October 22-23, 2012, the MSM Consortium decided
to focus its discussions on “Multiscale Modeling for
Precision Healthcare”3. Four sets of MSM Consortium
panelists were asked to read the NAS report on
Precision Medicine1, collect other sources of relevant
information, and comment on the following questions:
1) Data-driven models, physiological models andstructural models – can they be tailored to individuals
for precision healthcare? 2) Can computational
biology facilitate precision healthcare? 3) How can
we utilize clinical data to inspire MSM research for
precision healthcare? 4) How do we incorporate
verification, validation and uncertainty quantificationin MSM for precision healthcare? The panelists led
the MSM Consortium in an interactive discussion oneach of these questions, and that the views expressed
below are those of the MSM Consortium, not the NIH
or other government agency members of IMAG.
2. MSM and Precision Medicine
Though the NAS report on Precision Medicine1
does not refer to computational models, the MSM
Consortium concluded that multiscale models will
be necessary platforms to derive knowledge fromclinical and scientific data and to integrate knowledge
for the purpose of informing diagnoses, supporting
therapeutic decisions and predicting clinical outcomes.
At the same time the improved data resources, throughthe New Taxonomy, Information Commons, Disease
Knowledge Networks recommended in the NAS
report, will all be critical in enabling the applicationof multiscale models to precision medicine. A
recently published review article by Winslow et al.4
defines computational medicine as the means to use
quantitative multiscale models to understand alteredstructure and function in disease, and develop new
methods for disease diagnosis and treatment.

1.2. IMAG and the MSM Consortium
Since 2003, the Interagency Modeling and Analysis
Group (IMAG), a coalition of program officers from10 government agencies in the United States andCanada, has been promoting funding activities in
the area of modeling and analysis of biomedical,
biological and behavioral systems2, with a particularemphasis on multiscale modeling. Since 2006,
IMAG has facilitated the activities of the Multiscale
Modeling (MSM) consortium, which is made up
of investigators in the field. Each year, the IMAGMSM Consortium meets to discuss timely issues thatconcern the field of multiscale modeling. This year, on

The MSM Consortium also discussed the use ofcomputational (systems) biology models, and its
unique ability to link molecular scales with macroscale information. Though computational biology
models have already shown success in correlating
high-dimensional, high-throughput data with disease22

4th Biosupercomputing Symposium (December 3 - 5, 2012)

References

diagnosis, there is a need to increase the influenceof mechanistic hypotheses in the modeling cycle5.
The MSM Consortium also discussed the need forgood clinical data appropriate for models, as well asa need to better engage the clinical community in themodeling discussion. Finally, the MSM Consortium
conducted an extensive discussion on model
validation, verification and uncertainty quantification(VV&UQ). The Consortium identified the need todevelop model evaluation criteria at the beginning of
the model development process, and that VV&UQ
should be addressed in the context of a model’s use
case and level of criticality of the clinical decision
being made.

1. Committee on a Framework for Development a
New Taxonomy of Disease, National ResearchCouncil. “Toward Precision Medicine: Buildinga Knowledge Network for Biomedical Research
and a New Taxonomy of Disease”, The National
Academies Press (2011),
http://www.nap.edu/catalog.php?record_id=13284.
2. Peng, G. C.Y. “Editorial: What Biomedical
Engineers Can Do to Impact Multiscale
Modeling”, IEEE Transactions on Biomedical
Engineering, Volume 58, Issue 12, Part 2, Page(s):3440 – 3442 (2011)
3. 2012 Multiscale Modeling Consortium Meeting,
http://www.imagwiki.nibib.nih.gov/mediawiki/
index.php?title=2012_MSM_CONSORTIUM_
MEETING
4. Winslow, R. L., Trayanova, N., Geman, D., Miller,
M. I. “Computational Medicine: Translating
Models to Clinical Care”, Sci Transl Med, Volume
4, Issue 158 (2012)
5. An, G. “Closing the Scientific Loop: Bridging
Correlation and Causality in the Petaflop Age”,
Sci Transl Med, Volume 2 Issue 41 (2010)

In conclusion, multiscale modeling has the potential
to play a critical role in implementing the goals of
precision medicine. While most diseases present
clinically at tissue, organ and whole body scales; they
are frequently treated with molecular interventions.
Multiscale models are among the most rigorous tools
available to integrate a hierarchy of information at
different scales and predict the dynamic state of apatient and the progression of disease.

23

4th Biosupercomputing Symposium (December 3 - 5, 2012)24

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Multiscale simulations of tumor angiogenesis
with therapeutic applications
Aleksander S. Popel
Professor, Biomedical Engineering, Johns Hopkins University School of
Medicine

Profile:Aleksander S. Popel, Ph.D., is a Professor of Biomedical Engineering at the Johns Hopkins University Schoolof Medicine. He holds joint appointments as a Professor of Oncology in the School of Medicine, and Professorof Mechanical Engineering and Chemical & Biomolecular Engineering in the Johns Hopkins Whiting Schoolof Engineering. His research areas are systems biology, computational medicine & biology, angiogenesis andmicrocirculation. He published over 250 scientific papers in these areas. He served as a Visiting Professor at MIT
and Harvard University, and a Visiting Fellow at the Isaac Newton Institute, University of Cambridge, U.K. His
honors and awards include the Eugene M. Landis Award from the Microcirculatory Society, C. Forbes Dewey
Distinguished Lectureship in Biological Engineering at the Massachusetts Institute of Technology, A.C. Suhren
Lectureship at Tulane University, Robert M. and Mary Haythornthwaite Distinguished Lectureship at Temple
University, and Kawasaki Medical Society Lectureship in Japan. He delivered a keynote address for The Virtual
Physiological Human (VPH) European Union Physiome Project in 2010; he received the William H. Huggins
Excellence in Teaching Award from Johns Hopkins University. He is an elected Fellow of the American Instituteof Medical and Biological Engineering, American Heart Association, American Physiological Society, andAmerican Society of Mechanical Engineers, and an Inaugural Fellow of the Biomedical Engineering Society.
He has been a member of editorial boards of biological and biomedical engineering journals, such as Annals ofBiomedical Engineering, American Journal of Physiology, Microcirculation, Microvascular Research, Cellular& Molecular Bioengineering, Wiley Interdisciplinary Reviews on Systems Biology & Medicine, Frontiers in
Computational Physiology & Medicine, and American Journal of Cancer Research. He is Co-Chairman of the
Physiome and Systems Biology Committee of the International Union of Physiological Sciences and Co-Leader,
Multiscale Systems Biology Working Group for the Interagency Modeling and Analysis Group (including NIH).
He has served in an advisory role to biotech companies. He regularly serves on grant review boards and advisory
panels at the National Institutes of Health, National Science Foundation, and other US and international funding
agencies.
25

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Tumor growth and metastasis are complex dynamic
processes that so far have largely evaded attempts
to control and manage them; a grand challenge ofmodern biology and medicine is understanding theirmechanisms at the molecular, cellular and tissuelevels and learning how to control the disease [1].
Considering the enormity of the task and thecomplexity of the subject, systems and computational
biology are absolutely necessary for achieving thesegoals. To make a crucial progress in the understandingof cancer and the role of new vasculature that sustains
tumor growth and metastasis [2, 3], new approaches
are necessary that would allow us to modularize theproblems, and to identify and quantify the interactions
between the modules. We need to combine
reductionist approaches, down to the gene and protein
levels, with integrative or systems approaches. Whatshould make the problem feasible is its modularity,
ie the ability to define and study the modules
individually as a first approximation, then combine
them to study interactions between the modules, up
the hierarchy of biological scales and complexity.

including growth factor receptor-ligand
interactions (Vascular Endothelial Growth Factor– VEGF) [9-11], intracellular signaling [12], andenzymatic processes that involve cell-extracellularmatrix interactions specifically with matrix
metalloproteinases, MMP [13-15].
4. Mechanistic models of transcriptional regulation;
Hypoxia-Inducible Factor HIF1α [16].
5. Agent-based three-dimensional modeling ofangiogenesis and tumor growth [17, 18].
6. Multiscale image-based modeling of tumor bloodflow and molecular transport [19].
7. Molecular-based mechanistic models of therapeuticinterventions in cancer [20, 21].
8. Computational modularity and systems integration[22, 23].
It is important that every module be validated
individually, to the best extent possible, using
experimental methods, from molecular biology, to in
vitro cellular assays, to animal experiments. However,some of the pathways may interact (cross-talk) andeven exhibit synergistic behavior, e.g. [24]; in thiscase the hierarchical complexity of the interactions
needs to be built up gradually so as not to lose theability to validate at every step.
Although our studies and those of other laboratories
demonstrate growing technical ability to approach
the systems biology of cancer, admittedly they only
scratch the surface in our quest to quantitatively
understand cancer. Analyzing and quantifying themodularity of multiple signaling pathways in thetumor microenvironment, e.g. cancer cells, tumor
endothelial cells, and stromal cells, including
immune cells will be a challenging task. Importantand significant spatial heterogeneities exist within
individual tumors that are the hallmark of cancer;
these heterogeneities are manifested in genetic
variability, physico-chemical properties, and drug
response. These issues need to be confronted infuture work. We predict that computational biology
of cancer will remain a major research frontier; it will
require novel methodologies and is likely to stretch
the capabilities of modern supercomputers.

2. Unraveling the Complexity of Cancer by
Multiscale Multimodular Modeling of Tumorand Vascular Dynamics
What are the major processes involved and how
do we approach the problem? Below is a list of theprocesses where we have made initial progress.
To be specific, here we focus on breast cancer, but
the methodologies developed should apply to solid
tumors and their metastatic dissemination moregenerally. Computationally, these are presented as
modules describing different biological scales as wellas integration:
1. Bioinformatic analysis of the anti-angiogenic
motifs in the human proteome and computer-aided
optimization of anti-angiogenic peptides [4, 5].
2. Bioinformatic analysis of protein-protein
interactions (PPI) with the aim of classifying andorganizing proteins important for angiogenesis andidentifying novel drug targets [6-8].
3. Mechanistic models of cell signaling networks,
26

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Work supported by the National Institutes of Healthgrant R01 CA138264, The Safeway Foundation for
Breast Cancer, The Thome Foundation.

[16] A.A. Qutub, A.S. Popel, Molecular and cellularbiology, 28 (2008) 5106-5119.
[17] A.A. Qutub, F. Mac Gabhann, E.D. Karagiannis,
P. Vempati, A.S. Popel, IEEE engineering inmedicine and biology magazine, 28 (2009) 1431.
[18] A.A. Qutub, A.S. Popel, BMC systems biology, 3(2009) 13.
[19] E. Kim, S. Stamatelos, J. Cebulla, Z.M.
Bhujwalla, A.S. Popel, A.P. Pathak, Annals of
biomedical engineering, 40 (2012) 2425-2441.
[20] S.D. Finley, A.S. Popel, AAPS J, 14 (2012) 500509.
[21] M.O. Stefanini, F.T. Wu, F. Mac Gabhann, A.S.
Popel, Cancer research, 70 (2010) 9886-9894.
[22] G. Liu, A.A. Qutub, P. Vempati, F. Mac Gabhann,
A.S. Popel, Theoretical biology & medical
modelling, 8 (2011) 6.
[23] F. Mac Gabhann, M.O. Stefanini, A.S. Popel,
Simulating therapeutics using multiscale modelsof the VEGF receptor system in cancer., in: T.L.
Jackson (Ed.) Modeling Tumor Vasculature:
Molecular, Cellular, and Tissue Level Aspects
and Implications., Springer, 2012, pp. 37-53.
[24] J.E. Koskimaki, E. Lee, W. Chen, C.G.
Rivera, E.V. Rosca, N.B. Pandey, A.S. Popel,
Angiogenesis, (2012).

References[1] D. Hanahan, R.A. Weinberg, Cell, 144 (2011) 646674.
[2] P. Carmeliet, R.K. Jain, Nature, 473 (2011) 298307.
[3] R.K. Jain, P. Carmeliet, Cell, 149 (2012) 14081408 e1401.
[4] E.D. Karagiannis, A.S. Popel, Proceedings of
the National Academy of Sciences of the UnitedStates of America, 105 (2008) 13775-13780.
[5] C.G. Rivera, E.V. Rosca, N.B. Pandey, J.E.
Koskimaki, J.S. Bader, A.S. Popel, Journal of
medicinal chemistry, 54 (2011) 6492-6500.
[6] L.H. Chu, C.G. Rivera, A.S. Popel, J.S. Bader,
Physiological genomics, 44 (2012) 915-924.
[7] C.G. Rivera, L.-H. Chu, J.S. Bader, A.S. Popel,
Applications of Network Bioinformatics to Cancer
Angiogenesis., in: A.S. Azmi (Ed.) Systems
Biology in Cancer Research and Drug Discovery,
Springer, 2012, pp. 231-244.
[8] C.G. Rivera, S. Mellberg, L. Claesson-Welsh, J.S.
Bader, A.S. Popel, PloS one, 6 (2011) e24887.
[9] F. Mac Gabhann, A.S. Popel, Microcirculation, 15(2008) 715-738.
[10] F.T. Wu, M.O. Stefanini, F. Mac Gabhann, A.S.
Popel, Methods in enzymology, 467 (2009) 461497.
[11] F.T. Wu, M.O. Stefanini, F. Mac Gabhann, C.D.
Kontos, B.H. Annex, A.S. Popel, Journal ofcellular and molecular medicine, 14 (2010) 528552.
[12] W.H. Tan, A.S. Popel, F. Mac Gabhann, PloS
one, (2013).
[13] P. Vempati, E.D. Karagiannis, A.S. Popel, The
Journal of biological chemistry, 282 (2007)
37585-37596.
[14] P. Vempati, F. Mac Gabhann, A.S. Popel, PloS
one, 5 (2010) e11860.
[15] E.D. Karagiannis, A.S. Popel, Journal oftheoretical biology, 238 (2006) 124-145.
27

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Is eukaryotic cell biology ready forsupercomputing? A perspective from the
Virtual Cell project
Ion Moraru
Director, High Performance Computing Facility and Associate Professor,
Cell Biology, University of Connecticut School of Medicine

Profile:Dr. Moraru is the Director of the High Performance Computing Facility and an Associate Professor of Cell
Biology at the University of Connecticut School of Medicine. Prof. Moraru received an MD from the Bucharest
Institute for Medicine and Pharmacy in 1988 and a PhD in Cell Biology from Carol Davila University in 1992.
His experimental studies in signal transduction led him to mathematical modeling and computer simulations, andhe eventually dedicated himself full-time to computational biology. His current main interest is the developmentof methods to bridge the gap between quantitative, "bottom-up", spatially-resolved, detailed simulations ofcellular functions, and large-scale, pathway- and "omics"-derived systems biology models.

29

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Logic-derived modelling has been used to map
biological networks and to study arbitrary functional
interactions, and fine-grained kinetic modelling canaccurately predict the detailed behaviour of wellcharacterized molecular systems, however, at present,
neither approach comes close to unravelling thefull complexity of a cell. It is possible nowadays tocreate a whole-cell model of a simple prokaryotic
cell that accounts for all molecular components andtheir interactions, from electrolytes and metabolites
to proteins and ribosome assemblies. A huge effortof painstaking data mining is required to assemble
such models, but the mathematical representations are
straightforward and the computational cost is low: anentire cellular lifecycle can be simulated with as little
as 1 TFlop*hr.
In contrast, we are far from being able to simulate
eukaryotic cells at a similar level of detail: the tools,
data, and computational power are all still inadequate.
Non-linear dynamics of intracellular molecularinteractions (especially in signalling networks, but
often also in metabolic or gene regulatory networks)
cannot be ignored. Emergent properties are due not
only to network topology, but to the detailed kinetic
rate laws and quantitative parameters. Multiple
phosphosites can create a combinatorial complexityof regulatory actions, and generate difficulties inmapping functional states to measured observable
quantities. Compartments, scaffolds, and diffusion
create spatial inhomogeneities and microdomains,which have critical functionality in most eukaryotic
cells (often also in prokaryotes). And to top it alloff, there has been increasing evidence that partsof the cellular machinery employ fleeting, nonstoichiometric, pleiomorphic assemblies of molecules
to carry out vital processes.

integrates new microscope technologies for making
quantitative in vivo live cell measurements with newphysical formulations and computational tools thatwill produce spatially realistic quantitative models
of intracellular dynamics. The latter are being madeavailable for the use of researchers worldwide throughtheir integration into the public, web-accessible, VCell
software framework. The creation of models in VCell
can range from the simple, to evaluate hypotheses
or to interpret experimental data, to the complex,
where multi-layered models can be used to probe thepredicted behavior of large, highly non-linear systems.
VCell has been designed from the beginning towork with spatially-resolved models, and it has
a wide range of capabilities to specify 2D or 3D
geometries for compartments. Geometries can be
re-used from other VCell models (including sharedmodels from other users). New analytic shapes can becreated by defining analytic equations that describe
compartment distributions in space or by manually
defining compartments on a regular grid of points inspace (image painting). An important VCell feature
is a possibility to use geometries from external
resources. Rasterized 2D or 3D (Z-stacks) images
can be imported, as well as CAD specifications
such as stereolithography and constructive solid
geometry. Additionally, surface mesh files can beimported and sampled within a spatial domain to
define compartments. Users can build complex
models which specify compartmental topology andgeometry, molecular characteristics, and relevantinteraction parameters. VCell automatically converts
the biological description into a corresponding
mathematical system, deterministic (concentrations,ordinary and/or partial differential equations) or
stochastic (locations, process probabilitites). VCell
will then generate the code to perform and analyze
simulations of the system using numerical solvers.
The distributed version of VCell maintains a databaseof models and simulation data that can be private,shared with specific users, or public.
VCell has been continuously and rapidly growing
in features and complexity over the past several
years. Compartmental and spatially-resolved

2. The Virtual Cell Project
The Virtual Cell (VCell) is a unique computational
environment for modeling and simulation in cellbiology developed at the University of Connecticut’sCenter for Cell Analysis and Modeling, an NIHdesignated National Research Resource. The Center30

4th Biosupercomputing Symposium (December 3 - 5, 2012)

models of reaction, diffusion, membrane transport,
electrophysiology, advection/flow and events can be
simulated both deterministically and stochasticallywith a large choice of solvers, including hybrid
stochastic-deterministic solvers, and advanced
parameter scan and estimation tools are available.

– a promise brought by community standards formodel exchange that are currently being developed(such as SBML level 3 extensions). But will it becomputationally tractable to build and simulatecomprehensive quantitative whole-cell models of
eukaryotic cells? A good benchmark towards thatgoal (perhaps a ‘grand challenge’?) would be theability to predict the phenotype produced by random
mutations in the case of a simple multi-cellularorganism such as C. Elegans (much in the same way it
was recently done for the bacterium M. Genitalium ).
Based on known –omics data, we can estimate thata hypothetical simulation with currently available
algorithms of the development lifecycle (from egg to959-cell worm) would require at least 100 EFlop*hr.
This is not realistic right now, but not out of reach,either; it seems reasonable to expect that the overallrequired computational infrastructure (algorithms,
standards, software) would be available by the endof the decade. But will there be enough experimental
data to parameterize, constrain and validate suchmodels?

3. Bridging the Systems Biology Gap
New developments will be presented that expand
VCell’s capabilities towards the systems biology
domain, such as automated model building from
pathway databases, rule-based modeling, and networkfree simulations. In parallel, we will discuss recent
developments of systems biology modeling toolsthat expand their capabilities towards the detailed
biochemistry domain, such as automated differential
equations generation from Boolean networks andcombined data-driven and pathway-driven methods of
network inference.
We now need efficient ways for seamless interfacing
between these different tools in order to combine allof the required different mathematical formalisms

31

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Modeling Nanotoxicity: Large Scale
Molecular Simulations of Protein-Nanoparticle
Interactions with IBM BlueGene
Ruhong Zhou
Manager, Soft Matter Theory and Simulation, Computational Biology Center
IBM Thomas J. Watson Research Center, IBM

Profile:Ruhong Zhou is currently a Research Staff Scientist and Manager of Soft Matter Theory and Simulation
Group, Computational Biology Center, IBM Thomas J. Watson Research Center, and an Adjunct Professor at
Department of Chemistry, Columbia University. He received his Ph.D. with Prof. Bruce Berne in chemistry from
Columbia University in 1997. He joined IBM Research in 2000, after spending two and a half years working
with Prof. Richard Friesner (Columbia Univ) and Prof. William Jorgensen (Yale Univ) on polarizable force
fields. He has authored and co-authored 120 journal publications and 15 patents, delivered 150+ invited talks atmajor conferences and universities worldwide, and chaired and co-chaired many conferences in computational
biology, computational chemistry and biophysics fields. He is part of the IBM Blue Gene team who won the2009 National Medal on Technology and Innovation. He has won the IBM Outstanding Innovation Award in
2011, IBM Outstanding Technical Achievement Award (the highest technical award within IBM) in 2008 and2005, IBM Research Division Award in 2005, Columbia University Hammett Award in 1997 (for best graduates),and American Chemical Society DEC Award on computational chemistry in 1995. His current research interests
include development of novel methods and algorithms for computational biology and bioinformatics, and large
scale simulations for protein folding, ligand-receptor binding, protein-protein interaction, and protein nanoparticle
interactions. He was elected to AAAS Fellow (American Association of Advancement of Science) and APS
Fellow (American Physical Society) in 2011.
33

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Nanoscale particles have become promising materials
in various biomedical applications, such as cancer
therapeutics, diagnosis, neuro-imaging, drug delivery,
and biosensors. However, in order to stimulate andfacilitate these applications, there is an urgent need forthe understanding of the nanotoxicity and other risks
involved with these nanomaterials to human health.
Some of our recent molecular modeling works withIBM Blue Gene supercomputer will be presented inthis talk. We show that a pristine carbon nanotube, oneform of hydrophobic nanoparticles, can interact anddisrupt the structures and functions of many important
proteins such as WW domains, SH3 domains, andhuman blood serum proteins. In some extreme cases,
such as the WW domains, the carbon nanotube canunexpectedly plug into the hydrophobic core of theprotein to form stable complexes. This plugging of
nanotubes disrupts and blocks the active sites of WW
domains from binding to the corresponding ligands,
thus leading to the loss of protein functions. In othercases, nanotubes compete with ligands for the receptor
binding sites involved in the signaling and regulatory
pathways. Different adsorption capacities of humanserum proteins on carbon nanotubes, on the other
hand, result in different cytotoxicity. The hydrophobic
interactions between the carbon nanotubes andhydrophobic residues, particularly aromatic residues
through the so-called π−π stacking interactions, arefound to play key roles. In addition, the molecularmechanism of tumor metastases inhibition by
metallofullerenol Gd@C82(OH)22 (i.e. toxic tocancer cells) has been studied with both experimental
and theoretical approaches. These findings mightprovide a better understanding of “nanotoxicity” at
the molecular level and help design better therapies
with nanomedicine. In the following, two specificexamples are provided to illustrate these important
protein-nanoparticle interactions.

fluorescence spectroscopy, CD, SDS-PAGE, andmolecular dynamics simulations, to investigate theinteractions of SWCNTs with human serum proteins(BSA, γ-Ig, Tf, and BFg), and find a competitive
binding of these proteins with different adsorption
capacity and packing modes. The π–π stacking
interactions between SWCNTs and aromatic residues(Trp, Phe, Tyr) are found to play a critical role indetermining their adsorption capacity. Additional
cellular cytotoxicity assays reveal that the competitive
bindings of blood proteins on the SWCNT surface cangreatly alter their cellular interaction pathways andresult in different cytotoxicity for these protein-coated
SWCNTs. These findings have shed light towards thedesign of safe CNT nanomaterials by comprehensive
preconsideration of their interactions with humanserum proteins.

Figure 1. Molecular mechanism of pancreatic tumor
metastases inhibition by metallofullerenol Gd@
C82(OH)22
2.2 Pancreatic Tumor Metastases Inhibition
by Gd@C82(OH)22
Pancreatic adenocarcinoma is probably the mostlethal of the solid tumors and the fourth leading causeof cancer-related death in North America. Matrix

2.1 Competitive Binding of Carbon
Nanotubes to Human Serum ProteinsIn this example, we have used both experimental
and theoretical approaches, including AFM images,
34

4th Biosupercomputing Symposium (December 3 - 5, 2012)activity (see Figure 1). Both the suppression of MMPs
expression and specific binding mode make Gd@
C82(OH)22 a potentially more effective nanomedicine
for pancreatic cancer than the traditional medicines
which usually target the proteolytic sites directly,
but fail in selective inhibition. Our findings basedon a combination of in vivo, in vitro and in silico
approaches provide new insights for de novo designof nanomedicines for fatal diseases such as pancreaticcancer.

metalloproteinase (MMP) has long been targeted
as a potential anti-cancer therapy due to its seminalrole in angiogenesis and extracellular matrix (ECM)
degradation for the tumor survival and invasion. The
inhibition specificity to MMPs and the molecular level
understanding of inhibition mechanism, however,remains largely unresolved. Moreover, finding better
alternatives to traditional medicines with emerging
nanomaterials (nanomedicine) is of great current
interest.
In this example, we find that endohedral
metallofullerenol Gd@C82(OH)22 can successfully
inhibit the neoplastic activity of pancreatic cancer
with experiments at animal, tissue, and cellular levels.
Gd@C82(OH)22 effectively blocks tumor growthin human pancreatic cancer xenografts in nude micemodel. Enzyme activity assays also show Gd@
C82(OH)22 not only suppresses the expression of
MMPs but also significantly reduces their activities.
We then further applied large scale molecular
dynamics simulations to uncover the molecularmechanism by studying the Gd@C82(OH)22-MMP-9
interactions at atomic detail. Our data demonstrated
that Gd@C82(OH)22 inhibits MMP-9 mainly via
an exocite interaction while the well-known zinccatalytic site only plays a minimum role. Steered by
non-specific electrostatic, hydrophobic and specifichydrogen bonding interactions, Gd@C82(OH)22
exhibits specific binding modes near the ligand
specificity loop S1′, thereby inhibiting the MMP-9

References
1. C. Ge, J. F. Du, L. Zhao, L. Wang, Y. Liu, D. H.
Li, Y. Yang, R. H. Zhou, Y. L. Zhao, Z. F. Chai,
C. Y. Chen, Competitive binding of human serum
proteins on single-wall carbon nanotubes reduces
cytotoxicity, Proc. Natl. Acad. Sci, 108, 16968,
2011
2. G. Zuo, Q. Huang, G. Wei, R. H. Zhou, and H.
Fang, Plugging Into Proteins: Poisoning Protein
Function by a Hydrophobic Nanoparticle, ACS
Nano, 4, 7508-7514, 2010
3. S. G. Kang, G. Zhou, P. Yang, Y. Liu, B. Sun, T.
Huynh, H. Meng, L. Zhao, G. Xing, C. Y. Chen,
Y. L. Zhao, R. H. Zhou, Molecular Mechanismof Pancreatic Tumor Metastases Inhibition by
Metallofullerenol Gd@C82(OH)22 and its
Implication for de novo Design of Nanomedicine,Proc. Natl. Acad. Sci., 109, 15431, 201235

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Multi-scale modelling of the motor system
underlying goal-directed behaviour in avertebrate model organism
Sten Grillner
The Nobel Institute for Neurophysiology, Dept of Neuroscience,
Karolinska Institutet

Profile:Sten Grillner has unravelled the intrinsic function of the modular network organization underlying fundamental
aspects of our motor repertoire. His initial work defined the basic organisation of the mammalian locomotor
system in terms of supraspinal command systems, spinal networks coordinating the movements (CPGs), and thesensory control of the CPGs. To address the next level question - the molecular, cellular and synaptic design ofthese neuronal circuits - he developed a novel and simpler vertebrate model (lamprey). The different network
interneurons, their synaptic interaction (transmitters, receptor subtypes), and their membrane properties (ionchannel subtypes expressed) have been identified. The palette of different subtypes of ion channels expressedin different neurones is found to be of critical importance for network function. Through an interaction between
detailed multi-faceted experimentation and large scale modelling with biophysically realistic numbers of HodgkinHuxley neurons, the operation of this entire motor control system has been uncovered. This conserved system isunderstood at a cellular/molecular level, and the conceptual gap between gene/molecule and behaviour has beenbridged. We now understand the design of the neuronal subsystems coordinating goal-directed ambulation. In
focus are now the forebrain systems that are responsible for selection of different patterns of behavior.
Sten Grillner has also helped develop the OECD initiated International Neuroinformatics Coordinating Facility(INCF) with the secretariat in Stockholm that promotes both the development of interoperable databases andmulti-scale modelling. INCF has now 17 member countries extending from Japan, Korea and India to Europe andthe US.
37

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Multi-scale modelling of the motor system underlying goal-directed behaviour in avertebrate model organism
Sten Grillner, Nobel Institute for Neurophysiology, Dept of Neuroscience, Karolinska Institutet, SE-171
77 Stockholm
The vertebrate brain controls a great variety of
movements through dedicated networks like thosecontrolling eye movements, expression of emotions,
respirations and locomotion. These networks are to alarge degree conserved through the vertebrate phylum.
The neural mechanisms underlying the control
of goal directed movement will be in focus, bothbiological background information and multiscalemodelling from the subcellular and synaptic levelto the microcircuit, systems and behavioural levels.
The propulsive locomotor synergy is controlled fromcommand regions in mesencephalon, which in turncontrol central pattern generating (CPG) networks
in the spinal cord. My presentation, based on thelamprey CNS, will address the intrinsic functionof the adaptable CPG that can generate different
motor patterns. In addition, I will discuss the tectal
mechanisms underlying steering movements withthe retinotopic motor map, and the control fromthe output of the lamprey basal ganglia, which is of
critical importance for deciding which motor program
should be turned on at a given moment of time. The
mechanism by which different motor programs are
selected will be considered with special reference
to the basal ganglia – experiments and modelling.
Our recent findings establish that the structure andfunction of the basal ganglia have been conserved toa surprising degree throughout vertebrate phylogeny
over a period of more than 500 million years from theancient lamprey version to primates. This applies tothe input to striatum (pallium, thalamus, dopamine,
5-HT, histamine input), the pallidal structures(GPi, substantia nigra reticulata (SNr), GPe and thesubthalamic nucleus) and their output targets, thecellular properties of striatal and pallidal neurons
and the effects of an MPTP induced dopamine
denervation.

simulated these networks with an approximate numberof model neurons (around 15 000) corresponding
to that observed biologically. The model neurons
are of the detailed compartmental Hodgkin-Huxley
type with detailed membrane properties including
all subtypes of ion channels observed biologicallyand the different transmitter receptors. Each modelneuron is designed to be very close to its biologicalcounterpart including the variability observed in each
population. The synaptic interaction is modelled
as well as network properties. The variability inneuronal properties is important in that it provides
for stability in network properties over a large rangeof frequencies. These networks have also been used
to control a simulated body with the ability to movein a natural fashion and be able to perform orienting
movements, avoid obstacles and approach certain
goals. These simulations have been performed on an
IBM blue gene supercomputer.
Kozlov A, Huss M. Lansner A, Hellgren Kotaleski J,
Grillner S. (2009) Simple cellular and network control
principles govern complex patterns of motor behavior.
PNAS 24;106:20027-20032.
Grillner, S. (2006) Biological Pattern Generation: The
Cellular and Computational Logic of Networks in
Motion. Neuron 52; 751-766
Stephenson-Jones M, E. Samuelsson, J. Ericsson,
B. Robertson and S. Grillner (2011) Evolutionary
conservation of the basal ganglia as a common
vertebrate mechanism for action selection. Curr. Biol.,
21;1081-1091.
Stephenson-Jones M, O. Floros, B. Robertson, S
Grillner (2012) Evolutionary conservation of thehabenular nuclei and their circuitry controlling the
dopamine and 5-hydroxytryptophan (5-HT) systems.
Proc Natl Acad Sci U S A.

In order to simulate the basal ganglia – brainstem –spinal cord networks underlying different motor
patterns including visuomotor control, we have38

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Large-scale simulations of biomolecularmachines
Karissa Sanbonmatsu
Principal Investigator, Los Alamos National Laboratory

Profile:Dr. Sanbonmatsu has been a principal investigator at Los Alamos National Laboratory since 2001. She received
her BA in Physics from Columbia University in 1992 and PhD in Astrophysical, Planetary and Atmospheric
Sciences from University of Colorado at Boulder in 1997. Dr. Sanbonmatsu's research is focused the mechanismof non-coding RNAs, including the ribosome, riboswitches and long non-coding RNAs. She has pioneered largescale biomolecular simulations of nano-scale molecular machines such as the ribosome. She uses a combinationof high performance computing and experimental biochemistry to construct an integrated picture of large-scaleconformational changes of these complexes.

39

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Large-scale biomolecular complexes play a keyrole in many aspects of cellular activity, includingtranscription, translation, splicing, scaffolding,transport and metabolism. New supercomputing
technology has made it possible to study the dynamics
of these complexes in atomic detail using largescale molecular dynamics simulations. Our grouphas focused on a molecular machine called theribosome for the past decade.1 The ribosome playsa central role in all life forms and is responsible
for protein synthesis. The ribosome must read thegenetic information encoded on messenger RNA
and implement these instructions by producing
corresponding proteins. The ribosome performs the
only non-trivial information processing operationin cell by executing a 'look-up' table operation,
converting nucleic acid sequence into proteinamino acid sequence. In this sense, the ribosome isanalogous to the CPU of the cell and constitutes anano-scale information processor. Approximately
50% of the antibiotic drugs used in US hospitals
function by targeting the ribosome. Mechanistic
understanding of ribosome activity is important forthe design of new antibiotics that combat the growing
antibiotic resistant 'superbug' problem present intoday's hospitals. In addition, understanding how theribosome decodes genetic information will help lay
the foundation for the design of synthetic nano-scale
information processors.
1.2. The ribosome 'decoding problem'
By performing large-scale molecular dynamics
simulations of the ribosome, we have been able
to examine the inner workings of this molecularmachine. Our long-term goal has been understandingtranslocation, a large-scale massive rearrangementof the ribosome. As a first step, we focused on the
more tractable problem of decoding, the process bywhich the ribosome decodes genetic information (alsoknown as tRNA selection). Here, a key rearrangementthat occurs is called 'accommodation', where transfer
RNAs (tRNAs) carrying protein building blocks(amino acids) move into the ribosome. In our studyof this process, we identified a new functional regionof the ribosome ('the accommodation corridor')and predicted that certain parts of this corridor areimportant for ribosome function.2 Our predictions
have been validated in studies by several experimentalgroups.3-5 In an additional separate set of studies thatcombined explicit solvent simulations, reduced model
simulations, and single molecule experiments, a new
picture of ribosome function has emerged.6 Rather
than the ribosome machine parts moving in lock-step,
both simulations and single molecule experiments
show the tRNAs making large-scale reversible
excursions in a trial-and-error fashion. This picture isconsistent with a dynamic energy landscape view of
the ribosome.7
2. Translocation of tRNAs through theribosome
After studying the relatively tractable problem of'accommodation', we are now applying what wehave learned to the mechanism of translocation, amore complicated rearrangement where several large
conformational changes involving the entire ribosome
and tRNA occur. This motion allows the ribosome tomove exactly 3 nucleotides along the messenger RNA
to the next amino acid codon. Using microsecond
sampling in explicit solvent for the full ribosome, in
combination with experimentally measured rates,8
we estimate barrier heights for various motions
important for translocation. We have also used coarsegrained methods to simulate the various sub-steps oftranslocation.9,10 Our goal has been to use simulationsof the ribosome to produce detailed energy landscapesof translocation.
References
1. Sanbonmatsu KY. Curr Opin Struct Biol.
2012;22(2):168-74.
2. Sanbonmatsu KY, Joseph S, Tung CS.
Proc Natl Acad Sci U S A. 2005;102(44):15854-9.
3. Rakauskaite R, Dinman JD. RNA 2011,17:855864.41.
4. Burakovsky DE, et al. FEBS Lett 2011,585:30733078.
5. Burakovsky DE, et al. RNA 2010,
16:1848-1853.
40

4th Biosupercomputing Symposium (December 3 - 5, 2012)

6. Whitford, et al. RNA 2010, RNA. 2010;16(6):1196204.
7. Munro, et al. Trends Biochem Sci. 2009;34(8):390400.
8. Whitford, et al. J Am Chem Soc.
2010;132(38):13170-1.
9. Ratje, et al. Nature. 2010;468(7324):713-6.
10. Whitford, et al. PNAS. 2011;108(47):18943-8.

41

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Use of Modeling and Simulation in Drug
Discovery and Development: It is all about
the Questions
Sandra Allerheiligen
Vice President of Modeling and Simulation, Merck Research Laboratories

Profile:Sandra (Sandy) Allerheiligen, PhD is currently Vice President of Modeling and Simulation at Merck ResearchLaboratories. Prior to joining Merck, she held a variety of positions over her almost 20 years at Eli Lilly andCompany, including Global Sr. Director of Pharmacokinetics, Pharmacodynamics (PK/PD) and Trial Simulation,
Sr. Director of Drug Disposition, and most recently, Distinguished Fellow and Chief Scientific Officer of
Quantitative Pharmacology. Her research interests focus on study design and the application of mathematicalmethods to enable quantitative decisions for nonclinical and clinical development. She has applied PK/PD
modeling to oncolytic and endocrine agents. Her recent work is on the integration of biomarkers, PK/PD
modeling and trial simulation in non-clinical and clinical drug development, drug disease models and utilizationof quantitative and systems pharmacology approaches. Dr. Allerheiligen received a doctorate in pharmaceutics(with a specialization in PK/PD) from the University of Texas in Austin in 1985 and completed postdoctoral
fellowships at the University of Texas Health Center in San Antonio in 1986 and was a clinical assistant professorof Clinical Pharmacology from 1986 through 1990. Through her involvement in the American Associationof Pharmaceutical Scientists (AAPS) and other organizations, she has worked to expand the use of PK/PD
modeling and Quantitative Pharmacology methodologies in academia and across the industry. She is a Fellowof the American Association of Pharmaceutical Sciences (AAPS) and has served as Chair PPDM section, AAPS
Member-at-Large, and AAPS Treasurer. She co-founded the Population PK/PD Focus Group at AAPS andthe Midwest User’s Forum for Population Approaches in Data Analyses. She has been a member of Editorial
Advisory Boards for PharmSci, the AAPS Journal and the Journal of Pharmaceutical Research. She is an AdjunctFaculty member of the University of Osaka Medical School Department of Biostatistics and frequent lecturer onmodeling and simulation topics.
43

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
Impact of Quantitative and Systems Pharmacology (QSP) in Drug Discovery and Development: It is
all about the Question
Sandra R. B. Allerheiligen, PhD
The pharmaceutical industry must reduce costs while
delivering innovative medicines with improved
benefits for patients. With increasing computational
capability, scientists can leverage quantitative tools toanswer critical questions that influence the discovery
and development of important new medications.
These new approaches and tools have led to the
expansion of a specialized field, quantitative andsystems pharmacology (QSP, sometimes called
modeling and simulation).
Quantitative and Systems Pharmacology is an
integrative science incorporating relationships
between diseases, drug characteristics, and individual

variability to leverage existing knowledge and guide
future research. QSP helps define the questions andassumptions and highlights knowledge gaps. Throughthis approach it is possible to estimate drug efficacy,
understand safety, plan experiments, and informdiscovery/development strategies. Examples thathighlight the ability to answer such questions as 1)should this molecule be developed as a drug, 2) whatis the right dose to balance efficacy and safety, and
3) does this new drug increase bone strength will bepresented. The ultimate goal is not a mathematicalmodel but rather the ability to gain insight to optimize
new medications for patients.

44

ISLiM Talks

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Development of Multiscale Thrombosis
Simulator
Shu Takagi
Team Leader, Organ and Body scale Team, CSRP, RIKEN
Professor, Department of Mechanical Engineering, The University of Tokyo

Profile:Professional Position:
2010.4-current: Professor, Department of Mechanical Engineering, The University of Tokyo
2007.4- current: Team Leader, Organ and Body scale Team, CSRP, RIKEN
2002.1-2010.3: Associate Professor, Department of Mechanical Engineering, The University of Tokyo
1998.4-2002.1: Lecturer, Department of Mechanical Engineering, The University of Tokyo
1996.10-1998.3: Research Associate, Department of Mechano-Aerospace Engineering, Tokyo Institute ofTechnology
1995.4-1996.10: Research Associate, Department of Mechanical Engineering, The University of Tokyo
Education:
1995: Doctor of Engineering, Department of Mechanical Engineering, The University of Tokyo
1990: Bachelor of Engineering, Department of Mechanical Engineering, The University of Tokyo

47

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Thrombosis is regarded as one of the most importantdiseases, which cause the myocardial and cerebral
infarctions. It is very complicated disease affectedfrom molecular scale protein-protein interaction tocontinuum scale in blood flow. Initially, platelets
start to aggregate at the injured wall, where vonWillebrand Factor (vWF) molecules are attached. The
Glycoprotein, GPIb-α, on platelet membrane starts
showing ligand-receptor type interaction with thisvWF and platelets start aggregating around this spot.
From this stage, very complicated activated process ofplatelets and interactions with blood, vessel walls red
blood cells, fibrin etc. occur and they end up with theblockage of the vessels.
In the present study, we are developing the numericalmodel of the initial stage of thrombus
formation. To analyze this process, we have developed
the full Eulerian fluid-membrane coupling method [1],which is a further extension of the full Eulerian fluidstructure coupling method [2]. This method is coupledwith the stochastic Monte Carlo simulation on GPIb-αand vWF interactions. The basic concept of thismultiscale thrombosis simulator is shown in Fig.1.

vessels with large deformation constructing so-called
blood plasma layer, and that small platelets tends toflow in this plasma layer. It was observed that, oncethe platelets are coming closer to the vessel wall, they
stay there and it is rarely observed that they go backto the center of the vessel.

(a) t = 0.0 [s](b) t = 75.0 [ms]Fig.2 Time evolution of RBCs shape and location
Next, the effect of RBCs on the adhesion of platelets
is analyzed. The simulation was conducted throughthe coupling of continuum mechanical simulation
with Monte Carlo simulation of protein(GPIb-α) protein (vWF) interactions. The result shown in Fig.3
illustrates the adhesion of platelets by the binding
effect of proteins is achieved. This phenomenon is notobserved without Red blood cells. More details will
be discussed in the presentation.

Blood flow(continuum mechanics)

Ligand(wWF)-Receptor(GP1bα) binding(statistical mechanics)•Substance diffusion•Metabolic reaction and activation•Morphology change

Molecular interaction(molecular dynamics)

Fig.1 Multiscale modeling on the platelets adhesion

Fig.3 Platelet adhesion on Vessel Wall

2. Simulation Results
The pressure driven flows including the multiple
RBCs and platelets in a circular tube were carriedout. Fig.2 illustrates the time evolution of the shapeof RBCs and platelets. It is well-known that, in actual
blood flows, RBCs tends to flow in the center ofReferences
1. Ii, S. et al., Commun. Comput. Phys., 12, pp.
544-676. (2012).
2. Sugiyama, K. et al, J. Comput. Phys., 230, pp.
596-627 (2011)

48

49

Previous studies for platelet adhesion: （Small Vessel Diameter)・Tsubota et al.(2006) (Particle Method)
・Mori et al. (2008) (Stokesian Dynamics & Voigt Model),・Mody & King (2008) (BEM & Monte Carlo),
・Fogelson & Guy (2008) (Multiscale IBM)

Fixed-grid, Eulerian FSI-FD method coupled with Stochastic
Monte Carlo Method for recepter-ligand bindings

 dv O 10 3
dp

Require massively parallel computation
 development of suitable numerical methods

Coronary artery diameter (dv): 1-5 mm,
Size of Platelet (dp): 1-3 mSimulate thrombosis causing a heart attack

Objective

4th Biosupercomputing Symposium

Shu TAKAGIRIKEN,The University of Tokyo

Development of Multiscale
Thrombosis Simulator

Molecular interaction(molecular dynamics)

Solidby left Cauchy-Green
deformation tensor

・ Sugiyama, Ii et al. (2011) J. Comput . Phys., 230, 596.
・ Ii, Sugiyama et al. (2011) Int. J. Numer. Meth. Fluids, 65, 150.

by the displacement ofmaterial points themselves

1
0.5
0

s
by solid volume fraction
How is the solid deformation described?

by the boundary of mesh

Solid

vs.
Lagrangian
Eulerian
How is the two-phase distinguished ?
Fluid
Fluid

Full Eulerian approach for Fluid-Structure Interactions•Substance diffusion•Metabolic reaction and activation•Morphology change(statistical mechanics)

Ligand(wWF)-Receptor(GP1b) binding

Blood flow(continuum mechanics)

Multiscale modeling of initial stage of thormobosis

4th Biosupercomputing Symposium (December 3 - 5, 2012)solidvolumefraction s Bij for  s   min
B ij  
for  s   min
 0
1/2
=
B ij s  ij at t  0,

B  F  FT ,
 xi
Fij 
,
X j

left CauchyGreen
deformation
tensor50

Software availablle at http://www.islim.org/islim-dl_e.html

ZZ-EFSI achieved the actual speed of 4.5 PETA FLOPS!!!

Performance of the simulator (Program Name: ZZ-EFSI)

St. Venant-Kirchhoff model:

t=5.86t=1.17t=7.03t=2.34t=8.20t=3.52ー Zhao et al. (Lagrangian for solid)
 Present (full Eulerian) NxxNy=512x512* Ii et al., Commun. Comput. Phys., accepted (2011).

1  J s  1,  2 tr( B s )1
2Jsm Eb  κ   R P q  ( s  m )  P

bending model (Pozrikidis **)**Pozrikidis, J. Fluid. Mech., 440 (2001) 269.
Membrane stress： τ s  q s n v    0t
B s
B s   s v   s vT  B s
 v   Bst|  |
J s
 v   J s J s  s  v
t
 R
0
 v   R *Evans and Skalak, CRC press, Boca Raton, FL. (1980).
t



 v

 v   v   p    (  ( v   v T ))  |   |  s  ( τ s  q s n ) t  v 0

Extension to Eulerian fluid-membrane coupling model*

t=4.69t=0

s=f=Ls=Lf=1, =102, G=0.1

Comparison with the previous study
by Zhao, Freund & Moser (JCP 2008)

In-plane stress (Evans-Skalak*)Cs 1 
τs
B s   E s 1  C s 2P1  1 ( 1  1) 2 s σ s G ( s1 2 B  B   s1 B ).
(Linear) Mooney-Rivlin model:  s σ s G s1 B   2 c2 s1 2 tr( B ) B  B  B  .

Hyperelastic solid stress:

strain rate

1   vi  v j 


Dij

,
2  x j  xi  f ,ij  2 Dij ,

Ds
 0,
Dt
D B ij  vi vBkj  j B ki = 0,
Dt
xk
xk

 v
 vi
v 
p   m ,ij 0,  m  i  v j i   fi ,

 xi
x j 
 xi
x j
 t
 m ,ij  (1   s ) f ,ij   s s,ij ,Cauchy 's stress tensor

Basic Equations

4th Biosupercomputing Symposium (December 3 - 5, 2012)

51

Time history of platelets position（Platelets move along the wall,
once they are coming out near
the wall.）

Red blood cells:- Flowing in the core region of the channel- Inducing the velocity fluctuations in thewall-normal directions
Platelet:- Adhering on the wall due to the velocity
fluctuations caused by RBCs.

with RBCs

Results of the multiscale simulation about
the platelet adhesion on vWF-immersed wall

Simulation results

without RBCs
(no Adhesion simulated)

Confidential

Gaehtgens et al.(1980) Blood Cells., 6, 799.

Experimental Observations
・Cell free layer near vessel wall
・RBCs with slipper shape

Simulation results for the flow wit RBCs

*S.-Z. Luo, et al., Blood, 109 (2007) 603.

vWF binding site

Summary(Kim et al. 2011, Fox et al. 1988, Aryaet al. 2005)

・ Sugiyama, el al, J. Comput. Phys., Vol.230 (2010), pp.596-627.
・ Takagi et al. J. Appl. Mech. Vol.79 (2011), 010911.
・ Ii et al., Int. J. Numerical Method. in Fluids, Vol.65 (2011). pp.43-66.
・ Ii et al., Comm. in Cmuput. Phys, (2012), to appear.
・ Shiozaki et al. JBSE, Vol.7 (2012), pp.275-283.
・ Ii et al., J.Comput Phys, , Vol.231 (2012), pp.2328-2358.

References:Preparation for the comparison with in-vitro experimental data
Scale-up of a vessel diameter toward D=1 [mm] with activation of many
chemical reaction toward heart attack simulation.

On going study Initial stage of thrombosis is modeled. Platelet adhesion by the binding of
GP1ba and vWF molecules was simulated, through the coupling of thecontinuum scale simulation with the stochastic Monte Carlo simulation ofmolecular bindings.

 A Novel Method was developed for Fluid-Solid and Fluid-Membrane Interaction
Problems, which is suitable for the massively parallel computing.

f   p ( l  l0 )

l0 60


[nm], k r 0 3 [s  1 ]4
 p 10[N/m],  ts 0.9 p [N/m]( R f , Rr  [0,1] : Random numbers) Pr 
1  exp(  k r  t )  Rr  dissociationPf 
1  exp(  k f  t )  R f  binding( l  l0 ) 2 k r (l ) k r 0 ex p  ( p   ts )
2 k bT ( l  l0 ) 2 k f (l ) k f 0 exp   ts

2 k bT Stochastic model with energetic elasticity(Bell 1978, Dembo et al. 1988)

Ligand-Receptor bond model( Use of Transition State Theory)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Multi-scale, multi-physics heart simulator"UT-Heart" for heart research
Seiryo Sugiura, MD, PhD, FAHA
Professor, Graduate School of Frontier Sciences, University of Tokyo

Profile:Education:
Institution				DegreeYear
Field
University of Tokyo, School of Engineering B.S.
1975
Chemical Engineering
University of Tokyo, School of Medicine
M.D.
1982
Medicine
Medical Training and Employment:
1982-1984 Resident in Internal Medicine, University of Tokyo Hospital
1984
Staff Physician, University of Tokyo Hospital
1985-1987 Research Fellow, Department of Biomedical Engineering, Johns Hopkins University, Baltimore, MD
1987-1989 Staff Cardiologist, JR Tokyo General Hospital, Tokyo
1989-1990 Staff Cardiologist, University of Tokyo Hospital
1995-2002 Assistant Professor of Medicine, University of Tokyo
2002-present Professor, Graduate School of Frontier Sciences, University of Tokyo

53

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Cardiology is one of the earliest fields of medical
research where numerical simulation was recognized
as a powerful tool for promoting our understanding on
the disease processes. The fields of application coverelectrophysiology, cell metabolism, cardiac mechanics
and hemodynamics, but, so far, most of simulation
models have only dealt with either microscopic ormacroscopic aspect of the problem.
Today, we are facing the explosion of knowledge
gained by the progress in molecular and cellularbiology. To fully make use of these pieces ofknowledge for the solution of various health problems,
multi-scale modeling integrating the hierarchicalbiological system is required. To achieve this goalin heart research, we have developed a multi-scale,multi-physics heart simulator “UT-Heart”.
2. UT-Heart for clinical applications
Our heart simulator is based on the finite element
method and its morphology was created from themulti-detector CT or MRI data. In each element of thesimulator, mathematical models of cardiac excitationcontraction coupling process are implemented sothat each element behaves as a virtual myocyte. The
excitation starts at the pacemaker site propagates in
the heart tissue to trigger the synchronous contraction
and relaxation of the heart. Because the blood in theheart chamber is also modeled as the fluid elements
and its interaction with the cardiac tissue is solvedby the strong coupling method, we can also analyze
the hemodynamic parameters. In addition, wholecoronary arteries, capillaries and veins are modeled toshow the flow dynamics influenced by the myocardial
contraction at each level.
All these features of the model open the possibility ofin silico diagnosis and treatment of the virtual heart.
We have started a feasibility study by developing thepatient-specific heart and torso models. With theseheart models, we have successfully reproduced theactual ECG and UCG data of each patient. Predictionof therapeutic effect and optimization of the protocol

for the cardiac resynchronization therapy is anothertarget of the research. UT-Heart can also be applied
to the development of novel diagnostic or therapeutic
devices.
3. UT-Heart for basic science
The power of K-computer enables us to develop
a further detailed model in which seamlessintegration of the heart activities from microscopic tomacroscopic scales are realized. In this model, cycling
crossbridges of sarcomere in each cardiac myocyte
are modeled with its microscopic structures so thatthe pumping function of the heart is influenced by theforce generating activity of each myosin molecule.
Heart model with such detailed subcellular structure
can be used as the genetically modified human heartmodel thus will surely contribute to the basic science.
Researchers & Collaborators
Graduate School of Frontier Sciences, Univ. of Tokyo
Takumi Washio, Junichi Okada, Hiroshi Watanabe,Toshiaki Hisada
Dept. of Cardiovascular Medicine, Univ of Tokyo
Hiroshi Yamashita, Taro Kariya, Yasushi Imai, Ryozo
Nagai
Fujitsu Ltd.
Yoshimasa Kadooka, Akira Hosoi, Masahiro
Watanabe, Takao Hirahara, Takashi Yamazaki, Takashi
Iwamura, Machiko Nakagawa, Kohei Hatanaka,
Kazunori Yoneda
Acknowledgement
This research is supported by the Japan Society for
the Promotion of Science (JSPS) through its “Funding
Program for World-Leading Innovative R&D on
Science and Technology (FIRST Program)."
References
1. T Washio, J Okada, T Hisada SIAM J. SCI.
COMPUT. 30:2855–2881 (2008)
2. S Sugiura, T Washio, A Hatano, J Okada, H
Watanabe, T Hisada Prog. Biophys. Mol. Biol. inpress (2012)

54

55

EP modelCa

2+

IK1 IK INaK INa

JSR NSR……….. ICa,bICa INaCa Ica

INa,b

……..

Iion=INa+IK1+Ito+IKr+IKs+ICaL+INaCa+INaK+IpCa+IpK+IbCa+IbNa

TenTusscher K.H.W.J. et al Am J Physiol 2003
Courtemanche M. et al Am J Physiol 1998

Metabolism

V
Im

I ion

Monodomain model

Cm

Im

Bidomain model

Ising model by Rice J.J. et al. Biophys J 2003

Sarcomere model

Virtual tissue

From molecules to cell and tissue

Virtual myocyte

3Electricalactivity

Biochemistrymechanics

Channels

ATP

10-9m

K+

10-2m

Cell

In silico heart

100m

後壁側

Rice J.J. et al. Biophys J
2003
modified to ODE
ファイバー方向
under submission

Contraction model

Ventricle:
TenTusscher K.H.W.J. et al Am J Physiol 2003
Atrium:Courtemanche M. et al Am J Physiol 1998
Conduction system:
Stewart P et al. Philos Transact A 2009

EP modelHuman body4

Basic science
Pt. care

Organ

ECG

P-V loop

Action potential

Na+ Ca++Multi-scale

Molecularfunction/structurelactateglycolysis

GlucoseX-bridge

Tension

Cell models and boundary conditionsof knowledge

Integration

Multi-physics

4th Biosupercomputing Symposium

Seiryo Sugiura
Graduate School of Frontier Sciences
The University of Tokyo

Multi-scale, multi-physics heart simulator "UTHeart" for heart research

What is multi-scale, multiphysics simulator?

42

4th Biosupercomputing Symposium (December 3 - 5, 2012)CT data

I
新心電図

II

Ⅱ lead

V1

V2

V3

V4

56
BNP

CK-MB

Images troponin-THbA1c

Labo. data

V5

V6

Membranepotential

APD distribution
Transmural or apicobasal

Tailor-made simulation

III

ECG

Body surface map

Okada J et al. Am J Physiol 2011

In silico diagnosis

5

7

5

UCG

4th Biosupercomputing Symposium

1. Prediction of outcome of the treatment
2. Optimizing the treatment

Purposes of patient-specific simulation

In silico diagnosis

8

6

4th Biosupercomputing Symposium (December 3 - 5, 2012)

57

4th Biosupercomputing Symposium

Multi-scale simulation

4th Biosupercomputing Symposium

Seamless integration of the biological system

11

9

Na

Lumped modelCa

A.D.A.M Interactive physiology

Reproducing subcellular structure（Extended）Multi-scale model： seamless integration

Representavie Na,
Ca, ..channels
Contractile proteins
mitochodria

Standard model

DW Fawcett 1965

Excitation
Membranepotential
Contraction
Force

The power of “K-computer” extends the limit

10

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Optimization of Life-Science applications on
the K computer
Yousuke Ohno
Senior Scientist, Computational Science Research Program, RIKENProfile:March, 1996 Ph.D. in Science, Department of Astronomy, Graduate School of Science, The University of Tokyo
April, 1996 - March, 1999 Special Postdoctoral Researcher, RIKEN
April, 1999 - March, 2001 Contract Researcher, Computational Science Laboratory, RIKEN
April, 2001 - March, 2004 Research Scientist, Integrated Volume CAD System Research Program, RIKEN
April, 2004 - March, 2008 Research Scientist, Genomic Sciences Center (GSC), RIKEN
April, 2008 - current Senior Scientist, Computational Science Research Program, RIKEN
April, 2012 - current Senior Scientist, Quantitative Biology Center (QBiC), RIKEN

59

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction“Integrated Simulation of Living Matter (ISLiM)
Group" is developing "Grand Challenge Applications"in the life sciences that categorized five scales/layers
using the K computer and future super computers.
Massive parallelization is trend of high-performancecomputing. In the case of the K computer, the system
has 82,944 processors, 8 cores per processor, two
2-way SIMD floating-point multiply-and-add units percore. Therefore we needed to optimize applications
by hybrid parallel, 80,000 MPI process, 8 threads, 2
SIMD. To achieve such high-parallelization, we must
optimize application at deep layer such as structure ofcode and data, algorithm, etc.
High-performance Computing Team is assigned
supporting optimization of applications that aredeveloping other five teams and developing coreapplications/libraries for life science on the K
computer.
In each category, one or more applications achieved
parallelization of 24,576 processes which was
maximum number of process in early access to the
K computer. Three applications achieved full system
of K computer at special trial. High-performance
Computing Team is developing the molecular
dynamics core program, one of these applications.
We performed MD simulation of protein in waterwith 522 million atoms and achieved 4.4 PFLOPS of
calculation performance using 79,872 processors.
2. Molecular Dynamics Core Program
We developed the molecular dynamics core program
aimed for establishment of optimization technique,providing optimized code on the K computer.
2.1. Optimization on the K computer
In the MD code, the calculation of pair-wise force
and potential is most significant part of calculation.
The optimization of most inner loop of the force andpotential (kernel loop) is hardly depend processor andcompiler. On the K computer, SIMD and software
pipelining are most important optimization andsensitive to structure of data and loop. We changed
the data used in kernel loop from STL vector tosimple C style array because the compiler can’t

SIMD-optimize to calculation using STL vector. We
also rewrote “if” statement as mask operation because
the processor of K computer has SIMD instructionof masked operation for “if” statement. Applying
such “K” specific optimization, we achieved 54% of
SIMD ratio and 56% of efficiency at kernel loop. We
use cell index method to processor parallelization.
It is suitable for the TOFU, torus network of the K
computer. Almost all communication in the cutoff
method is local and well scales on the torus network.
We performed the cutoff method with cutoff length of
28 Å and observed that it completely weak scaled at
6500 atom/node (Table 1). We also measured strong
scaling. The scalability decreased 50% at the case of
number of atom per node was about 100.
We achieved 4.4 PFLOPS calculation performance,
43% of theoretical peak on the simulation of 522
million atoms with 28 Å cutoff using 79,872 nodes.
Table 1 weak scaling: time consumption (ms/step)against number of node
Node

64

512

4096

32768 79872 82944

Time

109.1 110.5 111.2 111.7 112.4 112.1

Part of the results is obtained by the K computer at
the RIKEN Advanced Institute for Computational
Science.
References
1. Akinori Yonezawa, Tadashi Watanabe, Mitsuo
Yokokawa, Mitsuhisa Sato, and Kimihiko Hirao.
Advanced Institute for Computational Science(AICS): Japanese national high-performancecomputing research institute and its 10-petaflops
supercomputer "K". In State of the Practice Reports,
SC '11, pages 13:1-13:8. ACM, 2011.
2. “Optimization of Molecular Dynamics Core
Program on the K computer.” Y. Ohno, R. Yokota,
H. Koyama, G. Morimoto, A. Hasegawa, G.
Masumoto, T. Narumi, and M. Taiji. In Proceedingsof JSST 2012 International Conference ofSimulation Technology, OS8-1, Kobe, Japan, 2012.

60

61

On-site Support : Number of Visitors

4th Biosupercomputing Symposium

Yousuke OHNO
RIKEN

0

5000

10000

15000

20000

25000

Jul-11

Sep-11

Nov-11

Jan-12Mar-12

May-12

Jul-12NEST
ZZ-EFSI
ZZ-HIFU
Platypus
IOSSIM
ParaHaplo
UT-Heart
cppmd

NEST

ZZ-EFSI

ZZ-HIFU

Platypus

IOSSIM

ParaHaplo

UT-Heart

cppmd

Parallelization : Number of Processors• Support other team
• Development core programs/libraries– MD core program• Kernel tuning• Parallel efficiency• Peak performance

Number of Nodes

Optimization of Life-Science applications
on the K computer

Abstract

4th Biosupercomputing Symposium (December 3 - 5, 2012)

62• SIMD, Software-pipelining– Loop structure• Simple loop
• “if” statement → mask operation– Data structure• Simple C style array• STL vector → C style array– SIMD 55 % of all instruction, 99.9% of FP
instruction– Efficiency 56 %, 72 GFLOPS / CPU

Kernel Optimization• Pair-wise Force/Potential– Main component of calculation– O (N2) → O(N) : cutoff method

MD Core Program

32,76879,872

82,944

5.312
0.003
39.0

Performance(PFLOPS)
Efficiency (%)

12.124

Communication
Other

91.622Force & Potential

38.5

0.025

5.313

13.69491.528110.535

38.3

0.201

5.574

14.01191.601111.186

38.1

1.599

5.332

14.01292.329111.672

37.9

3.871

5.332

15.82091.262112.414

38.0

4.031

5.376

15.06791.641112.085

109.058

4,096

418,707 3,349,656 26,797,248 214,377,984 522,546,336 542,644,272

512

Total (ms/step)

64

Number of Atom

Number of Node

Weak Scaling of Cutoff#pragma omp parallel forfor (i=0;i<N_mycell;i++){ // multi threadfor (l=0;l<nl[i];l++){// SIMD, software pipeline
j = list[i][l];dx = x[i]-x[j]; dy = y[i]-y[j]; dz = z[i]-z[j];r2 = dx*dx + dy*dy + dz*dz;
if(r2<cutoff2){mask=1.0}else{mask=0.0} // mask_r = 1.0/sqrt(r2);_r2 = _r * _r; _r6 = ...
e += (LJP(_r12,_r6,) + c[i]*c[j]*CLP(_r,))*mask;
dp = (LJF(_r14,_r8,) + c[i]*c[j]*CLF(_r3,))*mask;
fx[i] += dx*dp; fy[i] += dy*dp; fz[i] += dz*dp;
fx[j] -= dx*dp; fy[j] -= dy*dp; fz[j] -= dz*dp;}}

Kernel Loop

4th Biosupercomputing Symposium (December 3 - 5, 2012)

63

4.11113

Other

3.57428

2.792047.60168

2.12656

Communication

30.350729.170129.4784

4.10638

8.68757

4.120148.36386

42.8347

Force & Potential

41.9641

41.1913

4,096

418,707 33,49,656 26,797,248

512

Total (ms/step)

64

Number of Atom

Number of Node4.213417.401023.9061531.479343.0937

60,293,808

9,216

32,76873,7284.12298

8.21423

4.38367

31.712544.0497

4.15676

8.20732

4.24487

31.9837

44.3478

214,377,984 482,350,464

Weak Scaling of FMM

Strong Scaling of Cutoff

Part of the results is obtained by the K computer at the RIKEN Advanced
Institute for Computational Science.

• Our team had supported optimization of applicationon the K computer.
• We developed MD core program for massive parallel
computer.
• Our MD core program has good scalability.
• Our MD core program had performed 4.39 PFLOPS,
42.% efficiency on the K computer.

Conclusion• 79,872 node
• 522,546,336 atom
• Coulomb and LJ with shift function• 28 Å cutoff• Calculation of potential : every time step
• Time consumption : 116.4 ms / step
• Performance : 4.39 PFLOPLS• Efficiency : 42.9 %Best Performance

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Research activity of Cell Scale Simulation
Team
Hideo Yokota
Team Leader, Cell-scale Research and Development Team, RIKENTeam Head, Bio-research Infrastructure Construction Team, RIKENProfile:2012-present:Deputy Team Leader, Measurement information Laboratory, RCIIC, RIKEN
2011-present:Team Head, Bio-research Infrastructure Construction Team, ASI, RIKEN
2007-present: Team Leader, Cell-scale Research and Development Team, ISLiM RIKEN
2006:Team leader: Bio-research Infrastructure Construction Team, VCAD, RIKEN
2003:Team leader: VCAT development team, RIKEN
1999:Contract researcher: Computational biomechanics unit, RIKEN
1999: Doctor of Engineering degree from the University of Tokyo
1993:Researcher, Kanagawa Academy of Science and Technology
1993: Graduated from Nihon University of Agriculture and veterinary medicine(M.S.(agri.))65

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Cells are the smallest units of life. The ultimate
goal of our cell-scale simulation team is to recreate
intracellular and intercellular phenomena in computer
simulations. However, most intracellular phenomena
have yet to be elucidated, and remain a major areaof research in cell biology. It is therefore currently
impossible to develop a simulator that can recreate
all intracellular phenomena. Hence, our cell-scale
simulation team’s objective is to simulate intracellularand extracellular phenomena that are reasonably wellunderstood and for which mathematical models havebeen established. In addition, using the computational
power of a next-generation 10-petaflop “Kei”
supercomputer, our team aims to simulate intracellularenvironments that heretofore could not be calculated
due to high computational costs. Future advances incell biology, mathematical model construction, andsimulation technologies should markedly expand thescope of cell simulation. In the hope of contributing to
the development of cell simulation, a cell simulation
platform that can analyze coupled phenomena bylinking simulators for multiple phenomena (RICS)
was developed. The cell phenomenon simulation
developed by our team aims to recreate actual
intracellular phenomena, and not merely run computer
calculations. Experiments have been performed
concurrently to gather the necessary parameters forthe simulation and to verify the simulation. This
multi-pronged approach will allow the simulation
to discover unknown phenomena, rather than justrecreate phenomena.

medium inside each voxel. The medium represents
organelles, and the system was designed to deal
with organelle movements in the future. The system
was also designed to withstand large-scale parallel
processing by appropriately compartmentalizing
calculation spaces. In addition to designing the cell
simulation platform, metabolism simulation was
also carried out. In metabolism simulation, “E-cell3”
was used because of its track record in intracellular
metabolism simulation. Furthermore, since membrane
functions are important cellular functions, thefollowing cell membrane functions were added:"substance penetration", "penetration on/off switch(ion-channel)", "active transport (pump)", and "cellmembrane receptor function (internal enzymes
activate when substances bind outside)". Moreover,a function to analyze the membrane potentials ofneurons and other cells is being developed. These
multiple phenomena were examined by weak coupledanalyses based on sequential computations per timestep.
3. Conclusions
We have developed a system to enable spatiotemporal
simulation of the cell. The resulting RICS system
was designed not to only to simulate specific cell
phenomena, but to be a universal spatiotemporal
simulator of the cell. In addition to enabling
simulation of biochemical reactions, diffusion,transport, membrane function, membrane potential,
and advection within the cell, the RICS platform
can couple simulations of multiple phenomena.
Furthermore, with the ability to create geometricmodels that replicate the shapes of actual cells, it
becomes possible to create simulations that accountfor the physical geometry of individual cells and toverify simulation results against actual observations.
Because this approach is computationally expensive,
it would have been extremely difficult to do using
conventional computers. However, spatiotemporal
simulations of cells are now possible by utilizing thelarge-scale parallel processing power of the nextgeneration “Kei” supercomputer.

2. Integrated cell simulation platform
The basic design and methodologies of a common
cell simulation base that can couple intracellular
fields and multiple simulations were investigated.
Basic fields inside the solver were 3-dimensional
scalar fields consisting of the number of various
molecules and the volume ratio of the medium.
As a general rule, different simulators working
on the common base calculate temporal changes
in the number of molecules and the volume of

66

67

• Simulation of cells cluster• Integrated simulation above phenomenon.

ｍｍ• Blood flow (Computational Fluid dinamics)• Membrane potential（Hodgkin-Huxley）– ion channels, transporters, receptor, et al.

• Membrane Functions（Selective permeability）– Membrane transport• Membrane transport（Active transport）– concentration gradient• Diffusion（Passive transport）– E-cell• Biochemical reaction（Metabolism）

ｎｍ

1

RIKEN
Integrated
Cell Simulation Platform
Novel 3D
cell simulator

RICS
Targets of cellfunctions

4th Biosupercomputing Symposium

Hideo YOKOTA
RIKEN

Research activity of Cell Scale
Simulation Team

Olganera3

E-cellOlganera2

E-cellBiochemical reaction

E-cell olgnera1

Portal vein

Transport

t=t+ΔtMembrane

Product1

DifusionTranspoter

Product3

t=t+ΔtProduct2

Coupled simulaton in RICS

Red blood cell: homogeneous

Metabolic path way

Almost cell: Un-homogeneous

Microscopic ImageRIKEN BSI MiyawakiLab.

Liver

Hepatic lobules

Need :spatiotemporal simulation of the cell

2

Central vein

4th Biosupercomputing Symposium (December 3 - 5, 2012)

A

H+,Na+
Ca2+,Cl-

4096 voxels

Ions

Large
molecules

channelTranspoter

membrane

Integrated enzymes

Receptors

Region of biochemical reaction

68

Vmax = 1.0
Km = 0.5

Ez1

B
Vmax = 2.0
Km = 1.2

Ez2
C

extracellular area
Substance A：100 mM

Cell membrane
Enzyme E2：1 mM

CenterEnzyme E1：1 mM

Cytosol
Substance A: 0 mM
Substance B: 0 mM
Substance C: 0 mM

Simulation of membrane transportwith bioreaction

Pump

Amino acids
Glucose

H2O
EtOH

O2
CO2

Selective permeability

Functions of membrane

RICS
Preprocessor

Sectional
images©RIKEN：Live Cell Modeling Project

イメージングデータ

3D rebuilding Cell data frommicroscopic images

Virtual cell shape (CAD,CG)

Section of a cell in voxels

RICS cell model

RICS cell modelTransport of the substance through the membrane is the sum of the flow velocity foreach transport.

SW = regulatory function of each transporter

Func = {penetrate, channels, pumps, transporters, gap junctions}f m = ∑ Func SW Func f mFunc∂φ= ∇ • (C∇φ ) − ∇ • f m + [ BioChemical Re action]∂tFlux of membrane transporters

4th Biosupercomputing Symposium (December 3 - 5, 2012)

69

Summary

diffusion＋reaction（apparent diffusion rate 28μm/s）
Only diffusion（diffusion rate 32μm/s）

11Open source release : December 2012

: large-scale parallel
more than 24,000nodes

・ 京(kei) available: verify simulation results against actual observations

・Real cell shape simulation:biochemical reactions, diffusion, transport,
membrane function, membrane potential, and advection

・Couple simulations of multiple phenomena:Spatiotemporal simulation of the cell

・RICS (Riken Integrated Cell Simulation Platform)

reaction：24
material：10

Ca2+reaction modelＣa+ kinetics simulation （ HepG2）

9Perisinusoid

C

Tadao Shibasaki
Susumu Seino

Kobe Univ.

Shingo Murakami
Hiroshi Hibino
Yoshihisa Kurachi

Taiji Adachi

Koichi Takahashi

Yuki Tsujimura
Shinya Goto

Osaka Univ.

Takashi Koyama

Shigeho Noda

Shin Yoshizawa

Yasuhiro Sunaga

RIKEN

Yohei Nagasawa

Tokai Univ.

Makoto Suematsu

Ayako Yachie

Keio Univ.

Acknowledgement

Hepatic cords
VCAT

Hepatic lobules

Real shape

Central vein

Microscopic image

Hepatic arteryPortal vein

Bile duct

Capillary

Liver cell simulation

Liver

10

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Dissection of regulatory mechanisms formetabolic systems by quantitative imaging
mass spectrometry
Makoto Suematsu
Professor, Department of Biochemistry, School of Medicine,
Keio University, JST ERATO Suematsu Gas Biology Project

Profile:2009-present: Leader, Japan Science and Technology Agency, ERATO Suematsu Gas Biology Project
2007-present: Dean, School of Medicine, Keio University
2007: Leader, Global Center of Exellence for Life Sciences, Human Metabolomic Systems Biology from MEXT
2003: Leader, National Leading Project for Biosimulation by Ministry of Edication, Sciences and Technology
2001: Professor and Chair, Department of Biochemistry and Integrative Medical Biology, Keio University Schoolof Medicine
1991: Bioengineer Step IV, Institute for Biomedical Engineering, University of California San Diego (Supervised
by Professor Benjamin W Zweifach and Professor Geert W Schmid-Schoenbein)
1983: Graduated from Keio University School of Medicine (MD)71

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
Collection of the data in reality is absolutely
necessary to empower the capability of large-scalebiosimulation of metabolic systems, while it has
been difficult to collect quantitative data of manyof metabolites in tissues in vivo. We have appliedadvanced metabolome technology including imaging
mass spectrometry to this effect. Application of
metabolome analyses based on CE-MS allowed us topredict and demonstrate roles of hemoglobin allosteryin O2-sensing mechanism in human erythrocytes, andnow extending the technology to analyses of cancer
metabolism in vivo. Metastatic progression of cancer
does not only upregulate their glucose metabolism but
might utilize metabolic properties of host organs that
benefit cancer metabolism, although such a hypothesis
remains elusive. Newly developed microscopic
imaging mass spectrometry combined with MS2
analyses allowed us to collect micrographs of manydifferent metabolites in a single frozen section, andcombination with CE-MS data collected from theserial section provided semiquantitative informationof individual signals. The current method revealed
that human colon cancer xenografts metastasized
in livers of super-immunodeficient NOD/scid/γnull(NOG) mice deprives L-alanine to support their
metabolic demands for synthesizing glutathione
and nucleotides. In this model, hepatic metastasis
triggered regenerative responses of the host liverconcurrently with hypoglycemia and accumulationof glutathione and nucleotides in the tumor-bearing
liver. MS2 analyses under loading 13C3-L-alanine
provided evidence for earlier filling of glutathione
with 13C2-γ-glutamylcysteine structure in metastases
than surrounding liver parenchyma. The 13C3-Lalanine loading also caused an increase in 13C2-UDP
in metastatic foci and in the host liver. MS2 analyses
to assess the pathways for 13C incorporation revealed
that L-alanine not only undergoes gluconeogenesis in
the host to synthesize ribose but serves as a substrate
to supply glutamate and pyrimidine carbons fornucleotide synthesis occurring in the metastases. Our
results suggest that human colon cancer metastases
utilize gluconeogenic substrates of the host not only

through pentose phosphate pathway but throughglutaminolysis, supporting their metabolic demand of
glutathione and nucleotides to cause hypoglycemia.
In order to explore novel molecular targets that play a
crucial role for survival of cancer cells, combinationof human cancer xenograft models with large-scalecomputer simulation of metabolic systems mightserve as a powerful stratagem in future.
References
1. Soga, T., et al. Differential metabolomics reveals
ophthalmic acid as an oxidative stress biomarker
indicating hepatic glutathione consumption. J.
Biol. Chem. 281(24) 16768-16776, 2006.
2. Kinoshita, A., et al. Roles of hemoglobin allosteryin hypoxia-induced metabolic alterations inerythrocytes: simulation and its verification bymetabolome analysis J. Biol. Chem. 282(14),10731-10741, 2007.
3. Hattori K, et al. Paradoxical ATP elevation inischemic penumbra revealed by quantitative
imaging mass spectrometry. Antioxid Redox
Signal (News & Views), 13(8), 1157-1167, 2010.
4. Ishimoto T, et al. CD44 variant regulates redox
status in cancer cells by stabilizing the xCT subunitof system xc- and thereby promotes tumor growth.
Cancer Cell 2011, 19(3), 387-400.
5. Morikawa T, et al. Hypoxic regulation of thecerebral microcirculation is mediated by a carbon
monoxide-sensitive hydrogen sulfide pathway. Proc
Natl Acad Sci USA 2012, 109(4), 1293-1298.

72

Copyright © 2006 Keio University

ISLiM Symposium 2012

Makoto Suematsu, MD, PhD
Professor & Chair
Department of Biochemistry, School of Medicine, Keio University
Japan Science Technology Agency, ERATO Suematsu Gas Biology
gasbiology@z6.keio.jp~ Application to ischemic diseases and cancer models ~Dissection of regulatory mechanisms for metabolic
systems by quantitative imaging mass spectrometry

0

50

100

150

0

100

200

300

I = median

Ii

IMS

caudalPeak height of ATP mass spectrum

ATP

C’

CE-MS

bregmarostralIi

53 columns x 35 rows

from all the spots in the contralateral hemisphere

I: the median of maximum intensities of a metabolite

specified range at the ith spot

Ii: the maximum intensity among mass spectra in acontralateral hemisphere determined by the CE/ESI/MS

Peak height of AMP mass spectrum

AMP

I

C’

C’: the metabolite content of tissue from a

Ci =

Ii

+

+
+
+

+
+

+

+

+

+

+

+

+

+

+

+

MS

Agilent Technologies

Ionizedmetabolites

+

Laser

53 columns x 35 rows

N

N

Fe
N

N

HOOC

CO

HOOC

Heme oxygenase

HO-2

HO-1

0

2

4

6

8

10

0

O
N HO-1ON

Km =15 µM

2

50
100 150 200 250
Biliverdin
IX α
O (µM)Km 35 µM =20mmHg

N

N

HO-2

Fe (II)

Suematsu M, et al. BBRC 1994, J Clin Invest 1995
Goda N, et al. J Clin Invest 1998
Kyokane, et al. Gastroenterology 2001

Vasodilation↑ cGMP

Ishikawa M, et al. Circ Res 2005
Shintani T, et al. Hepatology 2009
Kajimura M, et al. ARS 2010
Morikawa T, et al. PNAS 2012

↓ H 2S

Vasoconstriction

Liver: soluble guanylate cyclase (sGC)
Brain: cystathionine β-synthase (CBS)Microvascular pericytes in sinusoids (Ito cells) Astrocytes in neurovascular units

Protoheme IX

HOOC

HOOC

+

Cytokines, ROS
cAMP, heavy metals, etc

Liver and Brain: Gigantic resources of CO with distinct gas-responsive receptors

O2

Matrix
http://biopedia.org/images/1/1e/Maldi-TOF-MS.jpg

matrix

Imaging data

Ionizedmetabolites

Laser

MS

Soft Ionization by MALDI Intensity map

MALDI-IMS

Content map

Hattori et al., Antioxid Redox Signal, 2010+

Quantitative data

Electrospray

Capillary
Electrophoresis

Complete Ionization by ESI

CE/ESI/MS

bregmaInitial velocity

Semi-quantification: apparent contents of a specific metabolite at theith spot of tissue (Ci) was estimated:

Frequency(x 10-3 min-1)

73

Frequency

Semiquantitative imaging mass spectrometry: IMS Combination with CE/ESI/MS

4th Biosupercomputing Symposium (December 3 - 5, 2012)74

cerebellum

HO2/GFAP

Mayumi Kajimura

Takayuki Morikawa(D1) Katsuji Hattori(D3)GFAP
CBS/GFAP

ml

gl

10 µm

0

2

4

6

6

8

8

10 0

10 0
6

Time (min)4** * * **** ** * ** *

2

8

10

0
4

*

6

*

*

8

10

Morikawa T, et al. PNAS 2012 Jan 9

2*

CBS null (n=3)

4** * * * **** *** *** **HO-2 null (n=10)

CSE null (n=6)

2*

WT (n=12)

10 % O2 inhalation

Pre-capillary arterioles

20

0

10 % O2 inhalation

Penetrating arterioles

40

60

80

0

20

40

60

80

(c)

Hypoxia-induced microvascular dilatation in brain depends on CO-CBS/H2S system

Hypoxia (g)HO-2 nullTime (min)xMovie: 1-min
stepwisez
21% O2
10% O2(e) Baseline

(d)

a
b
c
d

(b)

WT

CBS

gl

H2S-producing site

Hypoxia-induced vasodilatation seen through gas-tight thinned-skull preparation

10 µm

(a)

HO2

gl

Pl

P
l

* *

ml

ml

CO-producing site

% Dilation

*

*

*

*

*

CBS

Astrocyte

H2S

Endfoot

HO2

1

10

1

6

21

O2 fraction (%)

21

*

6

Morikawa T, et al. PNAS 2012 Jan 9

0

1

2

3

CBS

HO2

H2S

CO

Neuron↓CO from HO2

EC

Endfoot

Smallarteriole～ATP

Respiratory lock↑ ATP

unlocked

Neural PO2(sensing by HO-2 )

20 µm

CBS

CBS

H2S

CBS

~CO from HO2~Neural PO2

HO-2-null mice increase basal ATP under normoxia and are intolerant against hypoxia

EC
20 µm
vasodilation

CBS

10

*

O2 fraction (%)*

WTHO-2 nullHS- measurements in ISF-brain by LC-MS2
with bimane-assisted delivertization

0

5

10

CO measurements in brain slice

NeuronCO

Neural PO2(sensing by HO-2 ?)

**

**

Small
arteriole*** *Time (min)

**

**

*

WT (n=14)
CBS-KO (n=10)
HO2-KO (n=8)

Morikawa T, et al. PNAS 2012 Jan 9

Arteriolar diameter (% change)

Either HO-2 null or CBS null significantly diminishes hypoxia-induced microvascular dilation in brain slices

CO(pmol/mg brain)
SDB(pmol/mg brain)

HYPOTHESIS: Heme oxygenase（HO) senses a drop of O2 to inhibit COthat might facilitate CBS-derived H2S generation to dilate arterioles
Neuron produces CO
Astrocyte produces H2S from CBS

4th Biosupercomputing Symposium (December 3 - 5, 2012)

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

MS

Agilent Technologies

Ionizedmetabolites

MALDI-IMS

B

Citrate (µmol/g tissue)

0

0.4

0

0.4

0.8 0.8TT

i.sp 2wks

Pc
Pc

*

0
0

6.0
6

12.0
12

18.0
18

C
Ctrl

CitrateTTi.sp 2wks

Pc
Pc

*†

0

0.4

0.8

0

0.4

0.8

0

0

4.0
4

8.0
8

12.0
12

Glu

Laser

Soft Ionization by MALDI

TT*†

i.sp 2wks

Pc
Pc

*

0

6.0

12.0

18.0

0

6.0

12.0

18.0

0

0

4.0
4

8.0
8

12.0
12

GSH

C
Ctrl

TT

†

0

4.0

8.0

12.0

0

4.0

8.0

12.0

i.sp 2wks

Pc
Pc

*

matrixMatrix

Imaging data

Ionizedmetabolites

Laser

MS

0

0

1.0
1

2.0
2

33.0

CtrlC

GSSG

TT

*

i.sp 2wks

Pc
Pc

*

0

4.0

8.0

12.0

0

4.0

8.0

12.0

10.0
10

0

0

5.0
5

C
Ctrl

GSO3H

TT

*†

i.sp 2wks

Pc
Pc

http://biopedia.org/images/1/1e/Maldi-TOF-MS.jpg

Content map

Ctrl
C

Intensity map

Kajiwara Ogawa

GSH is highly accumulated in metastatic cancer fociMechanisms for detoxifying anti-cancer reagents ?

C
Ctrl

*

i.sp 2wks-F

A Control-F

Glutamate (µmol/g tissue)

Takeshita Aoki

Shimadzu: β-2 IMS: 5-micron resolution
Fine spraying for matrix coating

CE-TOFMS

Hattori K, et al., ARS 2010, Kubo A, et al. ABC 2011++

Quantitative data

Electrospray

Capillary
Electrophoresis

Complete Ionization by ESI

GSH (µmol/g tissue)←GSH 306.08←UMP 323.03←AMP 346.05←G-SO3H 354.07

350↓UDP 403.00

400
m/z

←ADP 426.02

450

x10

←UTP 482.96↓ATP 505.99←taurocholate 514.28

500←GSH 306.08

300↑UMP 323.03↓AMP 346.05

350↑G-SO3H 354.07↓UDP 403.00

400
m/z

Tumor-bearing liver (isp. 2w)
x10

450

Intensity A.U.

13C -γ-Glu-Cys
2

13C -Alanine (from host liver)
32-oxoGlu (TCA cycle)

Ishimoto et al. Cancer Cell 2011
Tamada et al, Cancer Res 2012
Yae et al. Nature Commun 201215N-cystine to 15N-cysteine in CD44(+) cancer

Pyruvate

glutamate

CD44, a cancer stem cell marker, is necessary to express xCT on the cancer cell membranethat facilitates GSH synthesis through extracellular resources

H2S

500

↓ATP 505.99

GSH and high-energy nucleotides:
major metabolites in tumor-bearing liver in vivo

300

Control NOG liver

Xenograft transplantation of human-derived cancer cells to super-immunodeficient NOG miceallows us to reconstitute solid tumor metastases in vivo

←ADP 426.02

Nakanishi Harada

GSSG (µmol/g tissue)↓UTP 482.96

CE-TOFMS for quantifying signals

GSH/GSSG

75
GSH+2GSSG(µmol/g tissue)←taurocholate 514.28

Semiquantitative imaging mass spectrometry: IMS Combination with CE/ESI/MS

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Closing the loop: simulation of the wholesensory-motor neural network in action
Kenji DoyaProfessor, Okinawa Institute of Science and Technology

Profile:KENJI DOYA took BS in 1984, MS in 1986, and Ph.D. in 1991 at U. Tokyo. He became a research associate at
U. Tokyo in 1986, U. C. San Diego in 1991, and Salk Institute in 1993. He joined ATR in 1994 and became thehead of Computational Neurobiology Department, ATR Computational Neuroscience Laboratories in 2003. In
2004, he was appointed as the principal investigator of Neural Computation Unit, Okinawa Institute of Scienceand Technology (OIST) and started Okinawa Computational Neuroscience Course (OCNC) as the chief organizer.
As OIST re-established itself as a graduate university in 2011, he became a professor and the vice provost for
research. He serves as the co-editor in chief of Neural Networks from 2008. He is interested in understanding thefunctions of basal ganglia and neuromodulators based on the theory of reinforcement learning. Contact: doya@oist.jp, 1919-1 Tancha, Onna, Okinawa 904-0495, Japan.

77

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Realistic simulation of the brain is a grand challenge
toward understanding how the brain works and also
in establishing how to simulate complex dynamicalsystems. Detailed computer simulations of the brain
have so far been limited to the cellular or local
circuit levels, such as the Blue Brain Project andCortical Columns in Silico. The aim of our projectis to reconstruct the whole brain network including
sensory inputs and motor outputs and to test how thebrain functions under dynamic interactions with theenvironments.
The use of super computers has mostly been limited
to batch operations in which huge jobs are put into aqueue and the results are given after indefinite hours
or days waiting. An important rationale for a locallyconcentrated super computer like K, as opposed toworld-wide distributed computing grids, is its use
for time-critical missions. Simulation of the brain ofa humanoid robot is a scientifically important andcomputationally demanding challenge. Despite big
progresses in movement control and appearances,current humanoid robots are far from human-likebecause of their limited behavioral responses, inother words, because of their poor brains. Humanlike behaviors require huge computations, fromsensory processing and motor control to language
and emotion. And natural human interaction requires
timely responses. Real-time simulation of the brain ofa humanoid robot is an ultimate test of what we know
(and do not know) about the brain function and also
an ultimate test in high-throughput computing.
As the first step in implementing a humanoid brain,
we focus on the oculomotor control loop, which isone of the best-studied model system in neuroscience.
2. Oculomotor Network Model
The mammalian oculomotor system has a hierarchicalorganization with the retina – superior colliculus –brain stem – eye muscles as the most basic pathway,
over which the cerebellum, the basal ganglia, and the
cerebral cortex impose more sophisticated control.
Here we focus on the basic pathway with the emphasison the neural circuit of the intermediate layer of thesuperior colliculus, which is regarded as the center forthe control of saccadic eye-movements.
The superior colliculus model [1] consists of thevisual input layer, burst neuron layer (output), buildup
neuron layer, and deep inhibitory layer. The neurons
were implemented by integrate-and-fire models andthe excitatory and inhibitory connections were set up
according to the anatomical literature and to reproduce
the activity features of different types of neurons.
The model was implemented using a spiking neuralnetwork simulation tool NEST [2]. The model of
one side of the superior colliculus consisted of about
100,000 neurons and the computing time was aboutfour times the simulated biological time on 800 nodesof the Intel cluster at RICC. The result replicated the
known velocity and duration features of saccadic eye
movements, known as ‘main sequence.’ The model
suggested the role of spreading activities of buildup
neurons in deciding when to stop the eyes [1].
3. Robot Eye Movement Experiment
We connected the left and right copies of the superior
colliculus model and the brain stem network modelwith a multi-simulation coordinator MUSIC [3]. The
integrated model ran at RICC was connected throughthe internet to the eye control system of a humanoid
robot CB-i at ATR [4]. The CB-i could interact witha human partner by tracking by its eyes a marker
moved by the partner. Although the K computer was
not available for such an interactive use, our result
demonstrates the feasibility and usefulness of highthroughput computing for real-time simulation of thebrain for live interaction with humans.
References
1. Moren J, Shibata T, Doya K (2010). Toward aspiking-neuron model of the oculomotor system.
International Conference on Simulation of Adaptive
Behavior, 2010, Paris. (Springer Lecture Notes in
Computer Science 6226, pp. 104-113.)
2. http://www.nest-initiative.org
3. http://software.incf.org/software/music.
4. http://www.cns.atr.jp/bri/en/robot/cb-i/

78

79

4th Biosupercomputing Symposium• On- line simulation• Off- line simulation• Superior colliculus model• Retina model

Outline

4th Biosupercomputing Symposium

Kenji Doya
Okinawa Institute of Science and Technology

Closing the Loop:
Simulation of the Whole SensoryMotor Neural Network in Action

3

4th Biosupercomputing Symposium• by Usui Group at RIKEN BSI– eye optics– photo receptors– interneurons– ganglion cells

Retina Model

4th Biosupercomputing Symposium

42• Connect a supercomputer with a robot– live human interaction• Simulate the wholesensorimotor loop
– oculomotor network

Goal of the Project

4th Biosupercomputing Symposium (December 3 - 5, 2012)

80• Activity

4th Biosupercomputing Symposium– near target: higher rate
– far target: longer duration– constant total spikes

Physiology

4th Biosupercomputing Symposium• The center of saccadic eye movement– superficial layer : sensory input processing– intermediate layer: motor pattern generation•retinotopic coordinate•spreading activation

Superior Colliculus (SC)75• Peak velocity

4th Biosupercomputing Symposium• Duration

4th Biosupercomputing Symposium• Visual input– Superficial SC– Quasi Visual layer• Burst neurons: output– inhibitory inputs•substantia nigra•reticular formation•SC deep layer• Buildup neurons– spreading activity

Anatomy8

6

4th Biosupercomputing Symposium (December 3 - 5, 2012)

81

9

th

• We built the first spiking neuron model ofthe intermediate SC based on experimental
data.
• The model reproduced:– burst neuron spike profiles– buildup neuron spreading activity– saccade main sequence• The model predicts:– spreading activity to stop the eye accurately• Real- time simulation of a robot’s brain– enables interactive
testing of the model
4 Biosupercomputing Symposium
11

Conclusion

4th Biosupercomputing Symposium• MUSIC– multi simulation
coordinator• Real- time simulation– at RICC• Connection over internet– TCP socket• Robot eyes
– CBi at ATR

Closing the Loop
• Outside view

4th Biosupercomputing Symposium• Robot’s eye view

In Action10

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Brain-scale neuronal network
simulations on K
Markus Diesmann
Professor, Director of Institute of Neuroscience and Medicine (INM-6)
Computational and Systems Neuroscience, Juelich Research Centre (Germany)

Profile:Markus Diesmann received his Diploma and PhD in Physics from the University of Bochum (Germany) in 1994
and 2002. The PhD work was conducted at the Weizmann Institute of Science, Rehovot (Israel) and continued
at the University of Freiburg (Germany). In 1999 he joined the Max-Planck-Institute for Dynamics and SelfOrganization, Goettingen (Germany) for a staff position. From 2004 to 2006 Markus Diesmann was appointed
Junior professor for Computational Neurophysics at the University of Freiburg. In September 2006 he moved
to the RIKEN Brain Science Institute, Wako (Japan) for the position of a Unit and later Team Leader. Since
March 2011 Markus Diesmann is Professor for Computational Neuroscience at the Medical Faculty of RWTH
Aachen University (Germany) and serves as the director of the Institute of Neuroscience and Medicine (INM6) Computational and Systems Neuroscience, Juelich Research Centre (Germany). His interests include thecorrelation structure of cortical networks and large-scale simulations.
www.csn.fz-juelich.dewww.nest-initiative.org

834th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
The human brain comprises about 1011 neurons thatare sparsely and specifically connected by about
10,000 outgoing synapses each, mediating electricalpulses (spikes). In computational neuroscience, thebottom-up approach starts from a mathematicaldescription of the single neurons and their interactions
in order to investigate the emergent network dynamics[1]. The NEST simulator [2] is tailored to thisresolution. Individual neurons with a single or few
compartments are represented as small systems ofdifferential equations, which interact by δ-impulses toform coupled networks of natural size and complexity.
For the neuron models frequently used in the field,
the time evolution of the (essentially linear [3])dynamical equations can often be integrated exactly,
requiring only a few floating point operations perneuron to implement precise update schemes [4]. The
mammalian cortex, which carries out higher brain
functions, is of particular interest. A network of about
105 neurons, called the local cortical microcircuit,is required to represent the layer- and cell-type
specific connectivity with all local synapses and canbe interpreted as an elementary unit of the cortex.
Fundamental dynamic properties, like the layerdependent firing rate are explained at this scale [5].
The top-down approach starts from an abstract
description of a particular brain function andinvestigates how this function is implemented at
the level of neurons and synapses (e.g. [6]). The
functional circuits of the cortex typically involve
several brain areas [7]. The well-studied visualcortex of primates comprises on the order of 100
million neurons, hierarchically organized in areasthat can exceed 10 million neurons and 100 billion
synapses. Previously, simulations were constrained
to models of the local cortical microcircuit or nonspiking macroscopic models. With the technologicaladvances of recent years, the representation of the fullmulti-scale connectivity of the brain has come withinreach. The memory demands of such simulations [8]are only met by distributed simulation software andsupercomputers, such as the K computer in Kobe.

2. NEST on peta-scale supercomputers
NEST is subject to incremental and iterative
development driven by the needs of the neuroscience
community. These efforts are coordinated by the
NEST Initiative (www.nest-initative.org), a non-profit
organization. The storage of synaptic connectionsof brain-scale networks dominates the memory
consumption of spiking network simulations, because
synapses outnumber neurons by a factor of 104. Only aparallel computer offers the required memory. In NEST
the synapses are stored on the machine where the target
cell is located. This scheme keeps the communicationload between processors low, as only the identifiers ofneurons that emitted a spike need to be communicatedto the other machines, employing collective MPI
communication.
The kick-off meeting of the Brain and Neural Systems
Team (BNT) took place on November 10th, 2008, at
the RIKEN Marunouchi office in Tokyo. First trials of
the simulation kernel developed for small clusters (2ndgeneration kernel, 2g) executed on the BG/L at RIKEN
Wako scaled to 1024 processors for a network of 105
neurons. At the 3rd BNT meeting October 6th, 2009,
we presented scaling data for 106 neurons and up to
32,768 core of the JUGENE BG/P computer at Juelich.
These data exposed memory consumption as thelimiting factor. To investigate the issue systematically,
we developed a mathematical model which describes
the memory usage per process as a function of the total
number of processors and network size [8]. Guided by
the model, we identified the data structures dominating
the memory usage of NEST. On supercomputers, the104 or more compute nodes outnumber the synapses
per neuron, such that on average each node stores
only one synaptic target per source neuron. Moreover,
only a small fraction of the total number of neurons
is local to a given node. By February 2010 we had
defined a roadmap for the K computer where along
3 milestones network size increases from 106 (1) to108 (2) and above (3). On September 28th, 2010, just
2 days before the 5th BNT meeting, the simulation
language interpreter compiled on a prototype machine.
On November 12th 2010 the complete NEST compiled

84network sizefor the first time as reported at the 6th meeting on
April 18th of the following year. On May 6th, 2011,
NEST compiled on the K test system and finally on
September 1st on K as reported at the 7th BNT meetingon the 29th of the same month. In the 3g kernel,
we employ data structures which account for thesparseness of synapses and enable the efficient storageof information about non-local neurons [8]. With these
improvements to the fundamental data structures thesecond milestone came into sight at the 8th meeting on
March 13th, 2012, and was reached in May [9] utilizing
just above 12,288 compute nodes which is less than
14 per cent of K. The simulation of realistic neuronal
networks requires the representation of different formsof synaptic plasticity [10,11]. The outgoing synapsesof a given neuron are stored in a set of homogeneous
data containers. Distinguishing different synapse types
therefore requires an intermediate vector-like data
structure. However, this distinction is inefficient forcomputers with 104 nodes, where typically for a source
neuron only a single synapse is stored. In the second
step of our redesign we therefore built an adaptive
framework which chooses the optimal container froma predefined set depending on number and type of thesynapses. The set of possible containers is created byrecursive template meta-programming [12]. In the third
step, we improved the memory layout of the synapse
objects. When a connection is established, it is checked
whether the involved neurons and synapse types match.
The improved handshake mechanism does not require avirtual function table pointer in the synapse. In addition
the alignment of synaptic parameters was improved.
The last two steps of the redesign are combined in
the 4g kernel presented at the 9th BNT meeting on
September 25th. Fig. 1 shows the network size whichjust fits on a given number of compute nodes and theruntime for the corresponding network simulation.
The 4g kernel enables simulations of 109 neurons on
K. The increase of runtime with increasing number ofnodes is due to the collective communication scheme.
Good scaling of network setup and simulation phase
are achieved with a hybrid code which combines finegrained parallelism employing OpenMP threads anddistributed parallelism based on MPI. The 3g and

109

106

108

105

107

104

106

103

64

6

24

runtime (s)4th Biosupercomputing Symposium (December 3 - 5, 2012)4
6
38 55316
6
number of compute nodes
25

10

96

40

Fig 1. Maximal network size (triangles) andsimulation time (dots) as a function of the number of
employed nodes for the 4g kernel.
4g kernels do not compromise on generality of NEST;
functionality and user interface remained the same [13,14].
Acknowledgements
Partly supported by early access to K at RIKEN AICS,
Next-Generation Supercomputer Project of MEXT,
Helmholtz Association: HASB and portfolio theme
SMHB, JARA and EU Grant 269921 (BrainScaleS).
References
1 Lansner and Diesmann (2012) Chap. 10 In
Computational Systems Biology, Le Novére ed.,Springer2 Gewaltig and Diesmann (2007) Scholarpedia,
2(4):1430
3 Morrison et al. (2008) Biol. Cybern. 98:459-78
4 Hanuschkin et al. (2010) Front. Neuroinform. 4:113
5 Potjans et al. (2012) Cerebral Cortex, in press
6 Potjans et al. (2011) PLoS Comput Biol 7(5):e1001133
7 Wagatsuma et al. (2011) Front. Comput. Neurosci.
5:31
8 Kunkel et al. (2012) Front. Neuroinform. 5:35
9 Helias et al. (2012) Front. Neuroinform. 6:26
10 Potjans et al. (2010) Front. Comput. Neurosci.
4:141
11 Pfeil et al. (2012) Front. Neuromorph Eng 6:90
12 Abrahams and Gurtovoy (2004) C++ Template
Metaprogramming, Addison-Wesley Longman
13 Eppler et al. (2009) Front. Neuroinform. 2:12
14 Gewaltig et al. (2012) Chap. 18 In Computational
Systems Biology, Le Novére ed., Springer

85

86

excitatory (80%) and inhibitory population (20%)
of neurons in each layercorrect distribution of spike rates
across layers



4th Biosupercomputing Symposium

Potjans TC & Diesmann M (2012) The cell-type specificconnectivity of the local cortical network explains prominent
features of neuronal activity. Cerebral Cortex, in press

asynchronous-irregular spikingof neurons



taking into account layer and neuron-type specificconnectivity is sufficient to reproduce experimentallyobserved:

105 neurons and 109 synapses

Local cortical microcircuit

4th Biosupercomputing Symposium

Brain and Neural Systems Team

Markus DIESMANN
Jülich Research Centre

Brain-scale neuronal network
simulations on K3post-synaptic neuron responds with small
excursion of membrane potential
inhibitory neurons (20%) cause negative
excursionemission rate of 1 to 10 spikes
per second

causing large fluctuations ofmembrane potential

4

building block for further studies: Wagatsuma N, Potjans TC, Diesmann M and Fukai T(2011) Layer-dependent attentional processing by top-down signals in a visual cortical
microcircuit. Front Comput Neurosci 5:31


4th Biosupercomputing Symposium

Potjans TC & Diesmann M (2012) The cell-type specific connectivity of the local cortical
network explains prominent features of neuronal activity. Cerebral Cortex, in press

Cortical flow of activity

2each neuron receives input fromup to 10,000 other neurons

4th Biosupercomputing Symposiumsupra-threshold value causes spiketransmitted to post-synaptic neuron

current injection into pre-synaptic neuron
causes excursions of membrane potentialFundamental interactions

4th Biosupercomputing Symposium (December 3 - 5, 2012)87

teaching in international courses



4th Biosupercomputing Symposium

collaboration of several labs (since 2001)

Collaboration

4th Biosupercomputing Symposium

produce public releases
under GPLsystematically publish new
simulation technology

Major goals:

 connectivity of local microcircuit
At least three levels of organization:  within-area connectivity with space constant long-range connections between areas

Architecture of the human cortex

75

only brain-scale model represents majority of inputs to a neuron
self-sustained activity

brain-scale models required to close functional circuits

6

4th Biosupercomputing Symposium

8

Kunkel S, Potjans TC, Eppler JM, Plesser HE, Morrison A and Diesmann M (2012) Meeting the memory
challenges of brain-scale network simulation. Front Neuroinform 5:35

memory overhead increases with cores
memory not simulation time limits network sizeuse full memory resources: maximum-filling scaling
analysis based on mathematical model of memory consumption:







scale 108 relevant because: larger than largest area (V1)
 enable visual cortex model respecting relative sizes
 larger networks: long delays, sparse macroscopic connectivity

scale-up on K guided by 3 milestones
1. port NEST software to K
(2nd generation simulation kernel)
2. scale of 108 neurons(3g simulation kernel)
3. towards full brain(4g simulation kernel)

Scale-up to networks of 109 neurons

4th Biosupercomputing Symposium

requires work on software technologybrain functions are distributed over many areas

underconstrained model limits predictive power

only 50% of incoming synapses of a neuron come from local microcircuit











Need for brain-scale networks

4th Biosupercomputing Symposium (December 3 - 5, 2012)

88

4th Biosupercomputing Symposium

How to measure scalability

4th Biosupercomputing Symposium data structures account forheterogeneity of synaptic dynamics
 for > 10,000 cores, neurons with few
local targets cause severe overhead novel adaptive data structure copeswith short target lists
 not compromising on generality

4g memory layout

911

Helias M, Kunkel S, Masumoto G, Igarashi J, Eppler JM, Ishii S, Fukai T, Morrison A, Diesmann M(2012) Supercomputers ready for use as discovery machines for neuroscience
Front Neuroinform 6:26

3g memory layout
 accounts for sparseness in neuronal and connection data structures
published as:

Memory layout of 3g and 4g kernel

How to measure scalability

4th Biosupercomputing Symposium

fixed number of neurons per coreideally: constant runtime
 increase due to communication
 again, faster element update leads toworse scaling




10

12

in case of limited memory resources:
 again, reduced memory consumption
improves scaling extreme case: same network can be
simulated faster on fewer cores (dotsindicate equally sized networks)

4th Biosupercomputing Symposium

improved solver
compromising on accuracy of the solution

Weak scaling

Faster integration of neuronal dynamics can be accomplished by:

faster element update leads to earlier
saturation (comm. dominates runtime
already on fewer cores) however, same network can be simulated
faster on fewer cores

fixed size of neuronal network
ideally: simulation time decreases
proportional to number of cores saturation due to communication

Strong scaling

How to measure scalability

4th Biosupercomputing Symposium (December 3 - 5, 2012)

89

4th Biosupercomputing Symposium

concepts for exa-scale machines•

next steps:

make 4g kernel available to partners

short run times enable use of K as discovery machine•no compromise on generalityprepare manuscript on 4g kernel

109 barrier in sight with 4g kernel (milestone 3)

•

108 milestone achievable with 3g kernel (milestone 2)

Summary

4th Biosupercomputing Symposium

neuroscientist interested in maximum use of resources
hold memory close to maximum available on node (16 GB)
neurons per core drop
simulation time increases due to increased communication











Maximum filling scaling of 3g on K

15

13



Abigail Morrison

Marc-Oliver Gewaltig

Hans Ekkehard Plesser

Jochen M Eppler

Susanne Kunkel

Moritz Helias

Yousuke Ohno

Gen Masumoto

Makoto Taiji

Jun Igarashi

Tomoki Fukai

Shin Ishii

4th Biosupercomputing Symposium

Acknowledgements

4th Biosupercomputing Symposium

wall clock time 1 biological second of the largest 4g simulation on K: 1520s = 25min

K largest number of MPI jobs tested 24,576 = 196,608 cores

Jugene largest number of MPI jobs tested 65,536 = 262,144 cores

end of dashed lines indicate machine size

Performance of 3g and 4g kernel
on Jugene and K

16

14

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Molecular machines and nuclear processes
studied by coarse-grained molecularsimulations
Shoji Takada
Associate Professor, Department of Biophysics,
Graduate School of Science, Kyoto University

Profile:1990
Graduate School of Science, Kyoto University (Master degree)
1991-1998 Institute of Molecular Science, (Technical officer)
1998-2007 Faculty of Science, Kobe University (Lecturer, Assoc. Prof)
2007Graduate School of Science, Kyoto University (Assoc. Prof.)91

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Living matter is intrinsically hierarchic. In molecularscales, atomic structures made revolution in the field,
but a gap between atomic resolution informationand cellular biological interest is still large.
Experimentally, large fluctuation involved in hugebiomolecular systems makes the structural analysis
harder. Thus, molecular dynamics simulations canaugment experiments hopefully filling a gap between
high resolution structural information and cellularbiology. Yet, atomic simulation cannot easily deal
with long time dynamics. To this end, for longtime simulations of huge systems, coarse-graining isuseful.
Given these situation, our purpose is, on the basisof atomic resolution structural data and simulation
technology, to1) systematically construct CG models,
2) develop a software CafeMol that implements the
CG models, and,
3) apply them to various biological phenomena atmolecular scale important in cellular biology.

2. Applications to various biologicalphenomena
In parallel to the development of CG models andthe software CafeMol, we have been applying
them to various biological phenomena in molecularscale important in cellular biology. In particular,
applications include a few molecular machines, suchas multidrug transporter AcrB and molecular motor
kinesin, and simulations of nuclear processes, suchas nucleosome dynamics and transcription factor p53
search dynamics.
References
1. H. Kenzaki, N. Koga, N. Hori, R. Kanada, W.
Li, K. Okazaki, X. Yao, & S. Takada, CafeMol:
A coarse-grained biomolecular simulator for
simulating proteins at work, J Chem TheoComp , 7: 1979-1989, 2011.
2. W. Li, P. G. Wolynes, & S. Takada, Frustration,
specific sequence dependence, and nonlinearity inlarge-amplitude fluctuations of allosteric proteins,
Proc Nat Acad Sci USA , 108: 3504-3509, 2011.
3. W. Li, T. Terakawa, W. Wang, & S. Takada, Energy
landscape and multiroute folding of topologicallycomplex proteins adenylate kinase and 2ouf-knot,
Proc Nat Acad Sci USA , 2012 in press.
4. N. Hori & S. Takada, Coarse-Grained StructureBased Model for RNA-Protein Complexes
Developed by Fluctuation Matching, J Chem
Theo Comp , 8: 3383-3394, 2012.
5. X. Yao, H. Kenzaki, S. Murakami, & S. Takada,
Drug export and allosteric coupling in a multidrug
transporter revealed by molecular simulations,
Nature Comm , 1:117 (8pages), 2010.
6. T. Terakawa & S. Takada, Multiscale ensemble
modeling of intrinsically disordered proteins: p53
N-terminal domain, Biophys J , 101: 1450, 2011.
7. T. Terakawa, H. Kenzaki, & S. Takada, p53
searches on DNA by rotation-uncoupled sliding
at C-terminal tails and restricted hopping of core
domains , J. Am Chem. Soc. , 134: 14555, 2012.

1.2. Development of CG models and asoftware CafeMol
In 2007, we started coding CafeMol, a genericsoftware for CG simulations of biomolecular systems.
Dr. Kenzaki primarily coded the major part of thesoftware, and others contributed to write some
portions. CafeMol is written in Fortran 90 and isparallelized with openMP and MPI. On K computer,
we tested it showing a good scaling up to 98000 coresand 33% efficiency in a single node.
CafeMol uses CG models that has a one- beadper- amino acid resolution for proteins, and threebead- per- nucleotide resolution for nucleic acids.
Lipid model is still premature. We developed somemultiscale protocols to derive CG parameters from the
all-atom force field and experimental structural data.

92

93

Summary– By CG models, it simulate long-time dynamics of huge biomolecular systems.

4th Biosupercomputing Symposium• Approach• Method– Clasical molecular dynamics with CG model• Simulations– Time propagation based on Langevin dynamics and others• Parallel programming– Neighbor list, replica exchange• Computer language and library– Fortran90, MPI, OpenMP• Code release– Freely downloadable from http://www.cafemol.org• Target size
– Simulations that correspond to milliseconds and longer.

•

CafeMol

4th Biosupercomputing Symposium

Shoji Takada
Dept Biophysics, Grad School of Science, Kyoto University

3

1

Molecular machines and nuclear processes
studied by coarse-grained molecular simulations•

PurposeOn the basis of atomic resolution structural data andsimulation technology, we１）systematically construct CG models,
２）develop a software CafeMol that implements the CG
models, and,
３）apply them to various biological phenomena atmolecular scale important in cellular biology.

•

AtomicCG

4th Biosupercomputing Symposium– Applications
Multidrug transporter, kinesin,
nucleosome, transcription factor (p53 etc)
– Experimental verification
Multidrug transporter, collaboration with S. Murakami (Tokyo Inst Tech)
Kinesin, collaboration with M. Tomishige (Univ Tokyo)

Current status of CafeMol– Completed
Implementation of CG models for protein, NA, lipid
Parallel computation up to 98000 cores on K
Single node efficiency up to 33% on K
– To doAccurate derivation of and more accuracy in CG models →Multiscale algorithm for CG parameter derivation from atomic level

CafeMol：Current status

4th Biosupercomputing Symposium

Background
Living matter is hierarchic
In molecular scales, atomic structures made revolution
in the field, but a gap between atomic resolution andcellular scale is still large. Large fluctuation makes the
structural analysis harder.
For long-time simulations of huge systems, coarsegraining is useful.

•

Background & Purpose

4

2

a.a.

4th Biosupercomputing Symposium (December 3 - 5, 2012)

2

IJ
nloc

[5(r

I

IJ
0

I
dih

IJ
0

/ r ) − 6(rIJ 12∑ ε (C / r )94“Entrance”

Murakami 2006

QuickTime™ and a
Cinepak decompressor
are needed to see this picture.

Access[[[ligandTransmembrane
domainPorterdomainTolCdocking
domain)rmsf bet. AA and CG

Fluctuation matching

pairwise contact energies

Energy decomposition

I > J +3

IJ 12

2w 2

(r IJ − r0IJ ) 2

non − nativeexp(−/ r ) ]+εIJ 10

J + 2≤ I ≤ J + 3

∑
IJ
loc

Local contact

Yao et al, Nature Comm 2010

Application 1) Multidrug transporter

Flexible local potential(Statistical potential)

I > J +3+ ∑εnativeI

Ia+ ∑V (θ I ) + ∑V (ϕ I ) +I

V = ∑ k b (r I − r0I )AICG2 model

Yao

Gohlke et al., 2003,
Chu and Voth, 2006
Li et al., 2010,

AICG1; Li, Wolynes, & Takada PNAS 2011
AICG2;Li, Terakawa, Wang, & Takada PNAS 2012

Deriving CG para from AA via a multiscale method

Parallel performance
REMD @RICC
parallel 99.995%
DNA：
Protein-DNA
Structure-based+ P-Chemcomplex

1-node efficiency 33％

InnermembraneTunnel 2A

B

Tunnel 1

Tunnel 2

Tunnel 1 , Long path from membrane surface
Tunnel 2 , Short path from periplasm

Drug uptake pathways

Protein:
Structure-based, switching

Source available for ver 2.0(proteins, DNA, & RNA), atwww.cafemol.org

Running on Kcomputer

Fortran 90,
MPI & openMP

Also for nucleic acids, membrane, and their complex

Simulating “protein at work”

CafeMol, software for CGMD

1000001000000
1 peta

1

100

10000

E

Tunnel 3

1000000

J. Comp.
Theo. Chem.
2011

Kenzaki#core
Membrane：
P-Chem10
10giga100

1000
1 tera

10000

FLOPS

Atomic interaction based CG model

4th Biosupercomputing Symposium (December 3 - 5, 2012)

95

Application 3)
Nucleosome dynamics

CG model of DNA
de Pablo et al

Assembly & disassembly of
nucleosome

Towards modeling of
Chromatin structure

Kenzaki

We founddrug-dependent uptake pathways

4th Biosupercomputing Symposium• 273 ligand molecules in a
160x160x300 (Å) box• 300K Langevin• 100 trajectories• Switching
ABE  BEA  EAB

Dynamic simulation of drug uptake and export
軸索[神経細胞]

微小管
キネシン
ダイニン

KanadaBPJ 2011
JACS 2012

Terakawa

In collaboration withexperimentalists (Tomishige, Kikkawa)

CG MD simulations together withhigh-resolution modeling of K-MT complex

Application 4) full-length p53

中心体

Application 2) Kinesin

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Biomolecular simulations under cellularenvironment
Yuji Sugita
Chief Scientist, RIKENProfile:Education and degree:
Ph. D from Kyoto University (March, 1998)
Professional record:
Postdoctoral fellow, RIKEN, April – August, 1998
Research associate, Institute for Molecular Science (IMS), September, 1998 – March, 2002
Lecturer, Institute of Molecular and Cellular Biosciences, the University of Tokyo, April, 2002 – March, 2007
Associate Chief Scientist, RIKEN, April, 2007 – March, 2012
Chief Scientist, RIKEN, April, 2012 - Present
Team Leader (concurrent post), RIKEN Advanced Institute for Computational Science, October, 2010 - Present
Team Leader (concurrent post), RIKEN Quantitative Biology Center, April, 2011 - Present

97

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
At the cytoplasm in cells, 70% of the volume fraction
is occupied by water molecules and in the rest of
the volume, a large number of proteins, RNA, orother metabolites exist. This environment is often
called as ‘macromolecular crowding’ in cell and has
been paid a lot of attention by many theoretical andexperimental researchers. One of the key questionsin this issue is whether protein is stabilized ordestabilized in such crowded environment. Mostof the previous studies suggest that protein is morestabilized in crowded environments rather than that indilute solution. However, recent experimental studies
suggest that protein stability significantly dependent
on the crowding environments. The hydrogenexchange experiments in NMR spectroscopy showsthat protein stability significantly depends on thetypes of crowding agents: In another experiment, incell NMR spectroscopy, Inomata et al. found thatubiquitin is destabilized inside of the cell [Inomata
et al. Nature (2009)]. These experiments suggest that
protein-protein interaction in crowded environment isessentially important and is not fully examined basedon the previous theoretical studies in which only
volume exclusion effect is considered as the major
crowding effect. We therefore performed atomistic
molecular dynamics (MD) simulations of proteins
under crowded conditions in explicit water.

specific protein-protein interactions between CI2 andlysozyme. CI2 interacts much less with BSA crowders
corresponding to a weak effect of crowding. Energetic
analysis suggests an overall favorable crowding freeenergy in the presence of lysozyme while weaker
interactions with BSA appear to be unfavorable.
3. Hydration and protein-protein
interactions in crowded conditions
In the next applications of crowding simulations,
we performed explicit solvent MD simulations of aseries of protein G and protein G/villin systems at
different protein concentrations. In this study, we
focused on hydration structure and protein stability.
Hydration structure was analyzed in terms of radial
distribution functions, three-dimensional hydrationsites, and preservation of tetrahedral coordination.
Analysis of hydration dynamics focused on selfdiffusion rates and dielectric constants as a functionof crowding. The results show significant changes inboth structure and dynamics of water under highly
crowded conditions. The structure of water is altered
mostly beyond the first solvation shell. Diffusion
rates and dielectric constants are significantly reduced
following linear trends as a function of crowdingreflecting highly constrained water in crowded
environments. The reduced dynamics of diffusion
is expected to be strongly related to hydrodynamic
properties of crowded cellular environments while thereduced dielectric constant under crowded conditions
has implications for the stability of biomolecules incrowded environments. The results from this study
suggest a prescription for modeling solvation insimulations of cellular environments.

2. Protein stability changes due to different
protein crowders
In the first applications of crowding simulations,
we examined the effect of cellular crowding byperforming MD simulations of chymotrypsin inhibitor
2 (CI2) in the presence of either lysozyme or bovine
serum albumin (BSA) crowder molecules. This is thesame protein and crowder proteins in which Pielak
et al recently studied by the hydrogen-exchange
experiment in NMR spectroscopy (Miklos, A. C.;
Sarkar, M.; Wang, Y.; Pielak, G. J. J. Am. Chem.
Soc. 2011, 133, 7116). The simulations confirm adestabilization and significantly slowed diffusionof CI2 in the presence of lysozyme and indicatethat this observation is a result of extensive, non-

4. Summary and perspectives
The results provided MD simulations with explicit
solvent seem to be complementary with experimental
data and show new insight on the crowding effecton proteins in cell. Using K computer, we expect
to simulate larger and more realistic systems of thecytoplasm to examine crowding and confinement
effect on biological molecules.

98

Dilutesolution(single protein)

Our study

Cellularenvironment（parts of cell）

99

folded

4th Biosupercomputing Symposium

Cell

Ultimate goal

unfoldedClassical view of crowding

4th Biosupercomputing Symposium

Crowdedenvironment(multiple proteins)

Current simulations

Simulations on K

Yuji Sugita
RIKEN

Biomolecular simulations under cellularenvironments

3Less stable

More stable

4th Biosupercomputing Symposium

A. Miklos, M. Sarkar, Y. Wang, G. Pielak: JACS (2011) 133, 7116-7120

NMR amide proton exchange experiments

Protein Stability in Crowded Environments

4th Biosupercomputing Symposium• Towards cellular-scale simulations– Hydration properties under crowded environments– Protein stability changes due to non-specific protein-protein interactions– NMR experiments on protein-protein interactions in crowded conditions• Hydration and protein-protein interactions in crowded conditions– NMR amide proton exchange experiments– MD simulations of crowded environments• Protein stability changes due to different protein crowders– How macromolecular crowding affects protein stability and function?

• Introduction

Topics42

4th Biosupercomputing Symposium (December 3 - 5, 2012)

100

A. Miklos, M. Sarkar, Y. Wang, G. Pielak: JACS (2011) 133, 7116-7120

Feig, Sugita: J. Phys. Chem B (2012) 116, 599-605

increased/decreased
fluctuations in MD

Residue fluctuations vs stability~120 ns MD~250 ns MD~160 ns MD

RMSF differences between crowded and dilute solvent

CI2 + 8 BSAs
835K atoms
253K H2O/136 Na+
104 g BSA/L solution
7% vol. fraction

CI2 + 8 lysozymes
184K atoms
56K H2O/64 Cl108 g lysozyme/L
7% vol fraction

CI2
18K atoms
6K H2O
infinite dilution

Feig, Sugita: J. Phys. Chem B (2012) 116, 599-605

MD simulations of crowded environments

7

5

4th Biosupercomputing Symposium

Feig, Sugita: J. Phys. Chem B (2012) 116, 599-605

lysozyme

-21.3

⊗⊗H

[kcal/mol]

Lysozyme -194.0

BSA

∑ =80

∑=65

∑=65

0.5

0.3

∑=65

4.3

4.2

-3.0

6.1

g

⊗⊗G80/65 ⊗⊗Gcrowdin⊗⊗Gbinding =⊗⊗H – T⊗⊗S + ⊗⊗Gsolv

Feig, Sugita: J. Phys. Chem B (2012) 116, 599-605

186.2

22.9

⊗⊗Gsolv -T⊗⊗S⊗⊗G8065

Free Energy of Crowding for CI2 in BSA/Lysozyme

BSA

CI2 interaction with crowded proteins86

4th Biosupercomputing Symposium (December 3 - 5, 2012)

PGVH1

PGVH2

PGVH3

101

1
MSD(t )
6tε=−132.361 f + 77.817--- : linear fitting

4th Biosupercomputing Symposium

Ds =−0.569 f + 0.308 [ Å 2 / ps]---: linear fitting＋: PG0, PG1-4, × : PGVH1-5

t →∞Ds = lim

PGVH5

11

9

NPT(1bar, 298K, 300ns)

PGVH4

PG1-4:
8 protein G
Crowding

Universal behaviors of diffusion coefficient anddielectric constant

4th Biosupercomputing Symposium

Harada, Sugita, Feig, J. Am. Chem. Soc. 2012, 134, 4842-4849

PGVH1-5:
4 protein G + 8 villin
Crowding

PG0
PG0:
1 Protein G
DiluteCrowding simulations of two different protein systems

crowdingdilutecrowding(Villin:PGVH1-5)

crowding(ProteinG:PGVH1-5)

10

Decreases in diffusion coefficient and dielectric constant were observed incrowded systems– Reduction of diffusion coefficient→ Modulations of rates of biochemical reactions– Reduction of dielectric constant→ Thermodynamic properties(directly protein-protein interactions ↑)4th Biosupercomputing SymposiumThis work has been carried out in collaboration with Prof. Michael Feig (Michigan
State University) and Dr. Ryuhei Harada (RIKEN AICS).
12

Crowding strongly altered hydration structure with protein volume fraction over
30%vol, reducing the number of bulk water keeping that of the first solvation shell.

Summary and perspectives

4th Biosupercomputing Symposium

Acknowledgement to collaborators

•

•

t'

(ProteinG:PG0-4)| r (t '+ t ) − r (t ') |2★ In crowded environment, diffusion of water and protein significantly slow down

crowding(Water:PGVH1-5)(Water:PG0-4)
crowding

MSD(t=)

Dynamics of water and proteins

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Supercomputing for Next-Generation Cancer
Research
Satoru Miyano
Professor of Human Genome Center, The Institute of Medical Science,
The University of Tokyo

Profile:Satoru Miyano, Ph.D., is a Professor of Human Genome Center, The Institute of Medical Science, The Universityof Tokyo. He received the B.S., M.S. and Ph.D. all in Mathematics from Kyushu University, Japan, in 1977, 1979
and 1984, respectively. He joined Human Genome Center in 1996. His research mission is to create computational
strategy for systems biology and medicine towards translational bioinformatics.

103

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
Computational systems biology is gradually seeping
in cancer research. We present a challenge foruncovering systems in cancer by supercomputer.
SiGN-L1 (NetworkProfiler) (based on L1regularizaon) is a method that will exhibit how gene
networks vary from patient to patient according
to a modulator, which is any score representing
characteristics of cells. We defined an EMT(epithelial-mesenchymal transition) modulator andanalyzed gene expression profiles of 762 cancer celllines. The computation took 3 weeks on 1024 CPU
cores. Network analysis unraveled global changes of
networks with 13,508 genes of different EMT levels.
By focusing on E-cadherin, 24 genes were predicted
as its regulator, of which 12 have been reported in
the literature. A novel EMT regulator KLF5 was also
discovered in this study. We also analyzed Erlotinib
resistant networks using 160 NSCLCs with GI50 asa modulator. Hubness analysis exhibited that NKX21/TTF-1 is the key gene for Erlotinib resistance in
NSCLCs. Our microRNA/mRNA gene network
analysis with Bayesian network method called SiGN
also revealed subnetworks with hub genes (including
NKX2-1/TTF-1) that may switch cancer survival.
The supercomputer was also applied for modelingdynamics in cancer cells from time-course gene
expression profiles and revealed dynamic networkchanges against anti-cancer drugs and network
differences between drug-sensitive and drug-resistantcancer cells. For dynamic system modeling, we
devised a state space model (SSM) with dimension
reduction method for reverse-engineering gene
networks from time-course data, with which we canview their dynamic changes over time by simulation.
We succeeded in computing a gene network withprediction ability focused on 1500 genes from data of
about 20 time-points. We applied this SSM model tohuman normal lung cell treated with (case)/without(control) Gefitinib, and we identified genes under
differential regulations between case and control. Thissignature of genes was used to predict prognosis forlung cancer patients and showed a good performancefor survival prediction.
References
1. Shimamura T, Imoto S, Shimada Y, Hosono Y,
Niida A, Nagasaki M, Yamaguchi R, Takahashi
T, Miyano S. A novel network profiling analysis
reveals system changes in epithelial-mesenchymaltransition. PLoS ONE. 6(6): e20804, 2011.
2. Tamada Y, Imoto S, Araki H, Nagasaki M, Print
C, Charnock-Jones DS, Miyano S. Estimating
genome-wide gene networks using nonparametric
Bayesian network models on massively parallel
computers. IEEE/ACM Transactions on
Computational Biology and Bioinformatics. 8(3):683 - 697, 2011.
3. Tamada Y, Yamaguchi R, Imoto S, Hirose O,
Yoshida R, Nagasaki M, Miyano S. SiGN-SSM:
open source parallel software for estimating gene
networks with state space models. Bioinformatics.
27: 1172-1173, 2011.
4. Yamauchi M, Yamaguchi R, Nakata A, Kohno T,
Nagasaki M, Shimamura T, Imoto S, Saito A, Ueno
K, Hatanaka Y, Yoshida R, Higuchi T, Nomura
M, Beer DG, Yokota J, Miyano S, Gotoh N.
Epidermal growth factor receptor tyrosine kinase
defines critical prognostic genes of stage I lungadenocarcinoma. PLoS ONE. 7(9): e43923, 2012.
5. Yamaguchi R, Imoto S, Yamauchi M, Nagasaki
M, Yoshida R, Shimamura, T, Hatanaka Y, Ueno
K, Higuchi T, Gotoh N, Miyano S. Predicting
differences in gene regulatory systems by state
space models. Genome Informatics. 21:101-113,
2008.

104

105

1/2 → Cancer
1/3 → Death

4th Biosupercomputing Symposium

Q: How is “My” gene-gene causal
network structured in My Cancer?

Modulator-dependent Graphical Model
Using Structural EQ

NetworkProfiler

SiNG-L1

4th Biosupercomputing Symposium“The Only Flower in the World”- Genetic Variations -Satoru MIYANOHuman Genome Center
The Institute of Medical Science, The University of Tokyo

Supercomputing for Next-Generation
Cancer Research

4th Biosupercomputing Symposium

Genetic variations/mutations inherited from parents (DNA)
Mutations accumulated in somatic cells (Cancer Genome)
Modifications of DNA and Histone (Epigenome)

Modulator A

Role of MODULATOR

4th Biosupercomputing Symposium

Gene 1

Change can be found if we define some
MODULATOR and look data along it•••

30-800 Mutatations on Genomes

Gene Networks

These variations affect malignancy of cancer, therapeutic
response/effect, and adverse reaction by anticancer drugs

2

4th Biosupercomputing Symposium (December 3 - 5, 2012)

106

High

http://gantoku3.umin.jp/topics/saito.html

4th Biosupercomputing SymposiumSome gene inducing EMT
have been known.
However, its mechanism isstill not well understood.

？

Signature-based hidden modulator
extraction algorithm
1. Selected 122 genes labeled “EMTUP”, “EMT DN”, “JECHLINGER EMTUP”, and “JECHLINGER EMT
DN”from Molecular Signatures
Database v2.5 ( [6];http://www.broadinstitute.org/gsea/msi
gdb/index.jsp).
2. Then, narrowed the set to 50
functionally coherent genes with p <10−5 by using the extraction of
expression module (EEM) .
3. Computed the first principal
component of these 50 genes as
hidden values of the EMT-related
modulator

4th Biosupercomputing Symposium

We selected coherent 50 genes from 122 EMT signature genes to define
the modulator for EMT (EEM, Niida et al., Bioinformatics, 2009)

Network Estimation 4th Biosupercomputing Symposium

Not always similar each otherα-th sample• Key developmental remodeling program, where cells
alternate between Epithelial-like and Mesenchymal-likephenotypes• Relate to tumor grade, metastasis, drug resistance

Gene Network

Low

modulator

Difficulty: Data for α-th sample is only one!!
Idea for solution: Data integration by sample weighting

Estimate coefficients for α-th sample βjkα

Modulator for EMT

4th Biosupercomputing Symposium

Samples

Genes

Epithelial-Mesenchymal Transition (EMT)

Co-expression Network

Node: gene transcript
Edge: conditional dependence(not equal to correlation)

Structural Equations

Modeling with Structural Equation

4th Biosupercomputing Symposium (December 3 - 5, 2012)

low

Low

107

4th Biosupercomputing Symposium

226 Patient Samples of NSCLC with survival data
Collaboration with National Cancer Center ResearchInstitute, Japan

Define “relapse risk score” by using relapse
makers; MODULATOR

Relapse Risk

4th Biosupercomputing Symposium

762 gene networks with 13,508 genes exhibited how networks
structurally change from low to high EMT modulator score

OutputInput• Transcriptome data of 762 cancer cells focused on 13,508 genes
13,006 mRNAs + 581 miRNAs• EMT Modulator

Run Supercomputer for 3 Months!

11

high

High

c

b

α-tubulin

VimentinE-cadherin

DAPI

si NC
E-cadherin

Mock

si NC

Replase risk

Compute each patient network

4th Biosupercomputing Symposium

Cancer Res 62(6): 1613-8, 2002

Oncogene 24(14):2375-85, 2005

Cancer Res 69(7):2783-91, 2009

Cancer Res 70(5):2115-25, 2010

PNAS 104(9): 3432-7, 2007

PLoS ONE 6(6): e20804, 2011

Nat Cell Biol 10(5): 593-601, 2008
Cancer Res 70(5):2115-25, 2010
J Cell Sci 122(Pt 7): 1014-24, 2009

Nat Cell Biol 10(5): 593-601, 2008
J Biol Chem 285(22):16854-63, 2010

Mol Cell 7(6): 1267-78, 2001

Cancer Res 72(9):2440-53, 2012asi KLF5-#1

TGF-β

si KLF5-#1

4th Biosupercomputing Symposium

12H

Shimamura T, Imoto S, Shimada Y, Hosono Y,
Niida A, Nagasaki M, Yamaguchi R, Takahashi
T, Miyano S. A novel network profiling analysis
reveals system changes in epithelialmesenchymal transition. PLoS ONE. 6(6):e20804, 2011.

Network of the lowest risk patient
Network of the highest risk patient
CTGF (Connective Tissue Growth Factor)
CTGF has been shown to be associated with tumor development and progression. CTGF may
regulate cancer cell migration, invasion, angiogenesis, and anoikis.

L

These 12 highlighted of genes are known to have regulatory effect among EMT
in the literature.

•

Nat Cell Biol 10(5): 593-601, 2008

Krueppel-like factor 5 (KLF5) from a list of the remaining candidate regulatorsand conducted in vitro validation experiments. As a result, we found thatknockdown of KLF5 by siRNA significantly decreased E-cadherin expression andinduced morphological changes characteristic of EMT (Validated by Takahashi
Lab, Nagoya University)•

EMT Regulators Predicted by Supercomputer
A New Regulator KLF5 of EMT

4th Biosupercomputing Symposium (December 3 - 5, 2012)

108

4th Biosupercomputing Symposium

Supercomputer

4th Biosupercomputing SymposiumHigh in poor prognosis
Low in poor prognosis

Expression ratio

3000

4000

Computatioanl Systems Biology

‡

‡

‡‡
‡

‡

‡

TTF-1/
NKX2-1‡High in TTF-1/NKX2-1-introduced
Low in TTF-1/NKX2-1-introduced: N.A. (miRNA) ‡ : Positive binding in E19.5 mouse lung
ChIP-chip assay data (Tagne et al., PLoS ONE, 2012)

Expression ratio (TTF-1-intorudced / Vector control)

‡

‡

‡

‡

After introducing TTF1/NKX2-1

‡

ROR1

NKX2-1/TITF1/TTF-1-Induced
ROR1 Is Required to Sustain EGFR
Survival Signaling in Lung
Adenocarcinoma. Cancer Cell
2012 by Takahashi Lab (Nagoya U)

TTF-1/
NKX2-1

Inferred subnetwork of TTF-1/NKX2-1
with good prognosis
ROR1

2000

Days after surgery

NKX2-1 and ROR1 appear clearly on the network andshows the consistency and new insights

Gene Network ofLung Adenocarcinoma

SiGN• Nonlinear Bayesian
Network Estimation
Softwarehttp://sign.hgc.jp/signbn/

Gene Network of 400 mRNA+32 microRNA
Computed from 124 Patient Expression Profiles

Cancer Biology

1000

14

subnetworks and their hub genes
which are significantly corelated withsurvival /death

4th Biosupercomputing Symposium0

P < 0.0001

High
Low

miR30c subnetwork

Survival Data

Searching for subnetworks with smaller p-valuesin Kaplan-Meyer diagram survival prediction andhub genes

0.00

0.25

0.50

0.75

1.00

Submnetwork B

Hub 遺伝子Subnetwork A

Functional
Analysis of
Hub Genes

4th Biosupercomputing Symposium

124 samples oflungadenocarcinoma
patients

Nonlinear Bayesian
Network

Finding genes that “switch”
good/poor prognosis in lung cancer

Bayesian Network and Nonparametric
Regression applied to a large set of lungcancer patient gene expression profiles

SiGN-BNmRNA and miRNA expression profile

Systems Biology approach unraveled key hub
genes that switch good/poor prognosis inlung adenocarcinoma

Overall survival

Collaboration withTakashi Takahashi (Nagoya U)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

109

4th Biosupercomputing Symposium

Survival prediction using the
139 genes is very successful
for stage I lungadenocarcinoma

Simulation↓
Identified 139 genes

・Simulation model based on
Gefitinib + EGF time-course
data・State Space Model withdimension reduction

・SAEC (lung epithelial cell)・Gefitinib +/・EGF stimulation
・19 time points in 48hYamauchi M*, Yamaguchi R*, Nakata A*, Kohno T*, Nagasaki M, Shimamura T, Imoto S, Saito A, Ueno
K, Hatanaka Y, Yoshida R, Higuchi T, Nomura M, Beer DG, Yokota J, Miyano S, Gotoh N. (2012)
Epidermal Growth Factor Receptor Tyrosine Kinase Defines Critical Prognostic Genes of Stage I Lung
Adenocarcinoma. PLoS ONE 7(9): e43923, 2012.

Duke U data

NCI data

19

The power of survival prediction for the NCI
validation data set, the Duke data set and NCC-Tokyo
data set using the 139 genes.

2

4th Biosupercomputing Symposium

Extracting the
Essense of System
Dynamics

００００）State Space Model (SSM) Approach Captured
Gefitinib Response Systems

xn Fxn −1 + vn
Gene =
Expression =ynSampleHx#n（～20）+ wn
Profiles（～4th Biosupercomputing SymposiumState Space
Model
The
Evil&of
Dimension
Dimensions
Reduction

SSM

State Space Model + Dimension Reduction

Gene

Discovery of hub genes that “switch”
survival/death for lung cancer

State Space Model for Inferring
Transcriptional Module Networks from TimeCourse Gene Expression Data

SiGN-SSM=
yn Hxn + wn ∈ R p

N obs << p ≈ 103∫

p(Y , X | θ )dX

xn p xn | y1 ,, yn −1 ,θ )dxn

∫ (R −1/ 2 ( yn − wn ) =ΨR −1/ 2 ( yn −1 − wn −1 ) + R −1/ 2 Hvn 18

Module-Based Gene Network Estimation:yn|n −1 = H

Gene Expression Prediction:

H T R −1 H = Λ = diag(λ1 ,, λk )L(θ ) =

dim(xn ) =k < dim( yn ) =p}

System Estimation with Dimension
Reduction:

Y = y1 ,, y N obs{

High-Dimensional Short Time
Series Data:

System Model
Observation Model=xn Fxn −1 + vn ∈ R k

State Space Model:

Modeling Cancer Systems with
State Space Model

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Large-scale protein-protein interaction
network prediction by an exhaustive rigid
docking system MEGADOCK
Yutaka Akiyama
Professor of Department of Computer Science, Graduate School ofInformation Science and Engineering, Tokyo Institute of Technology

Profile:Yutaka Akiyama is a professor of the Department of Computer Science, Graduate School of Information Scienceand Engineering, and also the Director of the Education Academy of Computational Life Sciences (ACLS), Tokyo
Institute of Technology. He is currently serving as the President of Initiative for Parallel Bioinformatics (IPAB), a
board member of BioSuperComputing Research Community (BSCRC), and a board member of Japanese Society
for Bioinformatics (JSBi).
He received B.E.(1984), M.E.(1986), and Dr. Eng.(1990) in Electrical Engineering from Keio University. His
doctoral thesis was about theoretical studies and VLSI implementation of a novel neural network model called“Gaussian Machines”. On that study, he received Best Paper Award for Young Researcher of IPSJ National
Convention (1988). After his enthusiastic activity to initiate world’s earliest bioinformatics WWW services(1995) on GenomeNet Japan when he was an associate professor at Kyoto University, he started his own
computational biology research group (1996) as a part of the Real World Computing national project with using
massively parallel computers and pioneering large-scale PC clusters. In 2001, he became the founding Directorof Computational Biology Research Center (CBRC), the biggest government-supported bioinformatics research
center in Japan under Ministry of Economy, Trade and Industry.
His research interest covers large-scale processing and acceleration techniques for computational biology andbioinformatics, including ultra-fast metagenome sequence analysis, protein-protein interaction prediction, virtual
screening of drug compounds, machine learning for pharmacokinetics prediction, and software system for mass
spectrometry analysis.
111

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
Protein-protein interaction (PPI) plays a core role incell functions. Massively parallel supercomputing
systems have been actively developed recently, thatenables us to solve large-scale biological problems
such as PPI network prediction based on tertiary
structures.
To challenge interactome level large-scale analysis
by fully utilizing protein tertiary structures, wehave proposed a large-scale PPI prediction system“MEGADOCK” based on exhaustive protein docking
and post-docking analysis [1][2]. We input protein
structure data to the system and get predictions ofpossible interacting pairs among them.
We have already applied our system to 44×44 (subset
of protein docking benchmark 2.0) and 89×89(structures of bacterial chemotaxis proteins) scale
analyses. In real biology problem, such as searching
drug induced pathway of EGFR signaling, about
200 proteins should be examined. In our preliminary
survey on the EGFR pathway and related proteinsdata, we found about 2000 structures corresponding
to these proteins. Therefore, the PPI network
prediction system needs to handle about 2000×2000
combinations of protein structures.

Speedup4

LinearObserved
1.97

2

1

3.76

1.00
6144

12288

#Nodes24576

Parallel scaling of MEGADOCK on K computer.
structure data including multiple structures for eachprotein species), human apoptosis pathway (57proteins, 158 structures) and human EGFR signaltransduction pathway (49 proteins, 497 structures).
The F-measure value was 0.44 when applied tochemotaxis pathway, 0.28 when applied to apoptosis
pathway and 0.36 when applied to human EGFR
pathway.
The proposed approach to computational PPI
detection is a promising methodology for mediating
between structural studies and systems biology byutilizing cumulative protein structure data for pathway
analysis.
References
1. Matsuzaki,Y.,Matsuzaki, Y., Sato T, Akiyama Y.
In silico screening of protein-protein interactions
with all-to-all rigid docking and clustering: anapplication to pathway analysis. J. Bioinform.
Comput. Biol ., 7:991-1012, 2009.
2. Ohue M, Matsuzaki Y, Akiyama Y. Dockingcalculation-based Method for Predicting ProteinRNA Interactions. Genome Inform ., 25:25-39,
2011.

2. Implementation
To solve such large-scale problems, highly efficientcomputing system is necessary. We implemented
MEGADOCK by using hybrid parallelization
where each docking processes are calculated inparallel within one node and number of docking
jobs are distributed among nodes. We also designed
MEGADOCK using a simple score model to reduce
calculations required for protein docking.

Acknowledgements
This work was supported in part by a Grant-inAid for Research and Development of The NextGeneration Integrated Life Simulation Software from
the Ministry of Education, Culture, Sports, Scienceand Technology of Japan (MEXT). Part of the result
was obtained by early access to the K computer at
the RIKEN Advanced Institute for Computational
Science.

3. Results and Conclusion
MEGADOCK showed almost linear scaling up to
24,576 nodes on K computer (RIKEN Advanced
Institute of Computer Science, Japan).
As an application to real biology pathway, PPI
prediction using MEGADOCK was performed on
the reconstruction of a canonical signal transduction
pathway of bacterial chemotaxis (13 proteins, 89
112

113

PPI prediction based on
tertiary structures

P1000

Where?

?

?

1:1 docking?

?

Our problem

Do not interact

Interact• Developing an ultra-fast docking tool MEGADOCK
feasible for exhaustive docking of million-order pairs• Questions asked by protein docking– 1:1 docking: How given proteins interact?
– Our challenge by exhaustive docking:
Which protein pairs interact?

P1000

P1
P2
P3
P4
P5
…

P1 P2 P3 P4 ….

Yutaka AKIYAMATokyo Institute of Technology

Large-scale protein-protein interaction
network prediction by an exhaustive
rigid docking system MEGADOCK

3

P2

P3

P4

P5

2+H
-45 -45 -45 -45 +iφ
5+H 2+H 1+H-45 -45 +iφ +iφ +iφ5+H 2+H 1+H-45 -45 +iφ +iφ +iφ2+H-45 -45 -45 -45 +iφ
2+H-45 -45 -45 -45 +iφ

3+H+iφ -453+H+iφ -453+H+iφ -453+H+iφ -452+H+iφ -45

Receptor

1+H 2+H 3+H 3+H 3+H 2+H 1+H+iφ +iφ +iφ +iφ +iφ +iφ +iφ2+H-45 -45 -45 -45 +iφ2+H+iφ -45

1+H 2+H 3+H 3+H 3+H 2+H 1+H+iφ +iφ +iφ +iφ +iφ +iφ +iφCompress three terms into
one complex number:• Shape complementarity• Hydrophobic interaction• Electrostatic interaction

1+iq 1+iq

Ligand

4th Biosupercomputing Symposium

Ohue M, Matsuzaki Y, Ishida T, Akiyama Y.
Lecture Note in Bioinformatics 7632,
178-187, Springer, 2012.

1+iq 1+iq

1+iq 1+iq 1+iq 1+iq 1+iq

1+iq 1+iq 1+iq 1+iq 1+iq

4

Convolution can becalculated fast by FFT(Katchalski-Katzirmodel)

2

Binding partners

A compact score function of
MEGADOCK

Objective:• Search novel PPIs in lung cancer pathway
proteins and cancer drug related proteins• The first large-scale PPI prediction challenge(2000 x 2000 = 4 million scale)• To contribute finding target proteinsof cancer drugs
th Biosupercomputing
4EGFR
signalingSymposium
network

PPI network prediction

P1

P1

Exhaustive docking & post-docking analysis

PPI network prediction by
exhaustive docking

4th Biosupercomputing Symposium (December 3 - 5, 2012)

1

2

4

3.88

#Threads

1.98

LinearObserved

8

7.33

1

2

4

12288
#Nodes

5

24576

3.76

…

7

An exhausting docking of 220 proteins
Size of FFT: 140

6144

1.00

1.97

LinearObserved

MPI parallelization

4th Biosupercomputing Symposium

One docking of FFT size: 140 proteins1

Save high scoring docking poses

Score calculation (Inverse FFT)

Convolution

Ligand FFT

Ligand voxelization

Loop for Ligand rotations

Receptor FFT

…

R3

R2

R1

R44

L1 L2 L3

・・・ L44

1

0

0.2

0.4

0.6

0.8

0

0.2

0.8

MEGADOCK
random1β-actin and profilin(PDB id : 2BTF)

6

8

Predicted docking pose
Crystal structure

Trypsinogen protein andtrypsin inhibitor(PDB id : 1CGI)PPI prediction by MEGADOCK
achieved better than random
performance on general benchmark
dataset (monomer pair from proteinprotein docking benchmark 2.0,
Mintseris et al, Proteins, 2005.)

0.6

FP rate

0.4

Matsuzaki et al., J Bioinform Comput Biol, 7:991-1012 (2009).
Ohue et al., submitted.

F-measure : 0.42

Binding partner prediction from
44x44=1936 dockings and postprocessing based on highest docking
score (Diagonal: interacting pairs)

No
PPI

PPI

4th Biosupercomputing Symposium

PPI prediction using generalbenchmark data

Ohue M., et al., Lecture Note in
Bioinformatics 7632,
178-187, Springer, 2012.

(single node, single thread)(Benchmark 4.0, 176 bound complexes)

Docking accuracy wascomparable to conventional tools

Performance of protein docking
by MEGADOCK

8.8 times speedup
compared to ZDOCK 3.0

・
・
・

2

4

8

LigandReceptor voxelization

Thread parallelization

Scalability on K computer

Worker
Node

WorkerNodeThread parallelization

Master
Node

WorkerNode

Receptor

MPI parallelization

Speedup114

Speedup

OpenMPTP rateImplementation withhybrid parallelization

4th Biosupercomputing Symposium (December 3 - 5, 2012)

115





13 proteins(101 structures)Precision 0.29

Recall0.47
F-measure 0.36

hsa05223

44 proteins（497 structures)

59

Negative

Counts are basedon protein species

131

53

Positive

11

747

No
Interaction

Interacting

Prediction

Achieved high PPI prediction performance

497 structures, all-to-all docking＝247,009 structure pairs

9

PPI prediction of about
250 thousand structure
pairs showed comparable
performance to theapplication to bacterial
chemotaxis (10 thousand
pairs).

Application to non-small celllung cancer pathway

4th Biosupercomputing Symposium

Completed large-scale exhaustive docking



CheD
A suggestion of CheY-CheC
interaction by using “Falsepositive” pair CheY-CheD as
a mediator

CheC

Asp

CheY

Matsuzaki et al., J Bioinform Comput Biol, 7:991-1012 (2009).

Acceptable performance was shown on areal biology pathway reconstruction problem.

F-measure : 0.44

：True Positive：False Negative：False Positive

Application to bacterial
chemotaxis pathway

364

88
96

Positive
Negative

Ohue, et al.,Tech. Rep. IPSJ SIG, 2012-BIO-32(14),1-8, 2012.

1105

No
Interaction
Interacting
Prediction

Ozbabacan et al., J Struct Biol, 2012).

PPI prediction by docking without
any other knowledge showedcomparable results to templatebased search of interaction
partners (F-measure 0.30,

F-measure : 0.28

44 proteins（497 structures)

EGFR pathway related tonon-small cell lung cancer

294 proteins（1424 structures)

Searching novel cancer related PPIs

4th Biosupercomputing Symposium

2,000 x 2,000= 4 million
docking

Proteins related to Gefitinib
estimated by Miyano lab, Tokyo
Univ. from microarray analysis

Interaction of cancer pathway
and proteins related to Gefitinib

False Negative

False Positive

True Positive

hsa04210

Ozbabacan et al, 201257 proteins(158 structures)

PPI prediction results (Apoptosis)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Supercomputing accelerates genomic
medicine
Tatsuhiko Tsunoda
Director, Research Group for Medical Informatics
RIKEN Center for Genomic Medicine

Profile:Tatsuhiko Tsunoda, Ph.D. (Medicine) & Ph.D. (Eng.), Group Director
1985-1989: B.S., Department of Physics, Faculty of Science, The University of Tokyo
1989-1991: M.S., Department of Physics (Elementary Particle Physics), The University of Tokyo
1992-1995: Ph.D., Department of Engineering, The University of Tokyo
Degrees: Ph.D. (Medicine) & Ph.D. (Engineering)
1995-1997: Assistant Professor, Graduate School of Engineering, Kyoto University
1997-1998: Research Associate, Human Genome Center, The Institute of Medical Science, The University of Tokyo
1998-2000: Assistant Professor, Human Genome Center, The Institute of Medical Science, The University of Tokyo
2000-present: Laboratory Head, Laboratory for Medical Informatics, RIKEN Center for Genomic Medicine(2000-2008 RIKEN SNP Research Center)2011-present: Director, Research Group for Medical Informatics, RIKEN Center for Genomic Medicine

1174th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
Whole-genome approaches accompanied withsupercomputing have opened new frontiers in medical
research. Genome-wide association studies (GWAS)
exhaustively explore disease-related single nucleotide
polymorphisms (SNPs) and genes in the humangenome. They also provide information for deciding
which drugs and dosages are adequate for individuals– personalized medicine. In addition, high throughput
next-generation sequencing (NGS) enables us to apply
whole-genome sequencing to individual genomes
with high quality. Supercomputing accelerates theprocessing of these huge datasets, and leads tonovel findings in genomic medicine by allowing forcomputationally intensive sophisticated mathematicalmodeling and analysis. These technologies have
revolutionized medical research as well as health care.

an efficient algorithm, ParaHaplo, for handling
haplotypes, combinations of alleles at adjacent locion the chromosome which are transmitted together,as markers [4-7]. Because the number of tests
significantly drops when using haplotype blocks
rather than single SNPs as units, we can expect anincrease of power. To calculate the empirical p-values
for significance, we apply permutation procedures andMarkov Chain Monte Carlo methods, which require
supercomputing power.
Another method for processing multiple markers
simultaneously is based on gene-gene/SNPSNP interaction detection. We have developed anefficient algorithm, ExRAT, which processes all SNP
combinations, a very large number of comparisons,and calculates the empirical p-values using theimportance sampling technique.
Beyond the aforementioned techniques for handling
multiple SNPs on current GWAS platforms, we are
exploring methods using other types of markers
to approach the missing heritability problem. We
are using genotype imputation to explore hidden
SNPs with lower allele frequency, as well as tocombine data from chips with different marker sets.
One promising approach for developing higher
quality reference haplotypes and for exploring
unknown variation is analyzing lower frequency
variations, e.g. single nucleotide variations (SNVs)
and copy number variations (CNVs), through NGS.
To develop our analytical pipeline for NGS data,
we sequenced a single genome at high coverage,
resulting in the first reported Japanese individual's
whole-genome sequence [8]. That work allowed us
to establish methodologies for detecting multiple
types of variations: SNVs, structural variations
including CNVs, and novel sequences. Based on
our methodology, we have constructed a pipeline
NGSanalyzer for analyzing personal genomes [8],cancer genomes [9,10], as well as whole-exome
analysis for common/monogenic diseases.

1. Introduction
In 2002, our center reported the world’s first GWAS
results [1]. In 2004, using gene-based data, I
constructed the world’s first linkage disequilibrium(LD) map, and found exotic patterns of natural
selection on genes [2]. Thereafter, we participated in
the International HapMap project to construct a more
extended LD map and select tagging SNPs, which has
been used for chips/arrays [3]. This resulted in a large
increase in the number of GWAS, further accelerated
by the BioBank Japan project, and the discovery ofmany genes related to common diseases, cancers, anddrug responses.
2. Approaches to missing heritability
problem
We are now facing the missing heritability problem:
current GWAS results are insufficient to explain the
expected heritability of common diseases. One ofapproaches to this problem is to increase the powerof GWAS. To this end, we are now enlarging sample
sizes using disease cohorts and performing metaanalysis through collaborations around the world.
Another approach is to examine multiple markers
simultaneously; typical GWAS looks at only single
common SNPs as markers. We have developed118

4th Biosupercomputing Symposium (December 3 - 5, 2012)

3. Supercomputing for genomic medicine
For analyzing huge GWAS and NGS datasets, we
had to implement massively-parallel pipelines
for ParaHaplo, ExRAT, and NGSanalyzer onsupercomputer systems. For each algorithm, we
have developed an efficient pipeline that can be runon the K computer, which allows for the processingof huge datasets of whole genomes for many cases/
controls within practical time: ~ days. In summary,statistical genetics and supercomputing enables us toanalyze many individuals’ whole genome sequenceswith high accuracy, speed, coverage, and precisenessfor applying specific therapies with each individual,
allowing the prevention of disease in an individual by
prediction – personalized medicine.

6. K.Misawa and N.Kamatani. ParaHaplo 2.0: aprogram package for haplotype-estimation andhaplotype-based whole-genome association study
using parallel computing. Source Code for Biology
and Medicine , 5, 5 (2010).
7. K.Misawa and N.Kamatani. ParaHaplo 3.0: A
program package for imputation and a haplotypebased whole-genome association study using hybrid
parallel computing. Source Code for Biology andMedicine , 6, 10 (2011).
8. A.Fujimoto, H.Nakagawa, N.Hosono, K.Nakano,
T.Abe, K.A.Boroevich, M.Nagasaki, R.Yamaguchi,
T.Shibuya, M.Kubo, S.Miyano, Y.Nakamura,and T.Tsunoda. Whole genome sequencing andcomprehensive variant analysis of a Japaneseindividual using massively parallel sequencing.
Nature Genetics , 42, 931-936 (2010).
9. The International Cancer Genome Consortium.
International network of cancer genome projects.
Nature , 464, 993-998 (2010).
10. A.Fujimoto, Y.Totoki, T.Abe, K.A.Boroevich,
F.Hosoda, H.H.Nguyen, M.Aoki, N.Hoshono,
M.Kubo, F.Miya, Y.Arai, H.Takahashi,
T.Shirakihara, M.Nagasaki, T.Shibuya, K.Nakano,
K.Watanabe-Makino, H.Tanaka, H.Nakamura,
J.Kusuda, H.Ojima, K.Shimada, T.Okusaka,
M.Ueno, Y.Shigekawa, Y.Kawakami, K.Arihiro,
H.Ohdan, K.Gotoh, O.Ishikawa, S.Ariizumi,
M.Yamamoto, T.Yamada, K.Chayama, T.Kosuge,
H.Yamaue, N.Kamatani, S.Miyano, H.Nakagama,
Y.Nakamura, T.Tsunoda, T.Shibata, andH.Nakagawa. Whole Genome Sequencing ofLiver Cancers Identifies Etiological Influences
on Mutation Patterns and Recurrent Mutations in
Chromatin Regulators. Nature Genetics , 44, 760764 (2012).

References
1. K.Ozaki, Y.Ohnishi, A.Iida, A.Sekine, R.Yamada,
T.Tsunoda, H.Sato, H.Sato, M.Hori, Y.Nakamura,and T.Tanaka. Functional SNPs in the lymphotoxinalpha gene that are associated with susceptibility to
myocardial infarction. Nature Genetics , 32, 650654 (2002).
2. T.Tsunoda, G.M.Lathrop, A.Sekine, R.Yamada,
A.Takahashi, Y.Ohnishi, T.Tanaka, andY.Nakamura. Variation of gene-based SNPs andlinkage disequilibrium patterns in the humangenome. Human Molecular Genetics , 13, 16231632 (2004).
3. The International HapMap Consortium. A second
generation human haplotype map of over 3.1
million SNPs. Nature , 449, 851-862 (2007).
4. K.Misawa, S.Fujii, T.Yamazaki, A.Takahashi,
J.Takasaki, M.Yanagisawa, Y.Ohnishi, Y.Nakamura,and N.Kamatani. New correction algorithms for
multiple comparisons in case-control multilocus
association studies based on haplotypes anddiplotype configurations. Journal of Human
Genetics , 53, 789-801 (2008).
5. K.Misawa and N.Kamatani. ParaHaplo: A program
package for haplotype-based whole-genome
association study using parallel computing. Source
Code for Biology and Medicine , 4, 7 (2009).

119

120

Japan (RIKEN CGM)

Extended the power!

The millennium project started
The first GWAS, Nature Genetics.
Biobank Japan Project
The International HapMap Project
The paper on Phase II published
550,000 tagging SNPs
Many GWAS papers published
Commercial platforms
Started genotyping Biobank
GWAS rush
samples with the platforms

550,000 tagging SNPs
Chip/array platforms

2007

2000
2002
2003

The world

Significant milestones in GWAS

4th Biosupercomputing Symposium

Tatsuhiko Tsunoda
Center for Genomic Medicine, RIKENSupercomputing accelerates
genomic medicine
DNA

DNA

AGCATGCTGGTCAGTCGATCGTAGCGAGCATGCTGGTCAGTCGATCGTAGCGAGCATGCTGGTCAGTCGATCGTAGCGAGCATGCTAGTCAGTCGATCGTAGCGAGCATGCTAGTCAGTCGATCGTAGCGAGCATGCTAGTCAGTCGATCGTAGCGNHGRI GWA Catalog
www.genome.gov/GWAStudies

Published Genome-Wide Associations through09/2011
1,596 published GWA at p≤5X10-8 for 249 traits

Test the whole genome:
Genome-wide association study

AG
6173 4209
Cases
Controls 4542 5538

Controls

Cases

Genome-wide Association Study (GWAS)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

#loci Heritability explained
5
50 %
32
20 %
6
15 %
*
180
10 % *
18
6 %
5.2%
7
2.8%
9
4
1.5%

121

Misawa K. et al., J Hum Genet, 53, 789-801 (2008). 7

Calculating the exact probability is very computational intensive.

Let us define the function f by f(…) = 1 if the test is significant at any locus, otherwise f(…) = 0.

Calculation of exact probability under theassumption of a multinomial distribution

Manolio TA et al. Nature, 461, 747-753 (2009)* Lango Allen H et al. Nature, 467, 832-838 (2010)

Age-related macular degeneration
Crohn’s disease
Systemic lupus erythematosus
HeightType 2 diabetes
HDL cholesterol
Early onset myocardial infarction
Fasting glucose

The results from GWAS are insufficient to explain
the heritability of diseases
SNP1 pos.

Misawa K. et al., J Hum Genet, 53, 789-801 (2008). 8This algorithm works much faster than the exact method,
but takes long time

An algorithm for sampling haplotypes under multinomial distribution.
1. As an initial state, arbitrary integer values of xji are given so that their sum is #samples.
2. j=0 or j=1 is selected in equal probability.
3. An integer value u is selected in equal probability from the integers from 1 to L.
4. If xjv > 0, an integer value v other than u (1 < v < L ) is selected.
5. New candidates x*jv = xjv + 1 and x*ju = xju – 1 are calculated.
6. Then, the following value is calculated, c= hvxju / hu(xjv+1) . If c > 1 then {xji } is updated
by substituting x*jv for xjv and x*ju for xju, the step is advanced, a significance test is
performed as described in (8) and the process then returns to (2).
7. If c < 1, then {xji} is updated by substituting x*jv for xjv and x*ju for xju, with probability c,else the state is kept invariant (probability 1 - c), the step is advanced, a significance test is
performed as described in (8) and the process then returns to (2).
8. A test of independence between the phenotype and alleles at each of the l loci is performed
using a chi-square test

Markov Chain Monte Carlo (MCMC) method

For calculating empirical p-values, perform
massively parallel permutation testing on the
K computer: ParaHaplo↓
Compare haplotype frequencies withineach block b/w case and control groups

Standard methods require genome-wide
significance α = 0.05 / #SNPs ~ 10-7
(Bonferroni correction)↓
However, this is too strict because
SNPs are in linkage disequilibrium

Haplotype-based analysis (ParaHaplo)

SNP2 pos.

Missing heritability problem

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Node 2

Node 3…

Node 80,000

122

Diseases/
drug responses

Association between diseases/drug responses
and combination of SNPs

Interaction?

Interaction?

Interlocus interaction

Node 1

11

Genotype Data download from WTCCC↓
Quality Check↓
Detection of LD blocks (Gabriel et al. 2002)↓
Haplotype phasing with Mach 1.0↓
Massively-parallel processing of haplotype-based GWAS

ParaHaplo analysis flow-chart

Rapid association test by importance sampling

Permutation testing: too computationally demanding

Bonferroni correction: too conservative

Genome-wide scan for gene-gene interaction
logit(P) = a1x1 + a2x2 + a3int(x1, x2)

Genome-wide gene-gene interaction

Example output of ParaHaplo

10

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Statistics

Statistics of candidate SNP combination

http://www.genome.gov/images/content/cost_per_genome.jpg

Approximate p-value = red area / total

Average ratio of sampling probability

Sampling probability of permutations from the red area

Frequency

Distribution of statistics based on permutation

Importance Sampling

Effect of
SNP B

Interaction
b/w SNPs123

NGSanalyzer

Reference
Identify polymorphism

Mapping toreference sequence

Sequence

Construct novel sequence

De novo assemble

AdaptorSequenced region

DNA fragmentation and library preparation

Sample DNA

Whole genome sequencing

Effect of
SNP A

logit
=
( P ) µ + α x A + β xB + γ x A xB

P: probability of being affected

Logistic regression model with a SNP-SNP interaction term

For linear regression models
and likelihood ratio test

Extended Rapid Association Test
(ExRAT)

4th Biosupercomputing Symposium (December 3 - 5, 2012)IndelHeterozygous
deletion
Gain

Rearrangement
Non-reference
human genome

Meyerson, M. et al., Nature Reviews Genetics, 11, 685- 696 (2010).

CNV(Copy number variation)

Homozygous
deletion

124

Nature(2010).
NatureGenetics,
Genetics,42,
42,931-936
931-936(2010).

First WGSWGS andand analysis
analysis ofof aa JapaneseJapaneseFirst
genomegenomeSNV(single nucleotide variant)

Genetic variation detection

NatureGenetics,
Genetics,44,
44,760-764
760-764(2012).
Nature(2012).

WGS and analysis of liver cancers

4th Biosupercomputing Symposium (December 3 - 5, 2012)

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Towards efficient improvement oftranscriptional circuit models by Life
Science Data Assimilation System (LiSDAS)
Tomoyuki Higuchi
Director-General, Institute of Statistical Mathematics

Profile:Tomoyuki Higuchi is currently Director-General of the Institute of Statistical Mathematics (ISM), Japan. He
is also Professor of ISM and of the Graduate University for Advanced Studies. He obtained the B.S. degree in
1984, the M.S. degree in 1986, and the Ph.D. degree in 1989 from the University of Tokyo. His primary researchinterests are in a Bayesian modeling and sequential Monte Carlo computation, in particular, data assimilation.

125

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Abstract
1. Introduction
It has been about 40 years since a computer simulation
is called as the third method of science. A recent
explosion of data, the so-called appearance of big data,
strengthens the research domain to study a method oftools for analyzing big data such as statistics, machinelearning, data mining, and visualization technologies.
This phenomenon is called the fourth paradigm
after a publication. The data assimilation (DA) is asynthesis technique based on the Bayesian filtering
method by embedding observation/experiment data
in a numerical simulation. DA is an emerging area inearth sciences, particularly oceanography. Its research
motivation is easily understood simply if we noticethat there are too many uncertainties in the modelsuch as the boundary condition, initial condition,
unknown parameters, and unknown dynamics. It
yields an accommodation ability to make a simulation
real, and the better initial and boundary conditions canbe automatically obtained.
Major objectives of DA are classified into thefollowing five aspects. The first is to find the best or
better initial condition for forecasting. It is actuallyrealized in the real weather forecast. The second
is to find the best or better boundary condition inconstructing a simulation model. This procedureincludes a setting of appropriate boundary conditions
necessary for dealing with the coupled phenomena.
The third is to optimize the configuration ofmacroscopic parameters which substitute a descriptionof complicated real phenomena involving different
temporal and spatial scales. A validation of theempirically given values is included in this problem.
The fourth is to inter/extrapolate physical quantitywith some numerical simulation model at times andlocations where no observation data is available.
This procedure is called “a generation of re-analysis
dataset (product)”. The fifth is to perform a sensitivity
analysis of virtual observation network in an attempt
to construct an effective observation network systemwith less budgetary cost and less consuming time.

2. Applications
We are studying the sequential DA methods basedon the Ensemble Kalman and particle filters, andconducting the DA experiments in several specificareas. We will give a brief explanation for thesequential DA and demonstrate a part of applications
carried out by our DA research group.
The understanding of transcription circuit is notenough from both the experimental and the theoreticalviews. Life Science Data Assimilation System(LiSDAS) has been developed as an application
dedicated for K computer, with the aim of proposing
improved models by assimilating observation data to aconventional model. LiSDAS describes transcriptional
circuits in terms of differential equations, which
include processes of protein-protein interactions,the activation or the repression of mRNAs bytranscription factors, and the degradation of mRNAs
and proteins. We will make a demonstration in ourtalk, employing two network models: a circadianrhythm network model and a statistically constructed
model for network of mRNAs in lung cancer cells.
Data assimilation works also in intracellular fluid
dynamics. A large flow is observed in a fertile egg
right before a cell division its mechanism is yet
uncertain. One hypothesis is that proteins “myosin”
in the cell wall drive physically such intracellularflow. We have performed data assimilation in order toobtain dynamical evidence supporting this hypothesis.
Our model obtained by data assimilation successfully
explains the observed cytoplasmic flow.
We have also applied DA to influenza epidemics
and estimated connection of cities in the context of2009 influenza pandemic. It is important to knowhow degree cities are coupled since such preventive
activities are conducted basically by individualautonomies.
References
1. T. Higuchi, Proceedings of 14th International
Conference Fusion (2011).

126

4th Biosupercomputing Symposium

Ryo Yoshida Hiromichi Nagao

Shin’ya Nakano

127

Collaboration with Miyano, S and Gotoh, N (HGC, Univ.Tokyo)

Gene regulation models for lung cancer cells withacquired resistance to anticancer drugs
Niwayama, R and Kimura, A (NIG)

3Cell mechanics for cytoplasmic flow in the
C. elegans embryo

Applications: Two types of inversion problems using biochemical and biophysics simulators

Developments: Data assimilation software implemented on K-computer

Members: Saito, M, Yoshida, R, Nagao, H, Nakano, S, Higuchi, T (ISM)

LiSDAS - Life Science Data Assimilation Systems

Masaya M. Saito

The Institute Statistical Mathematics
Tomoyuki Higuchi

Data assimilation in life-science=yn g (x n ) + noise– Relate the model output to given data=x n fG (x n −1 , q ) + noiseState variables

Data

D = {yn }ParameterModel

x = {x n }fG ( ⋅ q )GState Space Representation

Data Assimilation – Method of Integration for Simulation and Data

4th Biosupercomputing Symposium (December 3 - 5, 2012)

ModelParameter

State variables

DataG
fG ( ⋅ q )x = {x n }D = {yn }ω can be any of G, θ, x128

Variable: RNA, protein, metabolite, carbohydrate
Process: Production, binding, degradation, diffusion

7

KEGG: www.genome.jp/kegg/

Functions are programmed into Networks

••Cell - an integrated device made of several types of biochemical reactions involving
diverse molecules

Biological Networks

Provided by Dr. Kimura, A (National Institute of Genetics)p ( w D ) ∝ p ( D w) p ( w)

Bayesian Inference – Solving the Inversion Problem• 29 simulation variables• 7 observation variables• 116 Model parameters

Scale of the model

8

Example1: Gene regulatory systems for circadian clock controls

w ~ p ( w D ) ∝ p ( D w) p ( w)

Bayesian Supercomputing

4th Biosupercomputing Symposium (December 3 - 5, 2012)

129

Normal strain

Drug resistant Strain― final MAP solution― TF threshold

Periodic behaviour in data is reproduced

ISLiM成果報告会2011

All genes (17654)■: ／increasing■:- Stationary■:＼decreasing

9

Time (up to 48hrs)

TFs, their concentration being below/above
thresholds for all time, are interpreted as
removal from the regulatory circuit.
(Network learning)― intermediate MAP

Example 2: Whole genes’ transcription simulators○ dataCircadian clock gene regulatory simulation-2

0 0.5 1 2 3 4 5 6

EGF (100 ng/ml)( - ) : 0 hr ~ 24 hr
EGF : 0 hr ~ 24 hr
EGF + Gefitinib : 0 hr ~ 24 hr
Gefitinib : 0 hr ~ 24 hr

RNA Sampling(hours)7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

cell
#3

?
Hypothetical network

cell
#2

cell
#1

(

(

)

)

I ijgap (t ) = gijgap (V j (t ) − Vi (t ) )1 + exp − β V j (t ) − V jeq( (gijchem Eijchem − Vi (t )

))

(

)

dVi (t )1= − m Vi (t ) − EiLeak + ∑ nijchem I ijchem (t ) + nijgap I ijgap (t ) + I istim (t )dtRijI ijchem (t ) =Cim

Collaboration with Iino (Univ Tokyo), Ishibara(Kyusu Univ), Iwasaki (Ibaraki Univ)

x1
x2
x3
x4

Spatio-temporal image of nerve cells of C.elegans

4th Biosupercomputing Symposium

12

Nerve cell activities are measured along time via the 4D calcium imaging data.
LiSDAS would promote the deep understanding of the complex network system behind.

Simulation

cell
#4

Gefitinib

AND/OR

Post GC Application 1: Neuronal circuit model

gefitinib (0.5 µM)

Data

- 48

Starvation

EGF

Example 2: Drug response prediction for non-small cell lung cancer

4th Biosupercomputing Symposium (December 3 - 5, 2012)

(P)Data: Cytoplasmic flow from the anterior to the posterior

Schematic illustration of intercelluar fluid of C.elegans embryos
Collaboration with Dr Kimura (National Institute of Genetics)

QuestionIs the observed flow (velocity field) reproduced under
the hypothesis?

Hypothesis: Myosin on the membrane is a driving forceof the flow.

130

Initial compound - toluene

Pattern recognition andfeature extraction

Compound assay data

15

Targets - COX2 activating compounds

Alter the chemical structure by ‘mixing’ the identified key
molecules using the fragment-flip MCMC method• Collaboration of biophysics and statistical science• Random sampling of chemical compound graph using a kernel
machine and Bayesian computational techniques

Developments: Computational method for generating de novo chemical compounds

Members: Yamashita, H (GUAS), Iba, Y (ISM), Higuchi, T (ISM), Yoshida, R (ISM)

Post GC Application 4:
Data-driven de novo drug design using kernel machine learning

Myosin(A→P)-flow, unknown mechanism(A)

Post GC Application 2: Intercellular fluid dynamics

•

•

2010

Hokkaido

Mie

Miyagi

multimodality

2011

2012

Synchronous Abrupt changes

εij

4th Biosupercomputing Symposium• Contact us
http://daweb.ism.ac.jp/• Potential applications– Neural cell network– Intercellular fluid– Infectious transmission network– (Any whenever you have data and model)• Data assimilation system «LiSDAS» has been developed for
network system analyses for cells

Summary

Goal: Learn connection strength εij between prefectures fromfeatures seen in time courses.
Such a network model help us to incorporate transmission
dynamics into the forecast of influenza epidemic.

2009

Time shift in peak

The number of Influenza patients in recent three seasons

Post GC Application 3: Infectious transmission network

16

14

4th Biosupercomputing Symposium (December 3 - 5, 2012)

Contact:
RIKEN
Computational Science Research Program
2-1 Hirosawa, Wako, Saitama 351-0198, Japan
http://www.csrp.riken.jp