Particle systems with geometrical constraints
and applications
Antoine Diez
Kyoto University Institute for the Advanced Study of Human Biology (ASHBi)
Joint work with J. Feydy (Inria)

French-Japanese Laboratory of Mathematics and its Interactions
The University of Tokyo
23rd January 2024
1 / 29Collective motion of active particle systems(a) Crowd(d) Bacteria (E. Coli)[YouTube, M. Copeland](b) Bird flock[Wikimedia](e) Cell aggregation[Glover et al., 2017](c) Ant trails[Perna et al., 2012](f) Data sampling[Clarté, D., Feydy, 2022]
2 / 29

The kinetic theory framework

The Boltzmann equation: ∂t ft (x, v) + v · ∇x ft = Q(f, f ) (1872)
Hilbert 6th problem (1900)
Developing mathematically the limiting processes [. . . ] which leadfrom the atomistic view to the laws of motion of continua.

Microscopic scale.
Mesoscopic scale when N → +∞.
N identical particles in a space E. ft ∈ P(E) distribution of a typical particle.
N d-dimensional dynamical system. Compute the evolution of statistical quantities
3 / 29

Stochastic Mean-Field models
Let Xti , . . . XtN , be N symmetric E-valued stochastic processes initially i.i.d.
Foundations of Kinetic Theory (1956)
Probabilistic interpretation of the Boltzmann equation and themathematical notion of propagation of chaos:Mark Kac• ∀t > 0, ftk,N (dz1 , . . . , dzk ) → ft (dz1 ) . . . ft (dzk )PN
1• ∀t > 0, µ̂Nt := Ni=1 δXti → ft
Propagation of chaos for a class of non-linear parabolic
equations (1967)

Henry P. McKean

Extension of Kac ideas to diffusion in Rd and other stochastic
particle models, in particular:N
1 X
dXti =
K(Xti , Xtj )dt + dBti .
N j=1


Z
√
∂t ft (x) = −∇ · ft (x)K(x, y)ft (y)dy + 2∆ft (x)Rd

Review: [Sznitman, 1989], [Méléard,1996], [Chaintron, D., Kinet. Rel. Models. 15, 2022]
4 / 29

Particles with a body-orientation(a) Aligned(b) Upside-down(c) [Woolley et al., 2009]−→ In this case the phase-space is the manifold E = R3 × SO3 (R)Structure of the space of rotation matrices SO3 (R)• A ∈ SO3 (R) ⇐⇒ AAT = I3 and det A > 0.
• SO3 (R) is a submanifold in the space of 3 × 3 matrices.
• SO3 (R) is endowed with the dot product A · B := 12 Tr(AT B).
• SO3 (R) is a compact group endowed with the Haar measure
R
SO3 (R)R
R
φ(A)dA= SO (R) φ(BA)dA= SO (R) φ(AB)dA,
3

3

for all B ∈ SO3 (R) and measurable function φ.
5 / 29

Alignment process for body-oriented particles
Constant speed motion
dXti = c0 Ait dt

 P
√
N
1
dAit = PTAi ◦ N j=1 K(|Xtj − Xti |)Ajt dt + 2 dBtit

Projection on the tangent plane

Alignment

Matrix-valued Brownian motions

Each particle is represented by acolored tetrahedron

6 / 29Main result for body oriented particles
Theorem ([Degond et al, 2017], [D. 2020], [Degond, D., Frouvelle, Merino-Aceituno, 2020])
When N → +∞, there is propagation of chaos towards
!
Z
′
′
′
∂t ft (x, A) + Ae1 · ∇x ft = −∇A · ft PTA
A ft (x, A )dA + ∆A ft .
SO3 (R)On long space-time scales, the spatial density ρ and average body-orientation
A = (Ω, u, v) ∈ SO3 (R) (moments of ft ) satisfy:∂t ρ + ∇x · (c1 ρΩ) = 0,ρ(∂t + c2 (Ω · ∇x ))Ω + PΩ⊥ (c3 ∇x ρ + c4 ρr) = 0,ρ(∂t + c2 (Ω · ∇x ))u − u · (c3 ∇x ρ + c4 ρr)Ω + c4 ρδv = 0,ρ(∂t + c2 (Ω · ∇x ))v − v · (c3 ∇x ρ + c4 ρr)Ω − c4 ρδu = 0,
where c1 , c2 , c3 , c4 are explicit functions of ρ and c0 andr := (∇x · Ω)Ω + (∇x · u)u + (∇x · v)v,
δ := [(Ω · ∇x )u] · v + [(u · ∇x )v] · Ω + [(v · ∇x )Ω] · u.
−→ Phase transitions, travelling wave solutions . . .
7 / 29A more complicated problem: particles with a shape
In biology, particles cannot be reduced to a point. . .

