2010.11.4-6 第13回情報論的学習理論ワークショップ（IBIS2010）

１分子計測データに対するパラメータ推定
Parameter estimation on single-molecule time series

京都大学理学研究科 物理学専攻
宮崎牧人

Dept of Physics, Kyoto Univ
Makito Miyazaki

1Introduction
progress in statistical physics of non-equilibrium small systemsex.) Jarzynski equality, fluctuation theorem, Hatano-Sasa equality, Harada-Sasa equality, ...

angle

but ... the model parameters should be determined.

insight

invisible

visible

time

probe
F1-ATPase

substratum

model

development in a single-molecule observation & manipulation techniqueex.) optical microscope, FRET, TIRF, optical tweezer, AFM, ....

but ... observable degrees of freedom are limited.

problem: there are no general frameworks to derive system informationfrom restricted measurements.
2

General framework
invisible degrees of freedom： x = (x1 , x2 , . . . , xm )visible degrees of freedom： y = (y1 , y2 , . . . , yn )system parameter：

trajectory of x and y are denoted as [x] and [y]∗ : true value,
problem :simple idea :ˆ: estimated value

How can we estimate the values of Πon the bases of [y] WITHOUT referring to [x] ?
maximize P (Π|[y]) with respect to Π

Are there any general evidence that the idea actually works ?

3

General framework
simple idea :maximize P (Π|[y]) with respect to Π

P ([y]|Π)P (Π)(Bayes theorem)
P (Π|[y]) =P ([y])
P (Π) = const. ( noninformative prior )marginal likelihood:�P ([y]|Π) ≡Dx P ([x, y]|Π)sum up all possible

path probability : P ([x, y]|Π)displacement

displacement

P (Π|[y]) ∝ P ([y]|Π) (marginal likelihood)yx

trajectory of x

timetime
4General framework
The maximization of P (Π|[y]) is replaced with the minimization of“Hamiltonian” deﬁned as Hτ (Π)[y] ≡ −τ −1 ln P ([y]|Π)general proof :In the limit of τ → ∞ , “Hamiltonian” has aminimal at.

( Hij ≡ )(non-negative deﬁnite)... but in practice, the idea should work with FINITE observation data.
5

General consequence in the case of finite data

general proof :The error of the estimate decreases

proportional to the observation time length.

(ε ≡ τ0 /τ )6

General framework
Estimation of the system parameter: Π
Minimize “Hamiltonian”, Hτ (Π)[y] ≡ −τ −1 ln P ([y]|Π) with respect to Π .
We obtained the general consequence that the framework works.
Estimation of the motion of hidden part: [x]Maximize the path probability, P ([x, y]|Π) with respect to [x] .

(Euler-Lagrange eq.)

(Bayes theorem)

7

General framework
Bayesian inference

Techniques in physics

calculation of the path probability :path-integral approach

P ([x, y]|Π)minimization of the action functionalS([x, y]|Π) with respect to [x]calculation of the marginal
likelihood :�P ([y]|Π) ≡Dx P ([x, y]|Π)(solving the Euler-Lagrange eq. )δS([x, y]|Π) = 0
δ[x]“classical trajectory” [x̂]
WKB (Wentzel-Kramers-Brillouin)
approximation around [x̂]8Experimental data of F1-ATPase40

4000 fps

angle()

30bangle

a20
100-10-20

invisible
-30
-40

0

0.1

0.2

0.3

visible
0.4

0.5

0.6

0.7

time(sec)
timeprobeF1-ATPase

substratum

the data given by T. Nishizaka9eatum

Modelbangleaa) interaction potential between the central subunit

timeprobe
F1-ATPase

invisible
Estimating
parameters:visibleΠ∗ ≡ {k ∗ , h∗ , γ ∗ , Γ∗ }(magenta) and the surrounding subunits (blue andblack), k ∗b) friction coeﬃcient for the rotation of the central∗γ
subunit,substratumc) spring constant of the linker, h∗d) friction coeﬃcient of the probe particle, Γ∗modelbinvisible

visibleﬁxed : γ = k = kB T = 1
varied : 10−2 < Γ∗ < 102 , 10−2 < h∗ < 102∗∗“size of the probe”“stiﬀness of the linker”

10Outline of the numerical experiment
bangleaγ ẋ = −kx − h(x − y) + ξ(t)Langevin equations：Γẏ = −h(y − x) + η(t) timeinvisible

