Exploring Complex Reaction Pathways Mahmoud Moradi Department of - - PowerPoint PPT Presentation

Mahmoud Moradi

Department of Chemistry and Biochemistry University of Arkansas

"Hands-on" Workshop on Enhanced Sampling and Free-Energy Calculation at Urbana, IL September 14, 2018

Exploring Complex Reaction Pathways

Outline

• Introduction

– How to study large-scale conformational changes?

• Methodology

– Empirical search for good pulling protocols – Iterative combination of free energy calculation methods and path-finding algorithms

Outline

• Introduction

– How to study large-scale conformational changes?

IF apo IF bound in

• ut

in

• ut

OF apo

• ut

in

• ut

in OF bound

Large-Scale Conformational Changes in Membrane Transport Proteins

• Membrane transporters rely on

large-scale conformational changes between inward-facing (IF) and

• utward-facing (OF) states

(alternating access mechanism)

• Channels may require large-scale

conformational changes between their open/active and closed/inactive states.

• pen/active

closed/inactive

A Case Study: Proton-coupled Oligopeptide Transporters (POTs)

GkPOT (4IKV, 1.9 Å) ~100,000 atoms Conventional unbiased simulations performed: 8 conditions × 400 ns × 2 repeats = 6.4 𝜈𝑡

K Immadisetty, J Hettige, and M Moradi, Wha hat C Can a n and C nd Cannot nnot B Be L Learne ned f d from

• m

Mol

• lecul

ular D Dyna ynamics S Simul ulations

• ns of
• f B

Bacterial P Prot

• ton
• n-Coupl
• upled

d Oligope gopept ptide de Trans nspor porter GkP kPOT? J. Phys. Chem. B, 121:3644-3656, 2017.

N

• B

u n d l e C-Bundle

N-,C-Bundle Interdomain Angle H5,H11 Interhelical Angle H1,H7 Interhelical Angle L4,5-L10,11 Distance

R43 E310

Monitoring Global and Local Conformational Changes

Global

Local Conformational Changes

2 3 4 5 6 7 8 9

UP: E310 unprotonated

UP:apo(Set−1) UP:apo(Set−2) 3 4 5 6 7 8 9

B

R43−E310 Distance(Å)

UP:AA(Set−1) UP:AA(Set−2) 100 200 300 400

Time (ns)

P:apo(Set−1) P:apo(Set−2)

F

P:AA(Set−1) P:AA(Set−2) 100 200 300 400

Time (ns) P: E310 protonated

There is a clear distinction between different conditions.

10 15 20 25 30 35 40 45 50

C−,N−Bundle Interdomain Angle(°)

UP:apo (Set−1) P:AA (Set−1) UP:apo (Set−2) P:AA (Set−2) 10 13 16 19 22 25 28 100 200 300 400

L4,10−L5,11 Distance(Å) Time (ns)

UP:apo (Set−1) P:AA (Set−1) 100 200 300 400

Time (ns)

UP:apo (Set−2) P:AA (Set−2)

Global Conformational Changes

There is no clear distinction between different conditions.

Global Conformational Changes

“Structural basis for dynamic mechanism of proton-coupled symport by the peptide transporter POT.” PNAS 2013 | vol. 110 | no. 28 | 11343–11348.

Although a common practice, statements made about millisecond-level biomolecular events based on sub- microsecond level simulations are not reliable.

Global Conformational Changes

20 25 30 35 40 45

C−,N−Bundle Interdomain Angle (°)

Set−1 Set−2 45 50 55 60 65 70

H1,H7 Interhelical Angle (°)

Set−1 Set−2 15 20 25 30 35 40

H5,H11 Interhelical Angle (°)

Set−1 Set−2 12 15 18 21 24 27

L4,5−L10,11 Interloop Distance (Å)

Set−1 Set−2 UP:apo UP:AA UP:EE UP:FF P:apo P:AA P:EE P:FF

Simulation System

UP:apo UP:AA UP:EE UP:FF P:apo P:AA P:EE P:FF

Simulation System

There is no statistically significant distinction between different conditions.

• Introduction

– How to study large-scale conformational changes?

It is not reasonable to speculate about the conformational transition between two states based

• n fluctuations around one of the end points.

A B

How to study large-scale conformational changes?

Inward-Facing Outward-Facing

GkPOT

GkPOT How to study large-scale conformational changes?