Objectives(1) Modelling shapes: an optimal transport approach(2) In silico modelling: implementing biophysical laws
(3) Large-scale phenomena: a new class of PDE models (in progress!)
Assumption
At a given time, there are N particles with prescribed volumes v1 , . . . , vN .
8 / 29

1. Modelling shapes(Static framework)

9 / 29

Some existing computational (and analytical) models
1. Point-particle system with soft constraints or parametrized shapes
Relatively easy to implement (up to N = 106 ) and tocoarse-grain but little control on shape deformations.
2. Boundary deformation: Level-Set method, curvature flow. . .
Mostly “one-particle”, theoretically and numerically difficult.
3. (Multi) Phase-Field and free boundary methods: Allen-Cahn equation
Detailed description but based on soft constraints andexpensive for N ≫ 1 with many technical parameters.
4. Cellular automata: Potts model, [CompuCell3D], [Morpheus] . . .
Easy to implement with hard constraints but fully discretemodel with no obvious continuum limit.
5. Vertex models: [Tissue Forge], [VirtualLeaf ]. . .
Discrete model mostly appropriate for crowded situations
such as tissues.
6. Finite Elements: [NVIDIA OmniverseTM ]. . .
Fast and physically accurate, but it may need to be adapted
to biological systems (small-scale physics).
10 / 29

Optimal transport: discrete caseTransport problem:• N target points xi , M source points yj
PP• Equal mass:i αi =j βj = 1• Symmetric transport cost: e.g. c(yj , xi ) = |yj − xi |2PP• Discrete measures: α̂ =αi δxi , β̂ = βj δyj−→ Minimize the total cost. . .
Monge problem

Kantorovich problem

Matching N = M , αi = βj = N1 ,M
XGeneral case: find π = (πij ) such that
N X
M
N
M
X
X
X
inf
πij c(yj , xi ),πij =πij = 1.

inf

1σ∈SM Mc(yj , xσ(j) )j=1−→ Measure theory:
Z

 inf
c(y, T (y))β̂(dy)TE
T #β̂ = α̂, i.e. α̂(B)=β̂(T −1 (B))π

i=1 j=1

i=1

j=1−→ Measure theory: Wasserstein distanceZ Wc (α̂, β̂) := inf
c(x, y)π(dx, dy)
π



E×E

π1 (dy) := π(E × dy) = β̂, π2 (dx) = α̂Note: Kantorovich reduces to Monge when πij = δi,σ(j) or π = (Id × T )#β̂.

11 / 29

Optimal transport: continuous and semi-discrete cases
Continuous probability densities.
Let µ, ν be two probability measures on Rd .
Theorem (Brenier, Gangbo, McCann. . . )Let c(x, y) = ℓ(y − x) with ℓ strictly convex and C 1 .
If the source measure has a density w.r.t. theLebesgue measure, then the Monge problem has aunique solution T : Rd → Rd , T #µ = ν.
Semi-discrete optimal transport
PN
Let Leb be the Lebesgue measure on Ω ⊂ Rd and µ̂ = i=1 vi δxi be discrete.
Lemma
Any transport map T : Rd → Rd such thatT #Leb = µ̂ satisfies:• T (y) ∈ {x1 , . . . , xN } for all y ∈ Rd .
• |T −1 ({xi })| = vi−→ Transport is equivalent to a partition of the space with volume constraint.
12 / 29

Laguerre tessellation
Theorem [Kitagawa, Mérigot, Thibert, 2016], [Bourne, Schmitzer, Wirth, 2018]
Assume that c(·, xi ) ∈ C 1,1 (Ω) for all i. Then the semi-discrete Monge problem
PN
with source Leb and target µ̂ = i=1 vi δxi has a unique solution given by thetransport map T : Ω → {xi }i , x ∈ Li (w∗ ) 7→ xi where the sets
Li (w∗ ) = T −1 ({xi }) = {x ∈ Ω : c(x, xi ) − wi∗ ≤ c(x, xj ) − wj∗ ∀j}are called the Laguerre cells and the potentials wi∗ ∈ R are uniquely defined sothat they form a partition of Ω such that for all i ∈ {1, . . . , N }|Li (w∗ )| = vi .
Note that the partition {Li (w∗ )}i is optimal in the sense that:N Z
X
Wc (Leb, µ̂) = minc(x, xi ) dx.
∪Ui =Ω|Ui |=vi i=1

