###### tags: `optimal transportation`

# Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental Trajectories in Reprogramming

Geoffrey Schiebinger^1,11 Jian Shu^1,2 Marcin Tabaka¹ Brian Cleary^1,3 Vidya Subramanian ¹, Aryeh Solomon¹ Joshua Gould¹ Siyan Liu^1,15 Stacie Lin^1,6 Peter Berube¹ Lia Lee¹ Jenny Chen^1,4 Justin Brumbaugh^5,7,8,9,10 Philippe Rigollet^11,12 Konrad Hochedlinger^7,8,9,13 Rudolf Jaenisch^2,3 Aviv Regev^1,6,13 and Eric S. Lander^1,6,14

¹Broad Institute of MIT and Harvard
²Whitehead Institute for Biomedical Research
³Computational and Systems Biology Program, MIT
⁴Harvard-MIT Division of Health Sciences and Technology
⁵Cancer Center, Massachusetts General Hospital
⁶Department of Biology, MIT
⁷Department of Molecular Biology, Center for Regenerative Medicine and Cancer Center, Massachusetts General Hospital
⁸Department of Stem Cell and Regenerative Biology, Harvard University
⁹Harvard Stem Cell Institute
¹⁰Harvard Medical School
¹¹MIT Center for Statistics
¹²Department of Mathematics, MIT
¹³Howard Hughes Medical Institute
¹⁴Department of Systems Biology Harvard Medical School
¹⁵Biochemistry Program, Wellesley College, Wellesley