If you are not sure how good your collective variable is, start with nonequilibrium pulling.

How to study large-scale conformational changes?

• Introduction

– How to study large-scale conformational changes?

• Methodology

– Empirical search for good pulling protocols – Iterative combination of free energy calculation methods and path-finding algorithms

Sampling Ideas

– Free energy calculations require dimensionality reduction. – Traditionally, this is done by designing intuitive, ad- hoc, knowledge-based collective variables. – Another approach is to use data-driven collective variables using standard dimensionality reduction techniques (PCA, diffusion maps, etc). – Alternatively (or in combination with the above approaches), one can calculate free energy along a transition path (a 1D curve). – The path can be obtained from path-finding algorithms. – Since sampling is never perfect the procedure could be iterative to reach higher accuracies.

Sampling Ideas

• Reaction coordinates

– System-specific collective variables

• Searching for efficient pulling protocols

– An empirical approach to sampling

• Along-the-curve free energy calculations

– Free energy calculations combined with path-finding algorithms

• Iterative sampling

– A posteriori tests of self-consistency

Moradi et al., Proc Natl Acad Sci 106 20746 (2009) Moradi et al., Chem Phys Lett 518 109 (2011) Moradi et al., J Chem Phys 133 125104 (2010) Moradi et al., Int J Quantum Chem 110 2865 (2010) Moradi et al., Biophys J 100 1083 (2011) Moradi et al., J Phys Chem B 115 8645 (2011) Moradi et al., PLoS Comput Biol 8 e1002501 (2012) Moradi et al., Nucleic Acid Res 41 33 (2013) Moradi et al., Proc Natl Acad Sci 110 18916 (2013) Moradi et al., J Phys Chem Lett 4 1882 (2013) Moradi et al., Methods Mol Biol 924 313 (2013) Moradi et al., J Chem Phys 140 034114 (2014) Moradi et al., J Chem Phys 140 034115 (2014) Moradi et al., J Chem Theory Comput 10 2866 (2014) Moradi et al., J Phys Conf Ser 640 012014 (2015) Moradi et al., J Phys Conf Ser 640 012020 (2015)

Moradi et al., Nat Commun 6 8393 (2015) Fakharzadeh & Moradi, J Phys Chem Lett 7 4980 (2016)

Empirical search for practical collective variables for inducing the conformational changes involved in the transition. Systematic search for a practical biasing protocol by using different combinations of collective variables.

I.2 Optimizing the Biasing Protocols I.1 Defining Practical Collective Variables

Use all of the conformations available to generate the most reliable transition pathway:

• 1. Bayesian approach for combining the data
• 2. Post-hoc string method (analysis tool)
• 3. String method with swarms of trajectories
• II. Optimizing the

Transition Pathway III.1 Free Energy Calculations

Using the most relevant collective variables (from I.1), biasing protocol (from I.2), and initial conformations (from I.2).

III.2 Assessing the Sampling Efficiency

Detecting the poorly sampled, but potentially important regions, e.g., by using PCA.

Sampling Ideas

• Reaction coordinates

– System-specific collective variables

• Searching for efficient pulling protocols

– An empirical approach to sampling

• Along-the-curve free energy calculations

– Free energy calculations combined with path-finding algorithms

• Iterative sampling

– A posteriori tests of self-consistency

Work Reaction Coordinate

Empirical search for reaction coordinates and biasing protocols

Optimized Protocol F r e e E n e r g y C a l c u l a t i

• n

s P a t h

• R

e f i n i n g A l g

• r

i t h m s

Sampling Ideas

• Introduction

– How to study large-scale conformational changes?

• Methodology

– Empirical search for good pulling protocols – Iterative combination of free energy calculation methods and path-finding algorithms

Example: Glycerol-3-Phophate Transporter (GlpT)

• Major facilitator

superfamily (MFS)

• Secondary active

transporter

• Crystalized only in the IF

state.

• GlpT transports G3P using

Pi gradient.

• Pi:Pi exchanger (in the

absence of organic phosphate)

periplasm cytoplasm

Pi G3P

Moradi, Enkavi, and Tajkhorshid, Nature Communications 6 8393 (2015)

Pi

H7

Periplasm Cytoplasm

H1

IF

b

OF

b

Cytoplasm Periplasm

H5 H11

OF

a

IF

a

Transport Thermodynamics

Lemieux, et al., Curr. Opin. Struct. Biol. 14, 405 (2004).