visibleprobe

Under the true parameter values,
trajectories were simulated.

substratum

displacementΓ}estimating parameters： Π = {k, h, γ,F1-ATPasetime

On the basis of [y] alone, Π = {k, h, γ, Γ} were estimatedby the minimization of “Hamiltonian”, Hτ (Π)[y]11displacement

Path probability
P ([x, y]|Π)

yxtimeP ([x, y]|Π) = Pinit (x0 , y0 |Π) Ptr ((x0 , y0 ) → [x, y]|Π)√khinitial distribution : Pinit ((x0 , y0 )|Π) =exp[−βH((x0 , y0 ); Π)]
2πTk 2 hH(x, y; Π) = x + (x − y)22
2transition probability : Ptr ((x0 , y0 ) → [x, y]|Π) =1S([x, y]; Π) ≡4γ

� τ01
dt [γ ẋ + kx + h(x − y)] +4Γ2� τ0� √γΓ
4πT ∆t�Nexp[−βS([x, y]; Π)]T (k + h)τ
T hτ
dt [Γẏ + h(y − x)] −−2γ2Γ212P ([x, y]|Π)displacement

displacement

Marginal likelihoodyxP ([y]|Π) ≡time

�Dx P ([x, y]|Π)time1.) calculate the MAP (Maximum A Posteriori) estimator [x̂]∂lnPinit (x0 , y0 |Π) Ptr ((x0 , y0 ) → [x, y]|Π) = 0
Euler-Lagrange equation:∂xi2.) apply WKB method around [x̂]s([x̂, y]; Π)+Hτ (Π)[y] =kB T

�

φ2 + χ2ln(Γ/4πkB T ∆t) c(x̂0 , y0 ; Π)
−+2
2∆tτ

s([x̂, y]; Π) ≡ τ −1 S([x̂, y]; Π)

φ2 ≡ (k + h)2 /γ 2 , χ2 ≡ h2 /γΓ
13Error & sample averaged varianceoptimized solution set→ the besti : index number of the diﬀerent initial settings of the optimization process
j : index number of the data ofvar. �= 0s.a.var = 0error :ˆvar. �= 0s.a.var �= 0

sample 2→ characteristic parameter ofthe roughness of Hτ (Π)[y]: estimated value

total → ssample 1

sample averaged variance :∗ : true value,

total → rHτ (Π)[y]sample 3

Π14

Estimation results

abcb
10010

0 0

10-1

10101

error
variance

1010
-1

10 -1
10-2
10-2
10-3
10

102
10
2
10 1
10
1010
10
100-1
2
10-1
10
1
2
10-2
100
10
1
10
10-1
10-2
10
100-2-1
10
10
10-2
10-2
10-1

variance
sampleaveragedvariance

10

0
10
100
10-1
10-2-110
10-3-2
10
10-4
10-5

101
0
10 0
10-1
10-1
10-2-3-2
10
10-4
10-3
10-5
10
1022
10

1011
101000
10-1
10-1101022
10110110100010-1
10-1
1010-2-2-2 1010-2
10

The sample averaged variance is abruptly increase at a certain curve.

15

Estimation results
The error is small when

110a

1

1010100probe

1010-1
10-2

100-1

error

102b

F1-ATPase
-2

10-1

100

1

1010

2

10-2

divergent
phase
210b

precisephase

invisible

angle

a

is small.

non-precise
phase

100

visible

10-1

1

10

time

10-2

100

substratum

sample-averagedvarianceandis large

10-3

10-1

10-4-210-5

1010-2

10-1100

1

10102γ ∗ = k ∗ = kB T = 1, ∆t = 10−4 , τ = 102

16

Estimation results
is not enoughto follow the probemotion.
1

10

10-2

100

10-1

10error

10c

non-precise
phase

1

0-1

precise
phase

10-2

100

10-3

10-1

10-4-2

10-5

10

10-2

10-1

10

0

1

10

10

2

100

sample-averaged
variance

10

divergent
phase
102

classiﬁed into 3 distinct regions.

sample-averaged
varianceb

The parameter space can be

10-1
10-2
10-3
10-4
10-5

103

104

iteration number

105

17

τ and ∆t

Dependency on
10210210211

1101010

10

0

100

10

10-1

10-1

