USING TOPOLOGY TO MEASURE DYNAM- ICS OF BIOLOGICAL AGGREGATIONS - - PowerPoint PPT Presentation

USING TOPOLOGY TO MEASURE DYNAM- ICS OF BIOLOGICAL AGGREGATIONS

August 7, 2019

Lori Ziegelmeier, Macalester College

ICERM Workshop: Applied Mathematical Modeling with Topological Techniques

INTRODUCTION

BIOLOGICAL AGGREGATIONS

In many natural systems, particles, organisms, or agents interact locally according to rules that produce aggregate behavior.

3

CLASSIC WAY TO ANALYZE BIOLOGICAL AGGREGATIONS

Alignment Order Parameter: ϕ(t)

1 Nv0

• N
• i1
• vi(t)
• 7

4

CLASSIC WAY TO ANALYZE BIOLOGICAL AGGREGATIONS

Alignment Order Parameter: ϕ(t)

1 Nv0

• N
• i1
• vi(t)
• 7

4

CLASSIC WAY TO ANALYZE BIOLOGICAL AGGREGATIONS

Alignment Order Parameter: ϕ(t)

1 Nv0

• N
• i1
• vi(t)
• 7

4

TOPOLOGICAL DATA ANALYSIS

• 1. Envision data as a point cloud
• e.g. position-velocity for one snapshot in time
• 2. Create connections between proximate points
• build simplicial complex
• 3. Determine topological structure of complex
• compute homology (measure # holes)
• 4. Vary proximity parameter to assess different scales
• calculate persistent homology

5

COMPUTE PERSISTENT HOMOLOGY

b0 ε

6

COMPUTE PERSISTENT HOMOLOGY

b0 7 6 3 2 1

6

COMPUTE PERSISTENT HOMOLOGY

6

COMPUTE PERSISTENT HOMOLOGY

6

TOPOLOGICAL DATA ANALYSIS

• 1. Envision data as a point cloud
• e.g. position-velocity for one snapshot in time
• 2. Create connections between proximate points
• build simplicial complex
• 3. Determine topological structure of complex
• compute homology (measure # holes)
• 4. Vary proximity parameter to assess different scales
• calculate persistent homology

7

TOPOLOGICAL DATA ANALYSIS

• 1. Envision data as a point cloud
• e.g. position-velocity for one snapshot in time
• 2. Create connections between proximate points
• build simplicial complex
• 3. Determine topological structure of complex
• compute homology (measure # holes)
• 4. Vary proximity parameter to assess different scales
• calculate persistent homology
• 5. Measure topology as time evolves.
• Crocker plots

7

EVOLVE IN TIME

Compute the kth Betti number bk(ε, t),

8

EVOLVE IN TIME

Compute the kth Betti number bk(ε, t), CROCKER plot Contour Realization Of Computed K-dimensional hole Evolution in the Rips complex (CROCKER)

(Topaz, Z., Halverson 2015) Topological Data Analysis of Biological Aggregation Models 8

CROCKER AS EXPLORATORY TOOL

D’ORSOGNA MODEL Dynamical system describing motion of interacting point

particles in an unbounded plane in continuous time.

Model written as:

• xi
• vi

i 1, . . . , N m

• vi
• α − β|

vi|2 vi − ∇iUi Ui

• ji

Cre−|

xi− xj |/ℓr − Cae−| xi− xj |/ℓa

After nondimensionalization, we have 4 parameters: α, β, C Cr

Ca

, ℓ ℓr ℓa

(D’Orsogna, et al 2006) Self-Propelled Particles with Soft-Core Interactions: Patterns, Stability, and Collapse 10

SNAPSHOTS IN TIME α 1.5, β 0.5, C 2, ℓ 0.25 with N 500 particles

11

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• 12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Are the particles moving in the same

direction?

12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• 12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• Are the particles rotating (in the same

direction)?

12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• Absolute Angular Momentum:

Mabs(t)

• N

i1 |ri(t)×vi(t)|

N

i1 |ri(t)||vi(t)|

• 12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• Absolute Angular Momentum:

Mabs(t)

• N

i1 |ri(t)×vi(t)|

N

i1 |ri(t)||vi(t)|

• Are the particles rotating?

12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• Absolute Angular Momentum:

Mabs(t)

• N

i1 |ri(t)×vi(t)|

N

i1 |ri(t)||vi(t)|

• Average Nearest Neighbor

Distance: Dnn(t) 1

N N

• i1

min

1≤j≤N |xi(t) − xj(t)|

12

ORDER PARAMETERS TO SUMMARIZE COLLECTIVE BEHAVIOR

Polarization: P(t)

• N

i1 vi(t)

N

i1 |vi(t)|

• Angular Momentum:

Mang(t)

• N

i1 ri(t)×vi(t)

N

i1 |ri(t)||vi(t)|

• Absolute Angular Momentum:

Mabs(t)

• N

i1 |ri(t)×vi(t)|

N

i1 |ri(t)||vi(t)|

• Average Nearest Neighbor

Distance: Dnn(t) 1

N N

• i1

min

1≤j≤N |xi(t) − xj(t)|

How close are the particles?

12

D’ORSOGNA SIMULATION ANALYSIS

(A)

0.00 0.25 0.50 0.75 1.00 10 20 30 40 50

Simulation Time Order Parameters

Order Parameter P M Mabs

b0 ≥ 5 b0 = 1 (B)

0.0 0.5 1.0 1.5 10 20 30 40 50

Simulation Time t Proximity Parameter ε

level 1 2 3 4 5

b1 = 0 b1 ≥ 5 b1 = 1 b1 = 2 (C)

0.0 0.5 1.0 1.5 10 20 30 40 50

Simulation Time t Proximity Parameter ε

level 1 2 3 4 5 (A) t = 5

−0.6 0.6 −0.6 0.6

(B) t = 23

−0.6 0.6 −0.6 0.6

(C) t = 34

−0.6 0.6 −0.6 0.6

13

D’ORSOGNA SIMULATION ANALYSIS

(A)

0.00 0.25 0.50 0.75 1.00 10 20 30 40 50

Simulation Time Order Parameters

Order Parameter P M Mabs

b0 ≥ 5 b0 = 1 (B)

0.0 0.5 1.0 1.5 10 20 30 40 50

Simulation Time t Proximity Parameter ε

level 1 2 3 4 5

b1 = 0 b1 ≥ 5 b1 = 1 b1 = 2 (C)

0.0 0.5 1.0 1.5 10 20 30 40 50

Simulation Time t Proximity Parameter ε

level 1 2 3 4 5 (A) t = 5

−0.6 0.6 −0.6 0.6

(B) t = 23

−0.6 0.6 −0.6 0.6

(C) t = 34

−0.6 0.6 −0.6 0.6

13

PARAMETER IDENTIFICATION

QUESTION OF INTEREST

Forward Problem:

15

QUESTION OF INTEREST

Forward Problem: Inverse Problem: Research Question What is an informative way to summarize collective behavior for machine learning techniques?

15

PHENOTYPES IN THE D’ORSOGNA MODEL

Different patterns emerge by altering parameters C and ℓ

16

PIPELINE

time scale

Simulate Swarm Behavior

25 parameter sets x 100 realizations = 2500 simulations

time value Topology based Description Order Parameter based Description

ε/2

Predict Parameters & Patterns

Machine Learning

Construct Problem independent Problem dependent

Topology

Design

Order Parameters

Construct Feature Vectors

17

PHENOTYPES AND FEATURE VECTORS

Single Mill: Double Ring: Escape:

18

SUPPORT VECTOR MACHINE Supervised machine learning algorithm 5-fold cross validation, with each training simulation labeled

with (C, ℓ), 20% of data used for each test

Accuracy out-of-sample simulations with (C,ℓ) correct

• ut-of-sample simulations

19

CLASSIFICATION RESULTS

Summary Feature Dimension Accuracy P(t) 87 57.7% Order Mang(t) 87 34.4% Parameters Mabs(t) 87 68.0% DNN(t) 87 91.1% All 4*87 89.2% TDA b0 200*86 99.6% (time-delayed b1 200*86 99.3% position) b0 & b1 2*200*86 99.1%

20

CLASSIFICATION RESULTS

Summary Feature Dimension Accuracy P(t) 87 57.7% Order Mang(t) 87 34.4% Parameters Mabs(t) 87 68.0% DNN(t) 87 91.1% All 4*87 89.2% TDA b0 200*86 99.6% (time-delayed b1 200*86 99.3% position) b0 & b1 2*200*86 99.1% Summary Feature Dimension Accuracy P(t) (PCA) 3 46.7% Mang(t) (PCA) 3 30.0% Order Mabs(t) (PCA) 3 58.8% Parameters DNN(t) (PCA) 3 81.5% All (PCA) 3 68.6% b0 (PCA) 87 99.7% TDA b1 (PCA) 87 99.9% (time-delayed b0 & b1 (PCA) 87 99.7% position) b0 (PCA) 3 89.7% b1 (PCA) 3 82.8% b0 & b1 (PCA) 3 89.6%

20

CLASSIFICATION RESULTS

(Bhaskar, Manhart, Milzman, Nardini, Storey, Topaz, Z. forthcoming) Analyzing Collective Motion with Machine Learning and Topology 21

ASSESSING MODEL VALIDITY

MOVEMENT OF PEA APHIDS

23

MOVEMENT OF PEA APHIDS

23

MOVEMENT OF PEA APHIDS

23

PEA APHID MODELS Classifies each aphid according to two motion states, moving

• r stationary.

Transitions between these states are probabilistic, depending

• n distance d to each aphid’s nearest neighbor.

24

PEA APHID MODELS Classifies each aphid according to two motion states, moving

• r stationary.

Transitions between these states are probabilistic, depending

• n distance d to each aphid’s nearest neighbor.

Model written as:

• Probability of transition state

PMS(d) P∞

MS +

P0

MS − P∞ MS

• e−d/dMS

PSM(d) P0

SMe−d/dSM + P∞ SM

d d + ∆SM

• Step length ℓ

ℓ(d) ℓ∞ + ℓ0 − ℓ∞ e−d/dℓ

• Turning angle θ drawn from a wrapped Cauchy distribution

centered at zero, with parameter ρ(d) controlling the spread of the distribution.

24

PEA APHID MODELS Classifies each aphid according to two motion states, moving

• r stationary.

Transitions between these states are probabilistic, depending

• n distance d to each aphid’s nearest neighbor.

Two models: Interactive and Noninteractive (Control)

24

ASSESSING MODEL VALIDITY

Goal: Use order parameters and crocker plots to compare each of the interactive and control models to experimental data

25

ORDER PARAMETERS OF EXPERIMENT AND SIMULATIONS

Polarization

5000 0.5 1

Experiment

5000 0.5 1

Interactive

5000 0.5 1

Control

Angular Momentum

5000 0.5 1

Experiment

5000 0.5 1

Interactive

5000 0.5 1

Control

Absolute Angular Momentum

5000 0.5 1

Experiment

5000 0.5 1

Interactive

5000 0.5 1

Control

26

ORDER PARAMETERS OF EXPERIMENT AND SIMULATIONS

Average Nearest Neighbor Distance

5000 0.01 0.02

Experiment

5000 0.01 0.02

Interactive

5000 0.01 0.02

Control

Percent of Aphids Moving

5000 0.5 1

Experiment

5000 0.5 1

Interactive

5000 0.5 1

Control

27

B0(POS) CROCKER PLOTS OF EXPERIMENT AND SIMULATIONS

28

EXPERIMENT VS MODEL COMPARISON

Histogram of crocker difference between control/experiment and interactive/experiment

29

SUMMARY OF STATISTICAL TESTS COMPARING MODELS

30

SUMMARY OF STATISTICAL TESTS COMPARING MODELS

30

SUMMARY OF STATISTICAL TESTS COMPARING MODELS

(Ulmer, Z., Topaz 2018) Assessing Biological Models Using Topological Data Analysis 30

PERSISTENT CROCKER PLOTS

BENEFITS AND DRAWBACKS OF CROCKERS

Benefits:

Can be used to summarize time-varying metric spaces Displays topological information at all times simultaneously Can be vectorized and applied to statistical and machine learning

tasks

32

BENEFITS AND DRAWBACKS OF CROCKERS

Benefits:

Can be used to summarize time-varying metric spaces Displays topological information at all times simultaneously Can be vectorized and applied to statistical and machine learning

tasks Drawbacks:

Crockers are not stable. Perturbing the dynamic metric space only

slightly could produce changes of unbounded size.

32

STABILITY OF PERSISTENCE DIAGRAMS

Stability Ideal if tools for data analysis are stable with respect to small per- turbations of the inputs.

33

STABILITY OF PERSISTENCE DIAGRAMS

Stability Ideal if tools for data analysis are stable with respect to small per- turbations of the inputs.

33

DISTANCE BETWEEN PERSISTENCE DIAGRAMS

death birth

Definition

The bottleneck distance between two PDs B and B′ is given by d∞(B, B′)

inf

γ:B→B′sup u∈B

||u − γ(u)||∞,

ranging over all bijections between B and B′.

34

DISTANCE BETWEEN METRIC SPACES

Definition

If X and Y are two subsets of a metric space Z, then the Hausdorff distance between X and Y is dZ

H(X, Y) max{sup x∈X

inf

y∈Y d(x, y), sup y∈Y

inf

x∈X d(x, y)}.

The Gromov–Hausdorff distance between metric spaces X and Y is dGH(X, Y) inf

Z,f,g dZ H(f(X), g(Y)),

where the infimum is taken over all metric spaces Z and isometric embeddings f : X → Z and g: Y → Z.

35

PERSISTENCE DIAGRAMS ARE STABLE

Stability of PDs (Chazal, de Silva, Oudot 2013) PDs are stable (Lipschitz) with re- spect to the bottleneck metric d∞(PD(VR(X)), PD(VR(Y))) ≤ 2dGH(X, Y) for metric spaces X and Y.

36

RANK INVARIANT

Definition

For a persistence module V and ε < ε′, the rank of the map V(ε) → V(ε′), is the number of intervals in the persistence barcode that contain the interval [ε, ε′]. The collection of all natural numbers rank(V(ε) → V(ε′)) for all ε < ε′ is called the rank invariant.

1 2 3 4 5 6 7 8 9 10 11

Figure: rank(V(4) → V(8)) 2

The rank invariant is equivalent to a persistence barcode i.e. can

• btain one from the other (Carlsson and Zomorodian, 2009).

37

α-SMOOTHED CROCKER PLOT

Let X be a time-varying metric space, with Xt the space at time t. An α-smoothed crocker plot, for α ≥ 0, at time t and scale ε is equal

to rank(Hk(VR(Xt; ε − α)) → Hk(VR(Xt; ε + α))).

A standard crocker plot is a 0-smoothed crocker plot.

1 2 2 1 2 2 2 1

Persistence Diagram Persistence Intervals Vertical Slice of Crocker

1 1 1 1 1 1

Persistence Diagram Persistence Intervals Vertical Slice of

• persistence line

38

CROCKER VIDEO

Definition (Crocker Video) A crocker video is a sequence of α-smoothed crocker plots as α increases continuously from 0. It is an integer-valued function f : R3 → N f(t, ε, α) rank(Hk(VR(Xt; ε − α)) → Hk(VR(Xt; ε + α))).

39

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

CROCKER VIDEOS

A crocker video is a non-increasing function, f(t, ε, α) ≤ f(t, ε, α′) for α ≥ α′.

40

FEATURES OF CROCKER VIDEOS

The crocker video is equivalent to the vineyard (Cohen-Steiner,

Edelsbrunner, Morozov 2006) i.e. can obtain one from the other.

The crocker video inherits some nice interpretability properties of the

crocker plot.

All times are represented in each frame α, whereas in a vineyard

• ne only sees information about a single time.

Crocker videos satisfy a stability property similar to the stability of

vineyards.

41

CROCKER VIDEOS ARE STABLE

Theorem (Stability theorem for crocker videos) If X and Y are totally-bounded time-varying metric spaces, and if d∞

GH(X, Y) ≤ δ/2, then the crocker videos for X and Y are close in

the sense that for all t, ε, and α, we have

fX(t, ε, α + δ) ≤ fY(t, ε, α), and fY(t, ε, α + δ) ≤ fX(t, ε, α).

Here, by an abuse of notation, we let fX and fY denote fPH(VR(X)) and fPH(VR(Y)).

(Adams, Topaz, Xian, Z. forthcoming) Capturing Dynamics of Time-Varying Systems with Topology 42

CONCLUSION

CONCLUSION Crocker plots (0-persistence crockers) are useful summaries of

topological properties of time-varying systems.

• Can be combined with machine learning and statistical tools
• Provide means of identifying global behaviors, parameter

identification, and model selection

Crocker videos are stable.

• If two time-varying metric spaces are “close", then their

corresponding crocker videos will also be “close"

44

QUESTIONS?

Lori Ziegelmeier lziegel1@macalester.edu Department of Mathematics, Statistics, and Computer Science NSF:CDS&E-MSS-1854703 Join WinCompTop: email WinCompTop+subscribe@googlegroups.com

45