M ethode de Monte-Carlo et marche sur les sph` eres pour l - - PowerPoint PPT Presentation

m ethode de monte carlo et marche sur les sph eres pour l
SMART_READER_LITE
LIVE PREVIEW

M ethode de Monte-Carlo et marche sur les sph` eres pour l - - PowerPoint PPT Presentation

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation M ethode de Monte-Carlo et marche sur les sph` eres pour l equation de Poisson-Boltzmann lin eaire et non-lin


slide-1
SLIDE 1

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation

M´ ethode de Monte-Carlo et marche sur les sph` eres pour l’´ equation de Poisson-Boltzmann lin´ eaire et non-lin´ eaire de la dynamique mol´ eculaire

Nicolas Champagnat, joint work with

  • M. Bossy, D. Talay, L. Violeau, M. Yvinec (Inria Sophia),
  • H. Leman (CMAP, X), S. Maire (Univ. Toulon)

S´ eminaire de probabilit´ es, LJK, Grenoble, 20 novembre 2014

slide-2
SLIDE 2

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

Poisson-Boltzmann PDE of molecular dynamics

The Poisson-Boltzmann (PB) PDE of molecular dynamics is solved by the electrostatic potential around a biomolecular assembly, which is used to compute global characteristics of the system like

  • the solvatation free energy
  • the electrostatic forces exerted by the solvent on the molecule.