Ui−→ For a given set of points xi , the Laguerre tessellation is a generalized Voronoi
diagram with a volume constraint. In comparison, classical Voronoi partition
assumes c(x, y) = |x − y| and wi∗ = 0. It also provides an optimal transport pointof view on weighted Voronoi diagrams and power diagrams.
13 / 29

Examples: boundary and shape
Boundary between two Laguerre cells:

The extra void space: xN +1 with{x ∈ Ω, c(x, xi ) − wi∗ = c(x, xj ) − wj∗ }c(x, xN +1 ) = 0−→ The base shape is a level-set ofthe cost function.
−→ For any arbitrary shape definedby a polar curve r = r(θ),
c(y, x) =c(x, y) = |y − x|2
Straight lines

Bubbly diagram[Hayashi et al., 2022][Ishimoto, Morishita 2014]|y − x|r θ(y − x)

14 / 29

Implementation on a discrete grid−→ Given µ̂ =PN

i=1 vi δxi , how to solve fast the dual semi-discrete OT problem(Zc(x, xi ) − wi dx +Wc (Leb, µ̂) = supw∈RN

Li (w)

NX)
wi v i?

i=1▶ c(x, y) = |x − y|2 , fast (CPU) [Bourne et al., 2020], [sd-ot], [Kitagawa et al., 2016]. . .
▶ For arbitrary costs (and in 3D), the computation time explodes. . .
Solution (thanks to Jean Feydy): discrete-semi-discrete optimal transport
Simply replace the continuous measure by a discrete sample
M

Leb ≈

1 X
δy ,
M j=1 jM ≫ N,

Approximate Monge map T̃ : yj 7→ xj ∗ wherej ∗ := argmin {c(yj , xi ) − wi }.
i∈{1,...,N }and solve the fully discrete Kantorovich problem using GPU routines, namely an
adapted Sinkhorn algorithm [Feydy at al., GeomLoss], [Feydy et al., AIStat, 2019].
15 / 29

The original idea: material science and steel industry
Steel is composed of small grains (“crystals”):−→ The grain geometry affect the mechanical properties of the material.
−→ Construct mathematically a new material by specifying the grains size, shape,
positions etc and study its properties in silico.
“Geometric modelling of polycrystalline materials:
Laguerre tessellations and periodic semi-discrete
optimal transport”[Bourne et al., Mech. Res. Comm, 2023][Alpers et al., Phil. Mag., 2015]16 / 29

2. Implementing biophysical laws
(Dynamical framework)
1. Incompressibility, volume exclusion, deformability
2. Contact interactions: pressure, surface tension
3. Deformation-driven motion

17 / 29

Dynamical setting: incompressibility force
For any positions xi of the seeds, there exists a unique Laguerre tessellation.
−→ Introduce the incompressibility force: [Xin et al., SIGGRAPH Asia 2016]Zẋi = −∇xi Wc (Leb, µ̂) = −∇xi c(y, xi )dy
Li

Similar to a repulsion force, but which is usually a sum of binary interactions:ẋi = −NX

∇xi V (xj −xi ),

j=1

which does not take into account the global shape configuration of the particles.
Example: c(x, y) = |y − x|2−→ Convergence towards the preferredshape (disk) in order to minimize
Wc (Leb, µ̂) = min

N Z
X∪Ui =Ω|Ui |=vi i=1c(x, xi ) dx.

Ui
18 / 29

Particle systems with arbitrary shapes
The incompressibility force can be added to any classical particle system:ẋi = F⃗ − ∇xi Wc (Leb, µ̂) + noise
where F⃗ may also depend on the computed shape.

Run and tumble in a bounded domain

Spherocylinders in a periodic domain
19 / 29

Tuning the incompressibility: ẋi = F⃗ − τ ∇xi Wα (Leb, µ̂)
Deformabilityc(y, x) = |y − x|αInitial configuration

Mass (inertia)

αo = 1, αb = 2
τo < τbα=1

α=2

α = 1020 / 29