10-1-2-2

10-2

1010-210-1

100110102∆t =d10−4

10

10

10-1

100

11010

2

10102

102

102

11

10

10-1

10-2

0-2

10

10-3
10-4-2

10-1

100

11010

2

10-5

100
10-1

1010

10100

100

100

10-1

10-1

10-1

10-4-2-2-2

10-5

10

10-2

10-1

100

1

1010

2

1010-2

10-1

100

11010

2

1010

10-2
10-3-2

10-1

100

11010

2

sample averaged
variancebsample averaged
varianceτ = 1018

Dependency on
τ = 10a101100error

10τ and ∆tb-1

10c∆t = 10−4

1010-510-4

10-3derror

102

101

010The error is almost
independent of.

010-1
10-1

100

1

1010-22-21010102

102

101

101

The error decreases
proportional
to100

10-1-3101

10-2-21010-2100

101010-1

1-2

21010-2-3

10.

1100

--> consistent
with the theory-1-1
1010-1

100

101

102

10-2-2101010-1

100

11010-2
2-21010

119Comparison with the spectrum fitting
the power spectrum of [y] was ﬁtted by the analytical solution :2��2hk+h2ω++γΓ
γ
2kB T�
������
�
S(ω) =hkhΓ
2 − iω k+h + h − kh
ω 2 + iω k+h+−
ωγ
ΓγΓ
γ
ΓγΓ

Our method is more precise than the spectrum ﬁtting up to 100 times.
20

Estimation of the motion of hidden partIs the MAP estimatora good approximationof x(t) ? ba

only in the case ; h /k � 1

Analytical solution:∗angle∗invisible

visibletimeprobe
F1-ATPase

substratumΓ∗ = 0.12 :ﬁxed

Precise estimation of the hidden motion does not required to estimate
the system parameters.

21

Conclusion
We propose the general framework of the parameter estimation methodwith hidden degrees of freedom.

We obtained the general consequence that the error of the estimate decreases as τ−1.

The practical utility of our method is conﬁrmed by the simple model.
- We found “loss-of-precision transition” and “loss-of-estimate transition”.
For precise estimate, the “probe” should be small and tightly coupled to the hiddenpart.

In “precise phase”, the error of the estimate decreases as τwith the theoretical prediction.

−1, which is consistent

If we mistake the choice of the “probe” and the “linker”, we fail to obtain the precise
estimates even if we take the data over long times with super-high time resolution.

- The proposed method is more accurate than the conventional method.
The spectrum method cannot be extended to non-linear systems.
22

Estimation of the motion of hidden part
observed data

the motion of hidden part:

DRP (Dominant Reaction Pathway)

Model A10

potential

Model (Langevin eqs.):

0.5

position
23

MAP estimator in a multiwell potential
Task: estimate the MAP estimatorfrom.

Maximize the path probability

with respect to.

→ Solve Euler-Lagrange eq.:relaxation method : too much computational cost
cf. TDGL (time-dependent Ginzburg-Landau) eq.

Go-and-Back method :& perturbation expansion

solvefromto. (Go)solvefromto. (Back)24

Result
Go-and-Back method yields
reasonable estimates.

100.5

position

20

Model Bpotential

Model A

0.525

ResultGo-and-Back method is faster than the gradient descent ~104 times.
26

Reference & Acknowledgement
Miyazaki, M., Harada, T. (2010) “Bayesian estimation of the internal structure of proteins fromsingle-molecule measurements”, J. Chem. Phys., 133, to appear, November 21, 2010.
!
Miyazaki, M., Harada, T. (2010) “Go-and-Back method: Effective estimation of the hidden
motion of proteins from single-molecule time series”, manuscript in preparation.
Single-molecule techniques
Takayuki Nishizaka (Gakushuin Univ)
Sugawa Mitsuhiro (Gakushuin Univ)
Tetsuaki Okamoto (Kanazawa Univ)

Physics
Shoichi Toyabe (Chuo Univ)
Miki Matsuo-Yamamoto (The Univ of Tokyo)

Bayesian statistics
Tetsuya J. Kobayashi (The Univ of Tokyo)
Ryota Tomioka (The Univ of Tokyo)
特定領域研究「膜超分子モーターの革新的ナノサイエンス」
京都大学GCOEプログラム「普遍性と創発性から紡ぐ次世代物理学」
27