It takes the form of the Poisson equation div("(x)ru(x)) + 2(x) sinh u(x) = f (x), x 2 R3, where

  • "(x) is the permittivity of

the medium

  • 2(x) is called the “ion

accessibility parameter”. This is an “implicit solvent” equation, which means that the solvent is considered as a continuum.

slide-3
SLIDE 3

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

The geometry of the problem

The atomic structure of the molecule is modelled in details:

  • N atoms at positions

x1, . . . , xN in Ωi with radii r1, . . . , rN and charge qi

  • Ωi = [N

i=1B(xi, ri)

slide-4
SLIDE 4

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

Linear and non-linear PB equations

div("(x)ru(x)) + 2(x) sinh u(x) = f (x), x 2 R3, (PB) When the source term is not too big (i.e. when the biomolecule is not very charged), the electrostatic potential can be approximated by the solution of the linearized Poisson-Boltzmann equation div("(x)ru(x)) + 2(x)u(x) = f (x), x 2 R3, (LPB)

slide-5
SLIDE 5

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

Difficulties

  • The source term is singular

f :=

N

X

i=1

qixi. This difficulty can be removed by considering the solution of "i∆u = f u0(x) = 1 4⇡"int

N

X

l=1

qj |x xl| 8x 2 Ωint. Then v := u u0 solves PB equation with a smooth source term, if has compact support in Ωi and ⌘ 1 in the neighborhood of {x1, . . . , xN }.

  •  is discontinuous
  • The operator has divergence form with discontinuous coefficient ".
  • Nonlinearity.
slide-6
SLIDE 6

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

Arguments in favour of a probabilistic resolution method

  • Dimension 3
  • The free energy of the solution of this PDE has the form

1 2

N

X

i=1

qi(u u0)(xi)

  • Walk on spheres techniques allow us to deal with the geometry of

the molecule efficiently

slide-7
SLIDE 7

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Introduction

Arguments in favour of a probabilistic resolution method

  • Dimension 3
  • The free energy of the solution of this PDE has the form

1 2

N

X

i=1

qi(u u0)(xi)

  • Walk on spheres techniques allow us to deal with the geometry of

the molecule efficiently This talk deals with two aspects of the problem:

  • the probabilistic interpretation of the divergence-form operator

L := div("r) (in any dimension)

  • the construction of Monte Carlo approximation algorithms based
  • n walk on sphere, appropriate discretization, branching

Brownian motion and adapted to the geometry of the problem

slide-8
SLIDE 8

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation What is known in dim 1?

The case of dimension 1

The case of dimension 1 is quite well understood (Portenko 1979, Le Gall 1985, Ouknine 1990, several works of Etor´ e, Lejay, Martinez...). Consider the case "(x) = ( " if x < 0 "+ if x 0 Then

  • L is the infinitesimal generator of the solution of the SDE

dXt = p 2"(Xt)dBt + "+ " 2"+ dL0

t (X ).

  • There is strong existence and pathwise uniqueness for this SDE.
slide-9
SLIDE 9

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation What is known in dim 1?

The case of dimension 1 (2)

  • Two linear scalings of R+ and R transform X into a skew

Brownian motion of parameter =

pε+pε pε++pε .

  • Another piecewise linear transformation transforms X into a

standard SDE with discontinuous diffusion coefficient.

  • Exact simulation schemes are known for skew the Brownian

motion:

  • this is a standard Brownian motion until it hits 0
  • the probability to hit h before −h starting from 0 is (β + 1)/2 and

the hitting time has the same law as for the Brownian motion.

slide-10
SLIDE 10

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation multidimensional SDE with local time

Notation

We assume that Γ is a smooth (C 3) manifold in Rd. We introduce the notations

  • ⇡(x) for the orthogonal

projection of x on Γ

  • n(y) as the outward normal to

Γ for y 2 Γ

  • ⇢(x) as the signed distance

between x and Γ ⇢(x) := (x ⇡(x)) · n(⇡(x)).

slide-11
SLIDE 11

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation multidimensional SDE with local time

A martingale problem

Smooth functions in the domain of L must satisfy the transmission property "intrint'(x) · n(x) = "extrext'(x) · n(x). We say that (Px)x2Rd on (C, B) solves the martingale problem (MP) for L if, for all x 2 Rd, one has Px{w 2 C : w(0) = x} = 1, and, for all ' satisfying ' 2 C 0

b (Rd) \ C 2 b (Rd \ Γ),

"r' · (n ⇡) 2 C 0

b (N),

  • ne has

M ϕ

t (w) := '(w(t)) '(w(0))

Z t L'(w(s))ds is a Px martingale.

slide-12
SLIDE 12

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation multidimensional SDE with local time

The main result

Theorem (Bossy, C., Maire, Talay, 2010) The martingale problem for L is well-posed. In addition, there is weak existence and uniqueness in law for the SDE 8 < : dXt = p 2"(Xt)dBt + "ext "int 2"ext n(Xt)dL0

t (Y ),

Yt = ⇢(Xt), (1) and the corresponding laws solve the MP for L.

slide-13
SLIDE 13

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation multidimensional SDE with local time

Ingredients of the proof

  • We construct a smooth local straightenings of Γ defined on a

neighborhood U of x, s.t. 1 = ⇢ and, if (Xt) solves (1), then Zt := (Xt) satisfies dZ 1

t =

p 2"(Xt)dBt + "e "i 2"e dL0

t (Z 1) + (drift)(Zt)dt,

and there is no local time term in the SDE solved by Z 2

t , . . . , Z d t .

  • Girsanov’s formula allows one to remove the drift term, so that

Z 1

t solves a one-dimensional SDE

  • The SDEs for Z 2, . . . , Z d can then be solved as time dependent

SDEs for all paths of Z 1

  • Only weak existence.
slide-14
SLIDE 14

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation multidimensional SDE with local time

Ingredients of the proof (continued)

Lemma (Generalized Itˆ

  • -Meyer formula)

If X is a continuous semimartingale and if u is a function satisfying the assumption of the MP for L, then u(Xt) = u(X0)+ Z t rintu(Xs) · dXs+1 2

3

X

i,j=1

Z t @2u @xi@xj (Xs)dhX i, X jis + 1 2 Z t f (Xs)dL0

s(Y ),

8t 0 a.s., where f (x) = ✓"int "ext 1 ◆ rintu(⇡(x)) · n(⇡(x)).

  • The formula would be easily obtained from Itˆ
  • ’s and Itˆ
  • -Tanaka’s

formulas for a functions of the form u(x) f (x)[⇢(x)]+ approximation argument.

  • If (Xt) solves (1), the local time terms cancel.
slide-15
SLIDE 15

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Feynman-Kac formulas

Back to linearized Poisson-Boltzmann equation

Proposition (Classical Feynman-Kac) Let v be the solution of div("rv) + 2v = g, where g is a smooth

  • function. Then, for all x 2 R3,

v(x) = Ex Z +1 g(Xt) exp ✓

  • Z t

2(Xs)ds ◆ dt

  • .

This representation is not very convenient for simulation since

  • One needs to precisely discretize X everywhere where g is

nonzero.

  • In general, the computation of g is costly.

Since X has (scaled) Brownian paths away from Γ, it is better to have formulas only involving informations on the entrance time and position in small neighborhoods of Γ.

slide-16
SLIDE 16

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Feynman-Kac formulas

Other Feynman-Kac formulae

Fix h > 0.

  • We define the stopping time ⌧ = inf{t 0 : Xt 2 Γ}
  • Since ∆(u u0) = 0 in Ωi, for all x 2 Ωi,

u(x) = Ex[u(Xτ) u0(Xτ)] + u0(x).

  • For all x 2 Ωe,

u(x) = Ex h u(Xτ1)e¯

κ2τi

.

  • For all x 2 Γ, using the transmission property

"intrint'(x) · n(x) = "extrext'(x) · n(x). and the relations hrintu(x) · n(x) = u(x) u(x hn(x)) + O(h2) and hrextu(x) · n(x) = u(x + hn(x)) u(x) + O(h2), we get u(x) = "int "int + "ext u(x hn(x)) + "ext "int + "ext u(x + hn(x)) + O(h2).

slide-17
SLIDE 17

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Feynman-Kac formulas

Jump on the boundary of the molecule

The last equation u(x) ⇡ "int "int + "ext u(x hn(x)) + "ext "int + "ext u(x + hn(x)) can be interpreted as the expectation of u(x + ⇠hn(x)) where ⇠ is a r.v. in {1, 1} equal to 1 with probability

εext εint+εext

jump of the particle when it hits Γ.

slide-18
SLIDE 18

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Feynman-Kac formulas

The Feynman-Kac formula

  • We define the stopping times

⌧k = inf{t ⌧ 0

k1 : ⇢(Xt) = h}

⌧ 0

k = inf{t ⌧k : Xt 2 Γ}.

  • A recursive application of the last steps suggests the following

formula: Proposition (Bossy, C., Maire, Talay, 2010) u(x) = Ex "+1 X

k=1

  • u0(Xτk ) u0(Xτ 0

k )

  • exp

  • Z τk

2(Xt)dt ⌘# .

slide-19
SLIDE 19

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Feynman-Kac formulas

Algorithm (Mascagni and Simonov 2004)

8x 2 Ωi, u(x) = Ex[u(Xτ 0

1) u0(Xτ 0 1)] + u0(x),

(2) 8x 2 Ωi, u(x) = Ex h u(Xτ1)e¯

κ2τi

, (3) 8x 2 Ωi, u(x) ⇡

εint εint+εext u(x hn(x)) + εext εint+εext u(x + hn(x)). (4)

Idea of the algorithm: starting from x0 2 Ωi (say),

1 simulate a Brownian path starting from x0 until it hits Γ at x1

and score u0(x0) u0(x1) (by (2), this is OK in expectation),

2 jump at a distance h inside with probability εint εint+εext or outside

with probability

εext εint+εext (by (4), this is OK in expectation), 3 if the new position x2 is inside, continue as in step 1, otherwise,

go to next step,

4 simulate a Brownian path starting form x2, compute its hitting

time ⌧ of Γ, interpret e¯

κ2τ as a probability of killing the

  • particle. If the particle is not killed, continue as in step 2,
  • therwise end the algorithm (by (3), this is OK in expectation).
slide-20
SLIDE 20

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Basic Walk on Spheres algorithm

Walk on Spheres algorithm (1)

How to simulate the Brownian path outside the molecule and the probability of killing? WOS algorithms are efficient techniques for this problem.

slide-21
SLIDE 21

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Basic Walk on Spheres algorithm

Walk on Spheres algorithm (2)

The exponetial terms in the FK representation can be interpreted as probabilities of killing the particle in each step of the WOS algorithm: if t0 is the first exit time of the Brownian motion out of B(0, r0), we have E(exp(t0)) = r0 p 2 sinh(r0 p 2) , 8 0.

slide-22
SLIDE 22

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Basic Walk on Spheres algorithm

Walk on Spheres algorithm (3)

WOS algorithm: given y0 2 Ωe, := ¯ 2/2"e and " > 0, fix k = 0

1 Let S(yk, rk) be the largest open sphere included in D centered

at yk.

2 Sample yk+1 according to the uniform distribution on @Sk. 3 Kill the particle with probability 1 rk

p 2/ sinh(rk p 2), and END if killed.

4 IF d(yk+1, Γ)  ", THEN set exit(y0) as the closest point of Γ

from yk+1 and END. ELSE, set k = k + 1 and return to Step 1.

slide-23
SLIDE 23

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Basic Walk on Spheres algorithm

Uncentered Walk on Spheres algorithms

Inside the molecule:

  • since  = 0, we can use a WOS using only the spheres

corresponding to atoms of the molecule.

  • Even if the initial position is not the center of the sphere, the

distribution of the exit position is explicit: Poisson kernel r 2 |y|2 4⇡r|y x|3 d(y).

  • exact simulation of the first exit point out of the molecule.
slide-24
SLIDE 24

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Basic Walk on Spheres algorithm

The algorithm for (LPB)

Given x0 2 Ωint, set k = 0 and score = 0.

1 IF xk 2 Ωint, 1 THEN use the UWOS algorithm to simulate exit(xk) and set

score = score − u0(exit(xk)),

2 ELSE use the WOS algorithm to simulate exit(xk).

IF the particle has been killed, END.

2 Let xk+1 = exit(xk) + hn(exit(xk)) w.p. εext εint+εext and

xk+1 = exit(xk) hn(exit(xk)) otherwise.

3 IF xk+1 2 Ωint, THEN set score = score + u0(xk+1). 4 Set k = k + 1 and return to Step (1).

This algorithm can be improved using other replacement procedures

  • f higher order or better numerical cost (Bossy, C., Maire, Talay,

2010, Lejay, Maire, 2012, Maire N’guyen, 2014...).

slide-25
SLIDE 25

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

An example The molecule. Initial score= 0.

slide-26
SLIDE 26

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

An example score = score + u0(xk)

slide-27
SLIDE 27

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

An example score = score u0(xk+1)

slide-28
SLIDE 28

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

An example score = score + u0(xk+1)

slide-29
SLIDE 29

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

An example End of the algorithm.

slide-30
SLIDE 30

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation A picture of the algorithm

Bias estimates

If Γ is smooth enough, we have the following error estimate: Proposition Denoting by ¯ uh(x0) the expectation of the score given by any of the algorithm starting from x0, we have sup

x02Ωint

|¯ uh(x0) (u u0)(x0)| = O(h + "/h), Higher orders of convergence can be obtained using more precise jump distributions (Lejay-Maire, 2012, Maire-N’Guyen, 2014...).

slide-31
SLIDE 31

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical test on a true molecule

Practical implementation on a true bio-molecule [BCLMVY, 2014]

  • For an efficient Monte-Carlo computation, one needs to be able

to construct the largest sphere in the WOS algorithm efficiently, i.e. the closest atom from the position of the particle.

  • It is classical (in algorithmic geometry) to solve this problem

when the Euclidean distance is replaced by the power-distance (squared distance) Edelsbrunner’s power diagram allows to find efficiently the atom with smallest power.

  • To find the closet atom, one can simply look at the other close

atom with small power.

  • All the tools to do that are implemented in the CGAL library

(C++).

  • Parallelisation is easy.
slide-32
SLIDE 32

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical test on a true molecule

Comparison of the localisation methods

slide-33
SLIDE 33

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical test on a true molecule

Comparison with APBS Algorithm

APBS is a free software developped by Baker and McCammon which solves Poisson-Boltzmann equations with deterministic methods. —: Resolution with APBS algorithm —: Resolution with WOS algorithm (N = 40000, h = 0.1, ✏ = 104)

slide-34
SLIDE 34

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical test on a true molecule

Numerical test: order of convergence

Mean error vs. h (log scales). Computations are done with ✏ = 106, N = 106 and values for h ranging from 0.9 to 0.003, on a bio-molecule composed of 103 atoms.

slide-35
SLIDE 35

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical test on a true molecule

Numerical test: performance

Performance (error vs. CPU time). Parameters are ✏ = 106, N = 106 and values for h ranging from 0.9 to 0.003, on a bio-molecule composed of 103 atoms.

slide-36
SLIDE 36

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

The branching Brownian motion in Ωext

in Ωext, "ext∆u(x) + ¯ 2 X

k0

1 (2k + 1)!u2k+1(x) = 0, Consider a Brownian particle starting (Xt, t 0) from x0 2 Ωext and an independent exponential time ⌧ of parameter := ¯ 2/2"ext. Let ⌧Γ the first hitting time of Γ by X .

  • If ⌧Γ  ⌧, the particle is stopped at XτΓ
  • If ⌧ < ⌧Γ, the particle is killed at time ⌧ and replaced by K

particles located at Xτ, where P(K = 2k + 1) = 1 (2k + 1)!, k 1, P(K = 0) = 2 sinh 1 > 0.

  • The new particles behave independently of each other and follow

the same dynamics as the first one.

slide-37
SLIDE 37

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

Subcritical BBM

Note that the mean number of children of a particle is X

k1

1 (2k + 1)!(2k + 1) = cosh 1 1 < 1. Therefore, the branching Brownian motion has a.s. a finite number of particles.

slide-38
SLIDE 38

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

Notations

We use the notation of Neveu-Harris-Ulam for the individuals of the population:

  • Individuals are elements of H = {v 2 [n0Nn}.
  • The ancestor is denoted by ; (= N0).
  • The jth child of an individual v 2 H is denoted by vj.

For each individual v 2 H, we denote by

  • ✓v and v its time of birth and death (✓v = v = +1 if the

individual never lived),

  • Kv its number of children,
  • X v

t its position for ✓v  t < v.

slide-39
SLIDE 39

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

FK formula with branching diffusion

Theorem (BCLMVY, 2014) Let ⌧n be the first time where the number of births in the BBM is n. Then, for all x 2 Ωext, Ex(Yτn) = u(x), (5) where Yt = Y

v2H

h

  • {σvt, X v

σv 62Γ, Kv1}

+u(X v

σv ) {σvt, X v

σv 2Γ} + u(X v

t ) {θvt<σv}

i In particular, if Y1 has a finite expectation, u(x) = Ex 8 < :(1)#{v:X v

σv 62Γ}

Y

v:σv6=+1

h

X v

σv 62Γ

Kv1 + X v

σv 2Γ u(X v

σv )

i 9 = ; .

slide-40
SLIDE 40

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

Comments on the formula

u(x) = Ex 8 < :(1)#{v:X v

σv 62Γ}

Y

v:σv6=+1

h

X v

σv 62Γ

Kv1 + X v

σv 2Γ u(X v

σv )

i 9 = ; .

  • As soon as one of the term is 0, i.e. as a particle dies in Ωe with

no children, the quantity in the expectation is 0.

  • The expectation might be undefined if u is too big on Γ. It is

finite for example if |u(x)|  1 for all x 2 Γ.

  • The variance of this method can be also related to a non-linear

PDE in Ωext. Again, the variance can be infinite if u is too big on Γ.

  • Other probabilistic interpretations with different distribution of

children, rate of branching... are possible.

slide-41
SLIDE 41

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Branching diffusions

Idea of the proof

True proof: write the semimartingale decomposition of Yt (it is a martingale) and apply the stopping theorem. Intuitive proof (standard argument for branching processes): let ⌧ be the first branching time and h a small parameter, and apply the Markov property at time ⌧ ^ h to get u(x) = Ex h

  • σ;h

K;1 u(X ; σ;)K;

i + Ex h

σ;>h u(X ; h )

i = (1 eλh) X

k1

pkE[ˆ u(x + B;

E (h))k]

+ eλh ˆ u(x) + 1 2 Z h E[∆ˆ u(x + B;

s )]ds

! + o(h). Assuming u is smooth and arranging the terms we obtain when h ! 0, 1 2∆u(x) + u(x) X

k0

pkg(k)u(x)k = 0, x 2 Ωext,

slide-42
SLIDE 42

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation The algorithm

The algorithm

This leads to a similar algorithm as the first one, where

  • the motion of particles in Ωext is replaced by that of the

branching Brownian motions.

  • the score of a Monte Carlo run must be computed recursively :

starting from particles with the highest generation, the scores of brother particles must by multiplied and substracted to the score

  • f their mother particle.
slide-43
SLIDE 43

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation The algorithm

An example

slide-44
SLIDE 44

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation The algorithm

An example

slide-45
SLIDE 45

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation The algorithm

An example

slide-46
SLIDE 46

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical tests

Implementation

  • Since the distribution of children is independent of the position of

the particle, we

  • sample first the Galton Watson tree representing the genealogy of

the particles;

  • then simulate the motion of each particle (corresponding to each

branch of the GW tree) and compute their individual score (forward in time);

  • then compute the total score of the run by computing the product
  • f the scores of particles along the tree, backward in time.
  • Because of the potential variance problems when the solution u is

too large on Γ, we used a variance reduction technique: stratified sampling w.r.t. the Galton-Watson genealogical tree.

slide-47
SLIDE 47

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical tests

Two spheres: convergence

Mean error (log scale) for two atoms of radius 1 ˚ A at a distance of 0.2 ˚ A, with opposite unit charge. N = 105, ✏ = 105 and the values of h range from 0.003 to 0.9.

slide-48
SLIDE 48

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Numerical tests

Two spheres: variance increases with distance

Mean error (log scale) for two atoms of radius 1 ˚ A at a distance of 0.5 ˚

  • A. N = 105, ✏ = 105 and the values of h range from 0.003 to 0.9.
slide-49
SLIDE 49

Introduction Probabilistic interpretation of L MC algorithms for (LPB) MC algorithm for the nonlinear PB equation Conclusion

Conclusion

  • The performance of the linear algorithm is good and scales well

with the number of atoms of the bio-molecule.

  • The non-linear algorithm still suffers from variance problems

appropriate variance reduction techniques need to be tested (tuning of free parameters, importance sampling to avoid too big variance within each stratum, pruning of the Galton-Watson tree...)