OFb IFb IFa OFa

Free Energy barrier Free Energy Transition State Reaction Coordinate

Law, et al., Biochemistry 46, 12190 (2007).

Transport Thermodynamics

a: apo b: bound

Pi

H7

Periplasm Cytoplasm

H1

IF

b

OF

b

Cytoplasm Periplasm

H5 H11

OF

a

IF

a

Full Thermodynamic Cycle

Pi

H7

Periplasm Cytoplasm

H1

IF

b

OF

b

Cytoplasm Periplasm

H5 H11

OF

a

IF

a

the only available crystal structure

Pi

H7

Periplasm Cytoplasm

H1

IF

b

OF

b

Cytoplasm Periplasm

H5 H11

OF

a

IF

a

Step 1: OFaç

èIFa

about 100 simulations with different protocols

IFèOF Nonequilibrium Work

Empirical search for reaction coordinates and biasing protocols:

IFç

èOF transition induced by imposing rotational change on

helices TM1 and TM7

Number of water molecules per Å (averaged over a 1 ns window)

IF OF

10 ns IF equilibrium 20 ns nonequilibrium (IFèOF) 10 ns OF equilibrium

θ7 θ1

Example: MsbA Transporter

• ATP-binding cassette (ABC) exporter with three crystal

structures.

• Reaction coordinates:

α, β, γ (relative orientation of different domains)

Nucleotide-binding domains (NBD) Two IF conformations: IF-closed (IF-c) and IF-open (IF-o) different views

Ward A., Reyes C. L., Yu J., Roth C. B., Chang G.. PNAS 104 19005 (2007)

Moradi and Tajkhorshid, PNAS 110 18916 (2013)

Conventional Equilibrium Molecular Dynamics

5 10 15 20 25 30 35 10 20 30 40 50 60

b a

OF IF-c IF-o 5 10 15 20 25 30 35 20 30 40 50 60 70 80 90

b g

OF IF-c IF-o apo MsbA in explicit water/membrane (300 ns)

400

α

200

γ β

100 200 300 400 500

40 80 120 160 work (kcal/mol) t (ns)

(a)

→α→γ→β →γ→α→β → →α→γ→β →α→β→γ → →γ→α→β → →α→β→γ →β→ →α→γ →β→α→γ → →β→α→γ →β→γ→α →γ→β→α → →γ→β→α →β→ →γ→α → →β→γ→α

α→β β→α→γ (β,γ)→α

160

Nonequilibrium work along select trajectories Optimized protocol Moradi and Tajkhorshid, PNAS 110 18916 (2013)

400

α

200

γ β

100 200 300 400 500

40 80 120 160 work (kcal/mol) t (ns)

(a)

→α→γ→β →γ→α→β →d→α→γ→β →α→β→γ →d→γ→α→β →d→α→β→γ →β→d→α→γ →β→α→γ →d→β→α→γ →β→γ→α →γ→β→α →d→γ→β→α →β→d→γ→α →d→β→γ→α

} } }

α→β β→α→γ (β,γ)→α

160

Moradi and Tajkhorshid, PNAS 110 18916 (2013)

Steering along orientation quaternion

A quaternion is a 4-component mathematical object:

𝑟 = (𝑟,, 𝑟., 𝑟/, 𝑟0)

To find the best rigid body rotation for {𝒚𝑙} → {𝒛𝑙}, find a unit quaternion 𝑟 to minimize:

𝑟 8(0, 𝒚:)𝑟 8∗ − (0, 𝒛:) /

It turns out:

𝑟 = (cos 𝜄 2 , sin 𝜄 2 𝒗 E)

Biasing Potential: 𝜄 and 𝒗 E are the angle and axis of rotation target quaternion at time t

Interpolation of orientation quaternion in colvars At each time t+Δt: (1) Linear interpolation: (2) Normalization:

Interpolation of orientation quaternion in NAMD It turns out: that can be used for work measurements using:

• Introduction

– How to study large-scale conformational changes?

• Methodology

– Empirical search for good pulling protocols – Iterative combination of free energy calculation methods and path-finding algorithms

!!

!! = !!!!

Non-Parametric Reweighting

!! ! ∝ !(!)!!!(!)! !

! ! = ! ! + !! ! !

Unbiased potential Biased density Biasing potential Biased potential Unbiased density Biasing factor Known Unknown Known Partition function Unknown

