A Fokker-Planck Model in for Two Interacting Populations of Neurons - - PowerPoint PPT Presentation

A Fokker-Planck Model in for Two Interacting Populations of Neurons

J.A. Carrillo1, S. Cordier2, S.Mancini2

1 ICREA, Universitat Aut`

• noma de Barcelona

2 MAPMO, F´

ed´ eration Denis Poisson, Universit´ e d’Orl´ eans

Stochastic Models in Neuroscience CIRM, 18 January 2010

Outline

• The Wilson-Cowan and Deco-Mart

` ı model

• The stochastic ODE system and its Fokker-Planck equation
• Stationnary solution (equilibrium)
• Time dependent problem and Generalized relative entropy
• Numerical simulations
• Slow-fast behavior

Wilson-Cowan model

⊲ Two neurons populations/cells ⊲ Neglect spatial interactions ⊲ Time dynamics

E(t) = excitatory cells prop. activ/unit time at time t I(t) = inhibitory cells prop. activ/unit time at time t

Model : The value of E(t+τ) and I(t+τ) is proportional to number of cells which are sensitive (not refractory) and which also receive at least threshold, θ, excitation at time t.

[WC] H.R. Wilson, J.D. Cowan, Biophysical Journal (1972)

Wilson-Cowan, the model

⊲ Excit./Inhibit. sensitives prop. of cells (r refractory time): 1 −

t

t−r E(t′)dt′ ,

1 −

t

t−r I(t′)dt′

⊲ Sub-pop. prop. receiving at least θ excitation/time, as a fonction of the mean excitation x(t) => response fonction (sigmoide): S(x) =

x

θ D(θ)dθ ,

D(θ) = thresholds distribution ⊲ α(t) decreasing rate of stimuli effect, ci coeff.connexion, Ie(t), Ii(t) stimuli; E(t + τ) =

• 1 −

t

t−r E(t′)dt′

• Se

t

−∞ α(t − t′)[c1E(t′) − c2I(t′) + Ie(t′)]dt′

• I(t + τ) =
• 1 −

t

t−r I(t′)dt′

• Si

t

−∞ α(t − t′)[c3E(t′) − c4I(t′) + Ii(t′)]dt′

• x(t)

Wilson-Cowan, simplified

Time averaging ⇒:

t

t−r E(t′)dt′ → rE(t) ,

t

−∞ α(t − t′)E(t′)dt′ → kE(t)

First order expansion in τ = 0 ⇒ τ dE dt = −E + (1 − rE)Se(kc1E − kc2I + kIe(t)) τ′dI dt = −I + (1 − rI)Si(k′c3E − k′c4I + k′Ii(t)) ⊲ Hysteresis multiple loops and limit cycles. ⊲ Results are independent of the choice of the sigmoide.

