High-Performance Systems Biology and

Associated Combinatorial Scientific

Computing Problems

David Alber

Materials and Computational Science Center

National Renewable Energy Laboratory

June 11, 2008

NREL/PR-590-43765

Presented at the 2008 Workshop on Combinatorial Scientific Computing & Petascale Simulations held June 10-13, 2008 in Sant Fe, New Mexico.
Project Overview

Project Participants and Funding

Collaborators:

NREL SCG: Chris Chang, Peter Graf, and Kwiseon Kim*

NREL Photobiology: Mike Seibert*

Summer Student from CU Boulder: David Biagioni

Participating institutions:

National Renewable Energy Laboratory

Colorado School of Mines

Stanford University

Funding through SciDAC (OASCR and OBER)

David Alber CSCAPES Workshop 2008 2/18
Project Overview

Metabolism and Metabolic Modeling

Metabolism

Chemical reactions occurring in living cells

Reactions catalyzed by enzymes

Metabolic species (e.g., glucose, pyruvate) produced and

consumed

Reaction Modeling

Several models

Michaelis-Menten kinetics

David Alber CSCAPES Workshop 2008 3/18
Project Overview

High-Performance Systems Biology

In a nutshell

Model complete metabolism of

Chlamydomonas reinhardtii

Develop high-performance software to

explore metabolism kinetics

Example problems

Parameter estimation (data fitting)

Sensitivity minimization

Parameter space characterization

David Alber CSCAPES Workshop 2008 4/18
Project Overview

Metabolic Model

For metabolic reaction:

dy

dt = f(y, k, E)

y – vector of metabolite concentrations

k – vector of kinetic parameters

E – vector of enzyme concentrations

Kinetic parameters:

1 Time consuming to determine experimentally

2 Essential to understanding metabolic kinetics

David Alber CSCAPES Workshop 2008 5/18
Project Overview

Metabolic Model

For metabolic reaction:

dy

dt = f(y, k, E)

y – vector of metabolite concentrations

k – vector of kinetic parameters (mostly unknown)

E – vector of enzyme concentrations

Kinetic parameters:

1 Time consuming to determine experimentally

2 Essential to understanding metabolic kinetics

David Alber CSCAPES Workshop 2008 5/18
Project Overview

Metabolic Model

For metabolic reaction:

dy

dt = f(y, k, E)

y – vector of metabolite concentrations

k – vector of kinetic parameters (mostly unknown)

E – vector of enzyme concentrations

Kinetic parameters:

Time consuming to determine experimentally

Essential to understanding metabolic kinetics

David Alber CSCAPES Workshop 2008 5/18
D-glucose

D-glucose

Glucose

6-phosphate

Fructose

6-phosphate

Fructose 1,6-

bisphosphate

Dihydroxyacetone

phosphate

(D)-Glyceraldehyde

3-phosphate

1,3-Bisphosphoglycerate

3-Phospho-(D)-

glycerate

2-Phospho-(D)-

glycerate

Phosphoenolpyruvate

Pyruvate

phosphoglucose

isomerase

phosphofructokinase

aldolase

triosephosphate

isomerase glyceraldehyde 3-phosphate

dehydrogenase

+ NAD

NADH

Pi

phosphoglycerate mutase

enolase

pyruvate

kinase

ADP

ATP

ADP

ATP

ADP

ATP

glucokinase

phosphoglycerate

kinase

external

cytosol

Pyruvate

acetyl-CoA

citrate

isocitrate

a-ketoglutarate

succinyl-CoA succinate

fumarate

(S)-malate

oxaloacetate

aconitase

isocitrate

dehydrogenase

a-ketoglutarate

dehydrogenase

succinyl-CoA ligase

succinateubiquinone

oxidoreductase

fumarase

malate

dehydrogenase

citrate

synthase

CoA

+ NAD

NADH

CO2

+ NAD

NADH

CO2

CoA Pi

ADP ATP

O

|| U

OH

|U

+ NAD

NADH

pyruvate

dehydrogenase

CoA

CO2

+ NAD

NADH

cytosol

mitochondrial matrix

pyruvate

acetyl-CoA

acetyl

phosphate

acetate

phosphotransacetylase

acetate

kinase I

ADH1

CoA

CO2

PFORox

PFORred

CoA

ADP

ATP

Pi

CoA

hydrogenase

H2 2 Ferrox

2 Ferrred

NADH

+ NAD

ethanol

acetaldehyde

pyruvate

decarboxylase CO2

NADH

+ NAD

- HCO2

pyruvate

formate-lyase

CoA
Project Overview

Parameter Estimation

Find set of k such that species concentrations match target

values

Expressed as optimization problem:

min

k

g(k),