The original idea: incompressible flows à la Brenier
A pioneering idea: The Lagrangian formulation of incompressible fluid flows is asemi-discrete optimal transport problem.
“A combinatorial algorithm for the Euler equations of incompressible flows”[Brenier, Computer Methods in Applied Mechanics and Engineering, 1989]∂t u(t, x) + (u · ∇x )u = −∇x p⇔∇x · u = 0(

d2
dt2 ϕ(t, x) = −∇x p(t, ϕ(t, x))

ϕ#Leb = Leb,

Modern implementation: Using the notion of Laguerre tessellationP
Let ϕ(x) ≃ m(x) = i xi 1Li (x) andconsider a Moreau-Yosida scheme:1d2
m = − 2 ∇m W22 (Leb, m#Leb)2dt
2εwhich is equivalent tod2
1
xi = − 2 ∇xi W22 (Leb, µ̂)
dt2
2ε[Lévy, J. Comput. Phys., 2022], @BrunoLevy01[Gallouët, Mérigot, Found. Comp. Math., 2018]21 / 29

Contact interactions: cell sorting experiments

Cell sorting [Armstrong, 1989]

Partial Engulfment [Armstrong, 1989]

Checkerboard [Togashi et al., 2011]

The Differential Interfacial Tension Hypothesis [Brodland, 2002]

22 / 29

Pressure and surface tension forces in cell doublets
For c(x, y) = |y − x|2 , the Laguerre cells are spherical caps.
Pressure force:Zγbo⃗
Fj→i =⃗n dσ.
|xj − xi | Γij
Triple junction force:z − xi  z − xiF⃗σ = (⃗γbo + ⃗γow + ⃗γbw ) ·.
|z − xi | |z − xi |Equilibrium configuration: when γbw = γow =: γmcos θ =γbo.
2γm−→ Known as the Young-Dupré equation.
[Godard et al., Dev. Cell 55, 2020]23 / 29

Examples
Initial configuration

Checkerboard
γbo < 12 γbb

Separation
γbo > γbw + γow

Sorting and Engulfmentγbo < γbw − γow

24 / 29

Motion driven by random protrusions
Random variations of the cost (protrusions) and a “linear-spring cyto-squeleton”.
 P
(

P
+ p=1 ξp Gθp (θ)φ := − log c, φ(y, x) = − log |y−x|r0dξp = −ξp dt + dBtp , ξp ≥ 0
Similar to a level-set method with potential φ [Yang et al., BMC Sys. Biol. 2(68), 2008][Shafqat-Abbasi et al., 2016]25 / 29

3. Mesoscopic continuum description?

26 / 29Effective diffusion for Active Brownian Particles
Let us consider the particle system xi ∈ Rd defined by:ẋi = c0 ni − τi ∇xi Wα (Leb, µ̂),dni = Pn⊥ ◦ dBti ,for the Wasserstein distance Wα associated to c(x, y) = |x − y|α .
−→ Constant speed motion in the random direction ni ∈ Sd−1 (Brownian motion
on the sphere) coupled with an incompressibility force between the particles.

α = 0.5α=2

α=8−→ Is there a phase transition from infinitely deformable particles (i.e. point
particles) to hard-spheres (transition from a fluid to a glassy state)? [Bi et al., 2015],[Saito, Ishihara, 2023]
27 / 29

Continuum limit?
More generally, is there a mean-field limit for the McKean-Vlasov system
N
1 X
K(Xti , Xtj )dt − τN ∇Xti Wc (Leb, µ̂N )dt + dBti .
dXti =
N j=1
Case c(x, y) = |y − x|2
Assume µ̂N → f (x)dx and τN = O(N ). By Brenier polar factorization theorem,
τN ∇Xti W22 (Leb, µ̂N ) → ∇Φ(x)where the potential Φ satisfies the Monge-Ampère equation:det(I + ∇2 Φ) = f.
Thus the formal mean-field limit:

Z

√
∂t ft (x) = −∇ · ft (x)K(x, y)ft (y)dy + ∇Φ(x)+ 2∆ft (x)Tdwhich is a Fokker-Planck-Monge-Ampère system [Brenier, Loeper, 2004].
Goal: For power costs, find a family of PDEs which interpolates from
non-interacting point particles (α → 0) to hard-spheres (α → +∞).

28 / 29

Thank you for your attention!
ありがとうございます!

29 / 29