!! ! = !!!!(!)!

!! = !"!!(!)!!(!)!

Conventional WHAM

• S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, P. A. Kollman, “The weighted

histogram analysis method for free-energy calculations on biomolecules.”

• J. Comput. Chem. 13, 1011 (1992)

!! = !! !!!!

!!!! ! !

!

!! ≡ !!

! !

, !! ≡ !!

! !

!

!! = !!!!

! !

! !!

! ≡ !! !! !

If#the#!#space#is#discrete#

• r#“binable”#such#that#

!!(!) ≈ !!(!!)#

e.g.,#Ul#is#a#smooth#function#of# a#1D#coordinate#

!! ! ≈ !! ! ,# ! ≡ [! − !! ∆! ]#

Generalizations

• C. Bartels, “Analyzing biased Monte Carlo and molecular dynamics simulations.”
• Chem. Phys. Letters 331, 446 (2000)

!! = !!!!

! !

! !! = 1 !!!!

!!!! ! !

!

Every conformation sampled is a state. Solved iteratively.

Generalizations

!! = !!!!

! !

! !! = 1 !!!!

!!!! ! !

!

!

! = −!"#

!!

!

!!!!

!!!

! ! !

! ! !! = !!!!!

• M. R. Shirts, J. D. Chodera, “Statistically optimal analysis of samples from multiple

equilibrium states.”

• J. Chem. Phys., 129, 124105 (2008)

!! = !!

!

!!!!

!!!! ! ! !

!

Multi-state BAR (MBAR) equation

Combining path-finding and free energy methods

• Potential of Mean Force:

𝐻 𝜼 = −𝛾I. log 𝜀(𝝄(𝒚) − 𝜼) 𝜀(𝝄(𝒚) − 𝜼) = N 𝜀(𝝄(𝒚) − 𝜼)𝜍(𝒚, 𝒒) 𝑒0R𝑦𝑒0R𝑞

• Perturbed Free Energy:

𝐺

V = 𝐺 𝜼V = −𝛾I. log 𝑎V

𝑎V = N 𝑓IYZ[ 𝝄 𝜍(𝒚, 𝒒) 𝑒0R𝑦𝑒0R𝑞 = N 𝑓IY(\ 𝝄 ]Z[ 𝝄 )

• 𝑒_𝜊
• Non-parametric MLE estimates:

Shirts, Chodera, JCP, 129, 124105 (2008)

e−βF

i =

e−βUi (ξ t ) nje

−β(U j (ξ t )−Fj ) j

∑

all samples

eβF

i =

ni wte−βUi (ξ t )

t

∑

wt = 1 e−β(Ui (ξ t )−F

i )

i

∑

Bartels, CPL, 331, 446 (2000)

Pi

H7

Periplasm Cytoplasm

H1

IF

b

OF

b

Cytoplasm Periplasm

H5 H11

OF

a

IF

a

Moradi, Enkavi, and Tajkhorshid, Nature Communications 6 8393 (2015)

2 4 6 8 10 12 14 16 18 20 40 60 80 100 150 140 130 120 110 100

Free Energy (kcal/mol) Image Index Image Index

IFa IFb OFb OFa TSb TSa

→unbinding ←binding →binding ←unbinding

Distinct conformational transition pathways

Quaternion-based principal components (QPCs) represent different modes of concerted motions of transmembrane helices.

with substrate without substrate

Free Energy Calculation Path-Finding Algorithm

Iterative path-finding algorithms and free energy calculations

Bias-Exchange Umbrella Sampling (Free Energy Calculation) String Method with Swarms of Trajectories (Path-Refining Algorithm)

Iterative path-refining algorithms and free energy calculations

Bias-Exchange Umbrella Sampling (Free Energy Calculation) String Method with Swarms of Trajectories (Path-Refining Algorithm) Post-Hoc String Method (Analysis)

Iterative path-refining algorithms and free energy calculations

Moradi, Enkavi, and Tajkhorshid, Nature Communications 6 8393 (2015)

Post-hoc string method (PHSM)

• Suppose we have already sampled a particular “continuous”

region of configuration space (or some multi-dimensional collective variable space ) and estimated the weight of each sample .

• PHSM finds the principal curve in the space using available

samples as an approximate minimum free energy pathway. A B

{ξi} {wt} ξ

Moradi, Enkavi, and Tajkhorshid, Nature Communications 6 8393 (2015)