where g(k) = krk

2

2 and r is vector of differences between

target and simulated values

David Alber CSCAPES Workshop 2008 7/18
Project Overview

Sensitivity Minimization

Find set of k to minimize sensitivity

Applies to organism engineering issues and general

“nature of life” questions

Objective function to minimize:

h(k) = kJy(k)k

2

F

=

Xn

i=1

Xm

j=1



∂yi

∂kj

2

David Alber CSCAPES Workshop 2008 8/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

yi: simulated metabolite

concentration

y ̄i: target concentration

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

Sensitivity minimization:

h(k) = Xn

i=1

Xm

j=1



∂yi

∂kj

2

Single entry of ∇h:

∂h

∂k`

= 2X

i,j

∂yi

∂kj

∂

2yi

∂kj∂k`

Computing second derivative

term expensive

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

Sensitivity minimization:

h(k) = Xn

i=1

Xm

j=1



∂yi

∂kj

2

Single entry of ∇h:

∂h

∂k`

= 2X

i,j

∂yi

∂kj

∂

2yi

∂kj∂k`

Computing second derivative

term expensive

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Gradient-based Optimization

Parameter estimation:

g(k) = Xn

i=1

(yi − y ̄i)

2

Single entry of ∇g:

∂g

∂kj

= 2Xn

i=1

(yi − y ̄i)

∂yi

∂kj

No problems computing ∇g

(adjoint sensitivity analysis)

Sensitivity minimization:

h(k) = Xn

i=1

Xm

j=1



∂yi

∂kj

2

Single entry of ∇h:

∂h

∂k`

= 2X

i,j

∂yi

∂kj

∂

2yi

∂kj∂k`

Computing second derivative

term expensive

David Alber CSCAPES Workshop 2008 9/18
Project Overview

Second Derivative Computation

CSC Problem Slide #1

Finite differences

Automatic differentiation

Collaboration with Paul Hovland

Collaboration with Radu Serban (formerly at LLNL):

Combines adjoint sensitivity analysis and AD (via

Tapenade)

Parallelized

Received code two weeks ago (i.e., nothing yet to show)

David Alber CSCAPES Workshop 2008 10/18
Project Overview

Sensitivity Analysis Cost Substantial

≈1000 dimension parameter

space (ultimately)

Each evaluation of h(k)

inexpensive (similar to cost of

∇g)

Evaluation cost of ∇h adds up

David Alber CSCAPES Workshop 2008 11/18
Project Overview

Load Balancing

CSC Problem Slide #2

The hierarchical parallelism in this application generates

load balance issues

More of the scheduling problem variety

Want to hear more on this

David Alber CSCAPES Workshop 2008 12/18
Project Overview

Kinetics Complications

Many gaps to fill regarding metabolic kinetic rates

Requires much experimentally-obtained data

Relevant experimental data difficult (tedious?) and

expensive to generate

Switch gears and discuss another systems biology

(metabolomics) problem

David Alber CSCAPES Workshop 2008 13/18
Project Overview

Stoichiometric Matrix

Stoichiometry refers to numbers in chemical formulae

2A + B GGGA C + 3D

C + 2E GGGA F

Stoichiometric matrix stores stoichiometry for all reactions

S =





−2 0

−1 0

1 −1

3 0

−2 0

0 1





A

B

C

D

E

F

1

2

David Alber CSCAPES Workshop 2008 14/18
Project Overview

Stoichiometric Matrix and Flux Analysis

Previously said dy/dt = f(y, k, E)

Can now also say

dy

dt = Sv,

where v is vector of reaction

velocities

No kinetic parameters to get in

the way!

At steady state, dy/dt = 0

(almost)

Now looking at Sv = 0

David Alber CSCAPES Workshop 2008 15/18
Project Overview

Problem Definition

Say we want to maximize concentration of species yi:

max

v

yi

This problem trivial without constraints

Possible constraints:

Thermodynamics (an arbitrarily large reaction rate not

possible)

Sv = 0 at steady state (almost)

All vi ≥ 0 (reactions do not happen in reverse – not exactly

true)

David Alber CSCAPES Workshop 2008 16/18
Project Overview

Add Another Layer

Suppose goal is to dispose of some reactions or pathways

Which pathways are more critical?

Conceptual connections to electrical grid modeling?

David Alber CSCAPES Workshop 2008 17/18
Project Overview

Molecular Dynamics

Shameless plug for another SciDAC project

Application:

Study function of CBH I enzyme

in “digesting” crystalline cellulose

Processivity?

Challenges:

CHARMM does not scale well,

has everything needed for

application

Other packages scale well, do

not have everything needed

If CHARMM: load balancing

David Alber CSCAPES Workshop 2008 18/18