The Deco-Mart ` ı model

w+ excitation coeff. between neu- rons of the same population w− excitation coeff. between neu- rons of different populations wI inhibition coeff. between all neu- ronal populations Synaptic force between population i and j: wij =

  

w+ − wI i = j w− − wI i = j ⇒ w11 = w22, w12 = w21 ⊲ Unbiased : λ2 = λ1 ⊲ Biased : λ2 = λ1 + 0.1

[DM] G.Deco, D.Marti, Biological Cybernetics (2007)

The associated stochastic system

⊲ ν1 = ν1(t), ν2 = ν2(t) firing rates of the 2 sub-populations: ⊲ ξ = ξ(t) white noise of amplitude β (brownian motion with variance β2/2).

  

τ ˙ ν1 = −ν1 + φ

• λ1 +

j=1,2 w1jνj

• + ξ

τ ˙ ν2 = −ν2 + φ

• λ2 +

j=1,2 w2jνj

• + ξ

where φ(x) is the sigmoide (response fct.), strictly monotone and bounded: φ(x) = νc 1 + exp(−α(x/νc − 1)), α, νc ∈ R

The Fokker-Planck equation

Let dxt = f(xt)dt + g(xt)dBt a ods syst. and L = fi(x)∂xi + (1/2) dij∂xij the generator

• f T, then u(x, t) = Ttϕ(x) satisfies ∂tu = Lu .

⊲ f(t, ν1, ν2) probability distribution, t ≥ 0 and ν = (ν1, ν2) ∈ Ω

def

= R2

+

⊲ ∂tf + ∇ · (F f) − β2 2 ∆f = 0 (FP) Ff = (−ν + Φ(Λ + W · ν)) f

• Ff − β2

2 ∇f

• · n = 0

⊲ flux incoming F · n ≤ 0 (H1) ⊲ normalization

• Ω f dν = 1

(H2)

No explicit equilibrium solution

Remark : The flux F is not the gradient of a potential A : ∂ν2F1 = ∂ν1F2 ⊲ No explicit solution of the stationnary problem associated to (FP). ⊲ If one proves the existence, uniqueness and positivity of a solution f∞ of the stationnary problem, then it is possible to split F in the sum of a gradient and non-gradient terms: let A = − log f∞ (ie. f∞ = e−A), ∇ ·

• Fe−A − β2

2 ∇(e−A)

• = 0,

hence ∇ ·

     

• F + β2

2 ∇A

• G

e−A

     

= 0 and so: F = −β2 2 ∇A + G Remark : G is s.t. ∇·(Ge−A) = 0, but we do not have an explicit form for G

[ACJ] A. Arnold, E. Carlen, Q. Ju, Communication in Stochastic Analysis (2008)

Stationnary problem

We first consider the stationnary problem associated to (FP) Af = −β2 2 ∆f + ∇ · (Ff) = 0 ,

• Ff − β2

2 ∇f

• · n = 0

(S) Theorem: If (H1) and (H2) hold, then there exists an unique positive solution f∞(ν) to (S).

Proof: Based on the Krein-Rutman theorem :

• T : L2(Ω) → L2(Ω), s.t. ∀g ∈ L2(Ω), Tg = f, with f the unique solution of :

Af + ρf = g in Ω, (Ff − β2 2 ∇f) · n = 0

• n ∂Ω
• T : H2 → H2 is a compact operator, and T : K → K str. pos., with K = W 2,2

+ (Ω).

• KR th. =

⇒ r(T) > 0 and ∃ g > 0 s.t. Tg = r(T)g. We have, Af + ρf = λf, f = r(T)g > 0, λ = 1 r(T) and Af = (λ − ρ)f ⇒ (λ − ρ)

• Ω

f dx = 0 ⇒ ρ = λ ⇒ Af = 0.

Time dependent problem

We consider here the parabolic problem: ∂tf + Af = 0 ,

• Ff − β2

2 ∇f

• · n = 0

(P) with the initial condition: f0(·) ∈ L2(Ω) Theorem: (P) admits an unique solution f(t, ν1, ν2).

Proof: Let a(t, f, g) be the bi-linear form associated to A: a(t, f, g) =

• Ω

β2 2 ∇f · ∇g dν −

• Ω

fF · ∇g dν , ∀ f, g ∈ H1(Ω) , (a)

• a(t, f, g) is continuous,
• a(t, f, g) + ρ < f, g > is coercive for ρ ∈ R large enough.
• Remark : The maximum principle does not apply !

(F has negative divergence ∇ · F ≤ 0)

Generalized relative entropy

Theorem: Let f1, f2 > 0 solutions to (P), and g > 0 a solution of the dual problem:

  

∂tg = −F · ∇g − β2

2 ∆g,

in Ω × [0, T],

∂g ∂n = 0

• n ∂Ω

Then, ∀ H convex d dt

• Ω gf1H(f2

f1 dν

• Hg(f2|f1)

= −β2 2

• Ω gf1H′′(f2

f1 )

• ∇(f2

f1 )

• 2

dν

• Dg(f2|f1)

≤ 0 , (GRE)

Proof:

∂ ∂t [gf1H] = −∇ · [Fgf1H] + β2

2 ∇ ·

• g2∇

f1

g H

• =0 integr. by part

−β2

2 gf1H′′ |∇ (f2/f1)|2

• Lemme: The solution f to (P) satisfies:
• Ω f(t, ν) dν =
• Ω f0(ν) dν.

Proof: Integration of d/dt

Ω fg

• [MMP] P. Michel, S.Mischler, B.Perthame, J.Math. Pures Appl. (2005).

Corollaries

Corollary 1: If f0(ν) > 0, then f(t, ν) > 0 ∀t.

Proof: We choose in (GRE): g = 1, f1 = f∞, f2 solution of (P) with f2(0, ν) > 0, and H(f2/f1) = ε(f2/f1)−. Then: d dt

• Ω

f∞(f2/f1)−dν

• h(t)

≤ 0. h(0) = 0, h(t) is decreasing and h(t) is positive ⇒ h(t) = 0 ∀t

• Corollary 2: If (H1) holds and f is a solution to (P) with initial con-

dition f0, then lim

t→∞

• Ω |f(t, ν) − f∞(ν)|2 dν = 0

Proof: We choose in (GRE): g = 1, f1 = f∞, f2 = f and H(s) = s2/2; and applying the Aubin-Lions theorem

Numerical approximation - finite difference method

Let f k(i, j) = f(k∆t, ni, nj) with ni = (i + 1

2)∆N1, i = 0...N1 − 1 and nj = (j + 1 2)∆N2, j =

0...N2 − 1. Then, the Fokker-Planck equation is discretised by : f k+1(i, j) = f k(i, j) + ∆t

• F k(i + 1/2, j) − F k(i − 1/2, j)
• /∆N1

+ ∆t

• Gk(i, j + 1/2) − Gk(i, j − 1/2)
• /∆N2,

where F k(i + 1

2, j), Gk(i, j + 1 2) are the flux at the interfaces :

F k(i + 1/2, j) =

• −ni+1/2 + Φ(λ + w11ni+1/2 + w12nj)
• f k(i + 1/2, j)

− β2 2∆N1

• f k(i + 1, j) − f k(i, j)
• ,

Gk(i, j + 1/2) =

• −nj+1/2 + Φ(λ + w21ni + w22nj+1/2)
• f k(i, j + 1/2)

− β2 2∆N2

• f k(i, j + 1) − f k(i, j)
• .

and we choose linear interpolation for f at the interfaces: f k(i + 1/2, j) = f k(i + 1, j) + f k(i, j) 2 , f k(i, j + 1/2) = f k(i, j + 1) + f k(i, j) 2 . Remark: Adaptatif ∆t (gain factor 100) ⇒ for i, j s.t. f k(i, j) = 0 and Fk(i, j) = 0: ∆t = min

i,j

f k(i, j) 2|Fk(i, j)|

Computed quantities

Marginales of f(t, ν1, ν2) with respect to ν2 and ν1 : N1(t, ν1) =

νM

f(t, ν1, ν2)dν2 , N2(t, ν2) =

νM

f(t, ν1, ν2)dν1. First order moments: µi(t) =

Ω νif(ν1, ν2, t)dν1dν2,

i = 1, 2 Second order moments: γij(t) =

Ω νiνjf(ν1, ν2, t)dν1dν2,

i, j = 1, 2. Probabilities ρi(t) to belong to the domains Ωi, i = 1, 2, 3 : ρi(t) =

Ωi

f(ν1, ν2, t)dν1dν2. We have N1 = N2 = 200 discretisation points, and stop computations when we get a 10−10 difference between to successive iterations. β = 0.1, α = 4, νc = 20, λ1 = 15, r = 0.3, w+ = 2.35, wI = 1.9, τ = 10−2.

Stationnary solutions

1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10

Contour levels for the distribution f(ν1, ν2) at equilibrium. Left: unbiased, symmetry, concentrations of ”mass” in S1 = (1.3, 5.9) and S3 = (5.9, 1.3). Right: biased, no symmetry, concentration of ”mass” in S1 = (1.1, 6.6) and S3 = (5.5, 1.6)

Convergence to steady state

0.0 0.5 1.0 1.5 2.0 2.5 3.0 −5 −4 −3 −2 −1 1

Convergence towards the steady state, in the unbiased case, in logarithmic scale. Linear regression done on the last quarter of the curve has a slope of -0.08 with a standard deviation of 0.004.

After a small transition period, the convergence of the solution towards its stationary state has an exponential behavior.

Marginales time evolution

initial 0.56 sec 1.15 sec 3.75 sec 1 2 3 4 5 6 7 8 9 10 5 10 15 20 25 30 35 40 initial 0.56 sec 1.15 sec 3.4 sec 1 2 3 4 5 6 7 8 9 10 10 20 30 40 50

Initial condition: Gaussian centered in S = (3, 3), near the instable point S2. Left : unbiased, evolution towards a symmetric double picked distribution in S1 and S3. Right: biased, evolution towards a double picked distribution, a bigger one in S1 and a smaller one in S3.

Probabilities in subdomains Ωi

rho1 rho2 rho3 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 rho1 rho2 rho3 rho1 rho2 rho3 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Domains : Ω1 = [0, 2] × [5, 10], Ω2 = [2, 5] × [2, 5], Ω3 = [5, 10] × [0, 2] initial condition: Gaussian centered near the unstable point S2, dans Ω2. Left: unbiased, time evolution of probablities ρi, ρ1 = ρ3. Right: biased, ρ3 almost equal to zero.

Escaping time

⊲ 1 dim: Kramers law E(T) = exp(H/β2) , H= potential height, β = noise ⊲ dim > 1, no general formulation Let f(0, ν1, ν2) a Gaussian centered in S1 and β = 0.2, ..., 1. Let T be the exit time=time needed to move half of the mass from Ω1 to Ω3: ρ1(T) < 2ρ3(T). We remark that, in our problem, T has an exponential behavior :

β 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 T 12.91 3.33 1.69 1.12 0.80 0.59 0.49 0.37 0.30

Escaping time T as fonction of β in log scaling. Computing times become rather long....

Slow-fast structure of the dynamics

Realization of a trajectory for the ODEs, starting in (5, 5) Blue lines are numerical ap- proximations of the solution of the o.d.e

  

˙ ν1 = −ν1 + φ

• λ +

j=1,2 w1jνj

• ˙

ν2 = −ν2 + φ

• λ +

j=1,2 w2jνj

• they highlight the slow man-

ifold to which belong the equilibrium points

• f

the system ⇒ fast towards the manifold, slow on the manifold

[BG] N.Berglund, B.Gentz, Noise-Induced Phenomena in Slow-Fast Dynamical Systems. A Sample-Paths Approach. Springer, Probability and its Applications (2005)

One-dimensional reduction

The change of variables x = ν1 + ν2, y = ν1 − ν2, in the unbiased case, leads to the following system:

  

ε ˙ x = h(x, y) + ξ ˙ y = g(x, y) + ξ , where ε = ∂yg ∂xh|S2 Let x∗(y) be the solution of h(x, y) = 0, then we can write the problem in a unidimensional equation defined on the slow manifold: ˙ y = g(x∗(y), y) + ξ Again, we associate an one-dimensional Fokker-Planck equation. Now, the stationnary solution is given explicitely by: exp

• −G

β2

• ,

G = ∂yg(x∗(y), y) In this case, we can say exact things on the escaping time too.

Comparison : 2D vs. simplified model

Plot of the y-marginal of the 2D equilibrium state and the stationary solution of the reduced slow fast system in the unbiased case and β = 0.3

Comparison : simplified model vs. moment system

Left : ν1-marginal obtained from the Deco and Mart ´ ı moment system, for various value of β in [0.1, 0.6] in the biased case. Right : y-marginal obtained for the reduced system in similar conditions.

Conclusions - Perspectives

• We propose a kinetic model for the time evolution of two interacting

populations (of neurons), derived from a system of Wilson-Cowan type. This model describes, statistically, both equilibrium and evolution.

• We present theoretical results about existence, uniqueness, positivity and

trend to equilibrium for the F-P equation, which as a non-potential flux.

• Study (in progress) of the simplified 1D model using the slow-fast re-

duction, in agreement with [Deco et Mart ` ı] and 2D model for escaping time, trend toward equilibrium