[[Paper](https://www.cell.com/cell/fulltext/S0092-8674(19)30039-X)] [[Code](https://github.com/broadinstitute/wot/)]

**Reconstruction of developmental landscapes by optimal-transport analysis of single-cell gene expression sheds light on cellular reprogramming**
[[bioRxiv:10.1101/191056v1](
https://www.biorxiv.org/content/10.1101/191056v1)]


(関連)[TrajectoryNet: A Dynamic Optimal Transport Network for Modeling Cellular Dynamics](https://arxiv.org/abs/2002.04461v1)

https://www.nature.com/articles/s41576-020-0223-2?proof=true

## Pseudotime Cell Trajectories
https://broadinstitute.github.io/2019_scWorkshop/index.html

- [Comparison: Cell Trajectories](https://www.biorxiv.org/content/10.1101/276907v1)
- [Diffusion maps for high-dimensional single-cellanalysis](https://academic.oup.com/bioinformatics/article/31/18/2989/241305)
- [Diffusion pseudotime](https://www.nature.com/articles/nmeth.3971)
- [Slingshot Paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6007078/)
- [Optimal Transport](https://www.cell.com/cell/fulltext/S0092-8674(19)30039-X)
- [RNA Velocity](https://www.nature.com/articles/s41586-018-0414-6)

### Google Slides

<iframe src="https://docs.google.com/presentation/d/e/2PACX-1vQuzaq2kbvEEc3mrUwILcCHuovrKKZWU45EQVEzWISgRVgl3A5KYR1FuY1cS2w0DHG-0wO19zGtvaNj/embed?start=false&loop=false&delayms=3000" frameborder="0" width="760" height="569" allowfullscreen="true" mozallowfullscreen="true" webkitallowfullscreen="true"></iframe>

## ワディントン地形 (Waddington's epigenetic landscape)


## Reconstruction of Probabilistic Trajectories by Optimal Transport


- 分化する細胞集団を時間発展する確率分布 (i.e. 確率過程) $\mathbb{P}_t$ としてモデル化
- $\mathbb{P}_t$ から有限個の経験分布 $\hat{\mathbb{P}}_{t_i}$, $i=1,2,3$ をサンプリングすることで分化プロセスを推測


- 多数のシングルセル発現プロファイル（各点は細胞, 色は採集時間）があれば各時間における分布を近似できる
- 各点は独立にサンプリングされるので各時間間の発現の結合分布: *temporal coupling* $\gamma_{t_1,t_2}$, $\gamma_{t_2,t_3}$ は得られない
- 細胞の状態変化に制約条件が無い場合, temporal coupling を推測することはできない
- 短時間では細胞は変化しないとすれば, 最適輸送を用いて, temporal coupling を推測できる
- 未知である真の結合分布 (カップリング) $\gamma_{t_1,t_2}$ は最適輸送のカップリング $\pi_{t_1,t_2}$ としてみることができ, 経験分布 $\hat{\mathbb{P}}_{t_1}$ と $\hat{\mathbb{P}}_{t_2}$ から計算される $\hat{\pi}_{t_1,t_2}$ で近似できる

### 最適輸送 Optimal transport

- Monge
- Kantorovich

仮定
1. マルコフ過程
2. 分布の時間変化


- $t_{i+1}$ の*子孫分布*は $t_i$ におけるある細胞の集合$C$ を temporal coupling に従う輸送により得られる質量分布
- $t_{i-1}$ の*先祖分布*のは temporal coupling によって時間を巻き戻すことによって得られる質量分布
- 共通祖先は*先祖分布*の混合として表現

### 遺伝子発現空間上の発生過程

> #### 定義1：シングルセル発生トラジェクトリ
> 細胞 $x(0) \in \mathbb{R}^G$を考える. $k(t)\geq0$を$t$における子孫細胞とする ($k(0)=1$). シングルセルの発生過程のトラジェクトリは以下のような連続関数である:
> $$x: [0,T)\to\underbrace{\mathbb{R}^R\times\mathbb{R}^R\times\dots\mathbb{R}^R}_{k(t)\text{ times}}$$.
> これは$x(t)$がそれぞれ$\mathbb{R}^G$のベクトルで表現される$k(t)$組の細胞からなることを意味する.

$$
x(t) = (x_1(t),\dots,x_{k(t)}(t))
$$

$x_1(t),\dots,x_{k(t)}(t)$を$x(0)$の子孫細胞と呼ぶ.

$x_1,\dots,x_n$の細胞集団を確率分布として表現する.

$$
\mathbb{P} = \frac{1}{n}\sum_{i=1}^n \delta_{x_i}
$$

>#### 定義2: 発生過程
>発生過程$\mathbb{P}_t$は遺伝子発現空間において時間発展する確率分布である(確率過程)

>#### 定義3: temporal coupling
>$\mathbb{P}_t$を発生過程として$s<t$の2つのtime pointを考える.
>$X_t\sim\mathbb{P}_t$を$t$におけるランダムな細胞発現プロファイル, $X_s$をその$s$における起源の細胞の発現プロファイルとする.

Temporal coupling $\gamma_{s,t}$を以下のような同時確率分布(の法則??)として定義する

$$
\gamma_{s,t} = \mathcal{L}(X_s,X_t)
$$

同じく

$$
\int_{x \in A}\int_{y \in B}\gamma_{s,t}(x,y)dxdy = \textrm{Pr}\{X_s \in A, X_t \in B\}
$$

for any sets $A,B \subset \mathbb{R}^G$.

Temporal coupling $\gamma_{s,t}$は$\mathbb{P}_s$と$\mathbb{P}_t$のカップリングではない. なぜなら必ずしも$\mathbb{P}_s$や$\mathbb{P}_t$の周辺分布をもたない:

$$
\int \gamma_{s,t} (x,y) dx = \mathbb{P}_t(y)
$$
ではあるが,
$$
\int \gamma_{s,t} (x,y) dy \neq \mathbb{P}_s(x)
$$
である.

生物学的には, 細胞は異なる速度で成長する場合である.初期の段階で増殖する細胞は後で細胞の起源を探す際に過多に表現されてしまう. 発現プロファイルと一定時間後の子孫細胞の平均個数についての関係を説明する相対的な成長率関数を導入する.

>#### 定義4: 相対成長率関数
>temporal couplingに伴う相対成長率関数は以下を満たす関数$g(x)$
>$$
>\int \gamma_{s,t} (x,y) dy = \mathbb{P}_s(x)\frac{g(x)^{t-s}}{\int g(x)^{t-s}d\mathbb{P}_s(x)}
>$$

左辺の積分は$x$から任意の$y$へ向かって出ていく質量である. 右辺の$\mathbb{P}(x)$は発現プロファイル$x$の細胞の豊富さに対応する. $g(x)$は質量の単位時間あたりの指数関数的増加を表す.

### 最適輸送

>#### 例1.
>$X_0\sim\mathcal{N}(0,\sigma^2)$と$X_1\sim\mathcal{N}(\mu,\sigma^2)$を時間$0$と時間$1$におけるある粒子の位置を表す一次元のガウス分布に従う確率変数とする. 短時間で粒子があまり動けないとすると$(X_0,X_1)$という組の同時確率分布を決定するカップリング$\gamma$をどのように推定しよう? 時間$0$から時間$1$まで粒子が動く最小二乗距離を求めるシンプルなヒューリスティクスによって$\hat{\gamma}$を推定する
>
>$$
>\hat{\gamma}\leftarrow\arg\min_{\pi}\mathbb{E}_{\pi}\|X_0-X_1\|^2
>$$
>
>$\mathcal{N}(0,\sigma^2)$と$\mathcal{N}(\mu,\sigma^2)$の周辺分布の全てのカップリング$\pi$について最小化する. 最適な周辺分布が以下のような依存構造をもつ2次元ガウス分布であることが確かめられる.
>$$X_1 = X_0 + \mu$$

この周辺分布をカップリングするヒューリスティクスは**最適輸送**と呼ばれる!

$c(x,y)$を単位質量を$x$から$y$へ輸送するコストとすると, その輸送量は$\pi(x,y)$として, そのような輸送プラン$\pi$に従う質量を輸送する総コストは

$$
\int\int c(x,y)\pi(x,y)dxdy
$$

で与えられる.

この論文ではユークリッド距離の2乗で定義されるコストについてのみ考える. 適切な入力空間上で

$$
c(x,y) = \| x - y \|^2
$$

期待されるコストを最小化する最適輸送プランは周辺制約に従う

$$
\pi(\mathbb{P},\mathbb{Q}) = \textrm{minimize}_{\pi}\int\int c(x,y)\pi(x,y)dxdy \\
\textrm{subject to} \int \pi (x,\cdot)dx = \mathbb{Q}\tag{1} \\
\int \pi (\cdot,y)dy = \mathbb{P}
$$

- $\pi$における線形計画(目的と制約が$\pi$について線型)
- 最適目的値は$\mathbb{P}$と$\mathbb{Q}$間の最適距離を定義する(**Earthmover's distance** or **Wasserstein distance**)
- 他の分布を比較する方法(たとえばKLダイバージェンスやtotal variation)と違い, 最適輸送は基底空間の幾何を考慮する. 例えば, KLダイバージェンスは2つの確率分布の台(サポート)がdisjointであると無限大になってしまうが, 輸送距離は台の離れ具合separationに依存する.

### 最適輸送の原理

1. 従来の最適輸送では制約(1)に質量保存則が入っているが, 前処理として$\mathbb{P}_t$をリスケールすることで増殖を考慮する.
2. 従来の最適輸送で推定されたカップリングは決定論的であり各点は単一の点へと輸送される. しかし, 細胞運命は決定論的でないので, 真のカップリングはエントロピーの程度(degree of entropy)をもつべきである. 少量のエントロピーを入れるのは完全に決定論的な子孫細胞の分布をもつ細胞集団についても意味をもつ. $t_2$において有限個のたくさんの細胞をサンプリングしたとき, $t_1$の任意の細胞の真の子孫細胞は捉えてない. ゆえにエントロピーを輸送写像に導入するのは推定された子孫分布の統計的な不確かさを表すものとして使われる.

$\mathbb{P}_t$: 発生過程
$\gamma_{s,t}$: temporal coupling
$g(x)$: 相対成長関数
$\mathbb{Q}_s$：$\mathbb{P}_s$を相対成長関数でリスケールすることで得られる分布

$$
\mathbb{Q}_s(x) = \mathbb{P}_s(x)\frac{g^{t-s}(x)}{\int g^{t-s}(z)d\mathbb{P}_s(z)}
$$

$\pi_{s,t}(\epsilon)$ ：以下の最適化問題の解として定義される$\mathbb{Q}_s$と$\mathbb{Q}_t$のエントロピー正則化した最適輸送カップリング

$$
\pi_{s,t}(\epsilon) = \textrm{minimize}_{\pi}\int\int c(x,y)\pi(x,y)dxdy - \epsilon\int\int\pi(x,y)\log\pi(x,y)dxdy \\
\textrm{subject to} \int \pi (x,\cdot) dx = \mathbb{Q}_s \tag{2} \\
\int \pi (\cdot,y)dy = \mathbb{P}_t
$$

発生過程における最適輸送原理：
$$
s \approx t \Longrightarrow \pi_{s,t}(\epsilon)\approx \gamma_{s,t}
$$

### Temporal couplingをデータから推定する

>#### 定義5 developmental time series
>発生時系列は発生過程$\mathbb{R}^G$上の$\mathbb{P}_t$からのサンプル列である. これは$t_1,\dots,t_T \in \mathbb{R}$の時間で採集された集合の列$S_1,\dots,S_T \in \mathbb{R}^G$である. 各$S_i$は$\mathbb{P}_{t_i}$から独立にサンプリングされる$\mathbb{R}^G$上の発現プロファイルの集合である.

この入力データから経験的な発生過程をつくる. 各$t_i$において $x \in S_i$にサポートされた経験分布をつくる.

>#### 定義6 Empirical developmental process
>経験的発生過程 $\hat{\mathbb{P}}_t$ は発生時系列$S_1,\dots,S_T$からつくられる時間発展する確率分布:
>$$
>\hat{\mathbb{P}}_{t_i} = \frac{1}{|S_i|}\sum_{x \in S_i}\delta_x
>\tag{3}
>$$
>経験的発生過程は $t \notin \{t_1,\dots,t_T\}$ においては定義されていない

$t_1$から$t_2$へのカップリングを推定するため, 先ず成長率関数$g(x)$の最初の推定をつくる. 推定された関数$\hat{g}(x)$は遺伝子発現空間での誕生と死亡の過程の期待値を細胞増殖とアポトーシスに関与する遺伝子発現レベルに関する誕生率$\beta(x)$と死亡率$\delta(x)$で定義する

**unbalaced transport** (Chizat et al. 2018)
を使ってデータから誕生率と死亡率の推定を改善する

リスケールされた経験分布をつくり:
$$
\hat{\mathbb{Q}}_{t_1}(s) = \hat{\mathbb{Q}}_{t_1}(x)\frac{\hat{g}(x)^{t_1-t_2}}{\int\hat{g}(z)^{t_1-t_2}d\hat{\mathbb{P}}_{t_i}(z)}
$$

$\hat{\mathbb{Q}}_{t_1}$と$\hat{\mathbb{P}}_{t_2}$間の最適輸送写像$\hat{\pi}_{t_1,t_2}$を計算する


#### 隣接するtime pointのカップリングを推定する

- エントロピーがゼロなら($\epsilon=0$)線型計画
- [Cuturi, 2013] がエントロピー正則化した輸送問題について大規模な場合でも高速のアルゴリズムを提案している
- 強凸
- gradient ascent on dual can be reealised by successive diagonal matrix scalings : Sinkhorn iterations

等式制約が輸送プランの周辺分布のバウンドとして緩和される(KL-ダイバージェンス, total variation or general f-divergenceに関して)
これはとても魅力的である

1. $\hat{g}(x)$を誤推定しているかも. Unbalanced transportは入力された成長率を輸送コストを減らすように調整する. 成長率は自動的に学習される.
2. 成長率が完全に一様分布だとしてもランダムサンプリングは成長っぽいなにかを導入してくれる. 例えば全体の5%のレアなサブ集団があったとして, ある時点で, これらの細胞をランダムに少なくサンプリングしたとして全体の4%になってしまったとする. そして次の時点で6%にしてしまうとする. そうすると集団が50%増加したようにみえてしまう. Unbalanced transportはこのような見かけ上の成長を自動的に調整していくれる.

$$
\hat{\pi}_{t_i,t_{i+1}} = \arg\min_{\pi} \sum_{x\in S_i}\sum_{y\in S_{i+1}}c(x,y)
\pi(x,y)-\epsilon\int\int\pi(x,y)\log\pi(x,y)dxdy\\
+\lambda_1\textrm{KL}\left[\sum_{x\in S_i}\pi(x,y)\mid\mid d\hat{\mathbb{P}}_{t_{i+1}}(y)\right]+\lambda_2\textrm{KL}\left[\sum_{y\in S_{i+1}}\pi(x,y)\mid\mid d\hat{\mathbb{Q}}_{t_i}(x)\right] \tag{4}
$$

where $\epsilon$, $\lambda_1$, $\lambda_2$ are regularisation parameters.

行列変数$\pi \in \mathbb{R}^{N_i \times N_{i+1}}$における凸最適化問題である. ここで$N_i=|S_i|$は$t_i$でプロファイルされた細胞数である. 普通のノートPCで$N_i\approx 5000$を計算するのに5秒間くらい.

式(3)によって定められる経験分布の離散集合$S_i$上の確率密度は$d\hat{\mathbb{P}}_{t_i}(x)=\frac{1}{N_i}$であるが non-uniform な経験分布を使うこともできる (細胞の質についての情報を使うとか)

発現プロファイル列$S_1,\dots,S_T$が与えられた時, 連続するtime pointsの組$S_i$,$S_i+1$の最適化問題(4)を解く. $(t_i,t_{i+1})$について$\hat{\pi}_{t_i,t_{i+1}}$を得る. 十分なデータを得られれば, $\pi_{t_i,t_{i+1}}$に関するよい推定が得られる

$$
\lim_{N_i,N_{i+1}\to\infty}\hat{\pi}_{t_i,t_{i+1}} = \pi_{t_i,t_{i+1}}
$$

最適輸送原理:
$$
\pi_{t_i,t_{i+1}} \approx \gamma_{t_i,t_{i+1}}
$$

### Estimating long-range couplings
マルコフ性を仮定すれば長時間のインターバルにおいても推定可能である

>#### 定義7 Markov developmental process
>マルコフ発生過程$\mathbb{P}_t$は$\mathbb{R}^G$上で時間発展する確率分布で２つの時間の組間のカップリングで完全に定まる
>$s<t<\tau$において長時間カップリング$\gamma_{s,\tau}$は短時間のカップリングの合成と等しい
>$$
>\gamma_{t,\tau} \circ \gamma_{s,t} = \gamma_{s,\tau}
>$$

最適輸送写像$\hat{\pi}_{s,t}$はこの合成則を必ずしも満たさない. $s$から$t$と$t$から$\tau$への最適輸送は$s$から直接$\tau$への最適輸送と同じでない.

任意の$t_i$と$t_{i+k}$の組でカップリング$\hat{\gamma}_{t_i,t_{i+k}}$を合成則により推定する

$$
\hat{\gamma}_{t_i,t_{i+k}} = \hat{\pi}_{t_i,t_{i+1}} \circ \hat{\pi}_{t_{i+1},t_{i+2}} \circ \dots \hat{\pi}_{t_{i+k-1},t_{i+k}}
$$

- マルコフ過程としてみるのは議論の余地あり
- 遺伝子発現の歴史はクロマチン変異に影響を与えるがこれは観測された発現プロファイルには反映されない

### 輸送写像を解釈する

#### 先祖, 子孫, 軌跡

輸送写像の*押し出し* *pushing*に依存する

積空間$\mathbb{R}^G\times\mathbb{R}^G$上の確率分布として, カップリング$\gamma$は任意の組$A,B\subset\mathbb{R}^G$に$\gamma(A,B)$をアサインする

$$
\gamma(A,B) = \int_{x \in A}\int_{y \in B}\gamma(x,y)dxdy
$$

$\gamma(A,B)$は$A$から来て$B$へ行く質量.

特定の行き先を定めなければ, $\gamma(A,\cdot)$という量は$A$から出てくる全質量の分布になっている.

これを輸送プラン$\gamma$の*pushing* $A$と呼ぶ. もっと一般化すれば, 分布$\mu$をpush foward可能である.

$$
\mu \mapsto \int \gamma(x,\cdot)d\mu(x)
$$

逆の操作を*pulling 引き戻し* $B$と呼ぶ. これにより得られる分布$\gamma(\cdot,B)$は$B$にたどりつく質量を表す. 同じようにpull分布も可能

$$
\mu \mapsto \int \gamma(\cdot,y)d\mu(y)
$$

これを*back-propagating* the distribution $\mu$. (pushing $\mu$は*forward-propagation*)

>#### 定義8 descendants in Markov developmental process
>細胞集団$C\in \mathbb{R}^G$を考える
>$t_1$において細胞集団の一部として生きていてマルコフ発生過程$\mathbb{P}_t$に従って変化していく
>$\gamma_{t_1,t_2}$を$t_1$から$t_2$へのカップリングとする
>$t_2$における$C$の子孫とは$C$を$\gamma$において*pushing*することで得られる

>#### 定義9 ancestors in Markov developmental process
>細胞集団$C\in \mathbb{R}^G$を考える
>$t_2$において細胞集団の一部として生きていてマルコフ発生過程$\mathbb{P}_t$に従って変化していく
>$\pi$を$t_2$から$t_1$への輸送写像とする
>$t_1$における$C$の先祖とは$C$を$\gamma$において*pulling*することで得られる
>

#### トラジェクトリ
集合$C$への*ancestor trajectory* を先祖分布の列として定義する
*descendant trajectory* は集合$C$からの依存分布の列として定義する

### エントロピー正則化のパラメータを解釈する

$N$個の区別できない粒子が拡散係数$\epsilon$でブラウン運動をしている. $N$個の粒子の位置を時間$0$と時間$1$において観測するとする. 粒子は区別できないので時間$0$のどの粒子が時間$1$のどの粒子と対応付くかわからない. もし$N=1$ならば, これは問題ではないが, 始点と終点を結ぶことを*Brownian bridge*と呼ぶ.

$N>1$において, 始点と終点を結ぶ取りうる経路についての分布は

1. どの粒子がどこへ行くのかを定める粒子のカップリング（観測からは一意に定まらない）
2. マッチングが与えられたときに各マッチングされた組の経路上の分布がBrownian bridge

- カップリングは時刻$0$と時刻$1$の点におけるランダムな組み合わせである. このランダムな組み合わせの分布はブラウン運動の分散(または拡散係数)に従う
- 拡散係数が大きければ, 大きく位置をスワップしがち
- カップリングの期待値(平均値)はエントロピー最大化による最適輸送によって得られる
- これらのアイデアはSchrodinger 1932に遡る
- Cuturi 2013, Leonard 2014

>#### 定理
>- ブラウン運動の拡散係数がエントロピー正則化パラメータと等しいとき, エントロピー正則化最適化輸送はブラウン運動のカップリングについての分布の期待値を与える

### Gradient flow and Waddington'S landscape
最適輸送が遺伝子発現空間上(capturing cell-autonomous processes)あるいは分布空間(capturing cell-nonautonomous processes)の勾配流として解釈できる

最適輸送と勾配流については
Ambrosio et al. 2005
Santambrogio, 2015
を参照

Waddington地形は発現プロファイル$x$の細胞にエネルギーポテンシャル$\Phi(x)$をアサインするポテンシャル関数$\Phi$を定義する

この細胞は$\Phi$の勾配に沿って転がり, 以下の微分方程式を満たす軌跡$x(t)$を説明する

$$
\frac{dx}{dt} = -\nabla\Phi(x)
\tag{5}
$$

式(5)は各細胞の軌跡を支配していて細胞集団の分布の中で流れを作りだす

$$
\frac{d\mathbb{P}_t}{dt} = \textrm{div}[\nabla\Phi(x)\mathbb{P}_t]
\tag{6}
$$

式(6)は各微小体積空間の質量の変化を定め（左辺）質量のfluxのinとoutに等しい(given by the divergence on the right hand side)

最適輸送はこの種のポテンシャルに駆動されるダイナミクスを捉えることができる

式(5)により定められる真のカップリングは短時間スケールにおいての最適輸送カップリングに近い.

#### 定理2 Benamou and Brenier 2001
輸送問題の最適目的値(1)は以下の最適化問題の最適目的値と等しい
$$
\min_{\rho,v} \int_0^1\int_{\mathbb{R}^G}\|v(t,x)\|^2\rho(t,x)dtdx \\
\textrm{subject to } \rho(0,\cdot) = \mathbb{P}, \rho(1,\cdot) = \mathbb{Q} \\
\nabla\cdot(\rho v) = \frac{\partial \rho}{\partial t}
\tag{7}
$$

- $v$はベクトル場で分布$\rho$と$\mathbb{P}$から$\mathbb{Q}$へadvect移流し, 最小化される目的値は流れの力学的エネルギー(mass $\times$ squared veolocity)である.
- 我々の設定ではこの２つの分布は$\mathbb{P}_s$と$\mathbb{P}_t$という発生過程の2つの時刻でのスナップショットである
- 輸送写像$\pi_{s,t}$は連続時間の流れの最小作用のpoint-to-point summary
- ベクトル場がポテンシャル$\Phi$の勾配であるという特別な場合(Waddington地形)に定理はカップリング(5)が最適輸送コストを達成する
- 言い換えると, 最適輸送はポテンシャル駆動のダイナミクスを記述
- このベクトル場が時間変化し, 細胞の分布全体によるならば, 最適輸送が一般化された発生過程（細胞間相互作用等を含む）を記述

時間発展(6)がWasserstein勾配流の特別な場合であり, 以下の線型なエネルギー汎関数を最小化することを示そう

$$
E(\mathbb{P}) = \int\Phi(x)d\mathbb{P}(x)
$$

細胞間相互作用を捉える非線型な勾配流を記述する

勾配降下

$$
x_{k+1} = -\eta\nabla E (x_k) + x_k
$$

これは*proximal procedure*として書き直せる. これは$x_k$のproximityにおいてすべての$x$に対して$E$を最小化する

$$
x_{k+1} = \arg\max_x E(x)+\frac{1}{2\eta}\|x-x_k\|^2
\tag{8}
$$

分布空間で似たproximal procedureをとり, ユークリッドノルム$\|\cdot\|^2$をWasserstein距離で置き換えると,

$$
\mathbb{P}_{k+1} = \arg\min_{\rho} E(\rho)+\frac{1}{2\eta}W^2_2(\rho,\mathbb{P}_k)
$$

これは $\mathbb{P}_0,\mathbb{P}_1,\dots,\mathbb{P}_k$の反復された列をつくる. この勾配流はステップ幅$\eta$をゼロにもっていった極限で, Jordan et al. 1998によれば, 線型のエネルギー汎関数について

$$
E(\mathbb{P}) = \int\Phi(x)d\mathbb{P}(x)
$$

勾配流の極限は(6)式の解へと収束する

細胞間相互作用は相互作用エネルギー

$$
E(\mathbb{P}) = \int\int I(x,y)d\mathbb{P}(x)d\mathbb{P}(y)
$$

Santambrogio, 2015 に詳細がある

### 遺伝子制御のモデルを学習する
- ベクトル場をモデルとして解釈する
- 現時刻の転写因子の発現の関数として以後の時刻の遺伝子発現を予測できる
- ベクトル場があまり変化しないと仮定する

## Waddington-OT

### 最適輸送写像の計算
- PCAした空間でのユークリーッド二乗距離がコスト関数
- 時刻間内での中央値で正規化する
- 正規化のパラメータ
- 相対成長率の推定
### 先祖/子孫/軌跡の推定

$$
p(x) =
\cases{
1 & $x \in C$\\
0 & otherwise
}
$$

$C$の子孫の分布は$\hat{\pi}_{t_1,t_2}p$
先祖の分布は$p^{\textsf{T}}\hat{\pi}_{t_0,t_1}$

### Learning gene regulatory models

### Geodesic interpolation for validation

## A Dense scRNA-Seq Time Course of iPS Reprogramming


Reprogramming of secondary ($2^\circ$) MEFs from E13.5 embryos. Each dot represents a collection time point.


### Cell signatures
### Ancestor descendant flow

- Stromal/MET fate decision occurs at the early stage

## The Model Is Predictive and Robust
- 細胞系譜追跡からは細胞集団の起源の転写プロファイルについてはわからない
- $t_1<t_2<t_3$において, 提案法により$t_1$から$t_3$へのトラジェクトリを補完していくことで$t_2$の分布を推定する
- $t_2$の観測されたものと比較

## In Initial Stages of Reprogramming, Cells Progress toward Stroma


#### Stromal gene activity

#### Stromal
#### MET

## iPSCs Emerge through a Tight Bottleneck from Cells in the MET Region




## Development of Extra-Embryonic-like Cells during Reprogramming


## RNA Expression Reveals Genomic Aberrations in Trophoblast-like and Stromal Cells


## Neural-like Cells Also Emerge from the MET Region during Reprogramming in Serum


## The Developmental Landscape Highlights Potential Paracrine Signals


## Experimental Validation Confirms that TranscriptionFactor Obox6 and Cytokine GDF9 Enhance Reprogramming

- Obox6, Tdgf1 というreprogramming関連因子を特定
- Oct4-EGFP overexpreessionにより検証

### Obox6 ↑ regulation in iPSC-fated cells
### Oct4-EGFPイメージング w/ Obox6 supply
### GDF9

## Discussion


>多時間single-cell RNA-seqでの子孫関係推定を最適輸送問題として定式化することで分化軌道を推定する論文を読んだ。MEFのiPS細胞リプログラミングに適用することで寄与する新たなパラクリンシグナルや遺伝子も発見している。Waddington地形への執着が生物学者らしくて好き

>MEFからiPS細胞への誘導(reprogramming)における>315000件のsingle cell RNA-seqのデータから発生に関わる多様な経路を明らかにした