density of $\mu_{0}$. The randomness of the above SPDE only comes from the
innovation process
$\mathrm{d}\nu_{t}=\mathrm{d}y_{t}-\mu_{t}\left(H\right)\mathrm{d}t$ as a
finite-dimensional white noise in time.

At last, if the functions on the right hand sides of (A1) satisfy the
Ornstein-Uhlenbeck processes with matrix coefficients, that is,
$F\left(u\right)=Fu+f_{t}$ and $H\left(u\right)=Hu+h_{t}$. With Gaussian
initial condition, $u_{0}\sim\mathcal{N}\left(\hat{u}_{0},C_{0}\right)$, the
conditional distribution $\mu_{t}=\mathcal{N}\left(\hat{u}_{t},C_{t}\right)$
given $\mathcal{G}_{t}$ in (A2) becomes a multivariate normal distribution,
where $\hat{u}_{t}=\mathbb{E}\left[u_{t}\mid\mathcal{G}_{t}\right]$ and
$C_{t}=\mathbb{E}\left[\left(u_{t}-\hat{u}_{t}\right)\left(u_{t}-\hat{u}_{t}\right)^{\intercal}\mid\mathcal{G}_{t}\right]$.
The filtering equations for $\hat{u}_{t}\in\mathbb{R}^{d}$ and
$C_{t}\in\mathbb{R}^{d\times d}$ are given by the Kalman-Bucy filter [28]

$\displaystyle\mathrm{d}\hat{u}_{t}$
$\displaystyle=\left(F\hat{u}_{t}+f_{t}\right)\mathrm{d}t+K_{t}\left[\mathrm{d}y_{t}-\left(H\hat{u}_{t}+h_{t}\right)\mathrm{d}t\right],$
(A4a) $\displaystyle\dot{C}_{t}$
$\displaystyle=FC_{t}+C_{t}F^{\intercal}-K_{t}K_{t}^{\intercal}+\Sigma\Sigma^{\intercal}.$
(A4b)

with the Kalman gain matrix $K_{t}=C_{t}H^{\intercal}$. Above, (A4a) is an SDE
coupled with the deterministic Riccati equation (A4b).

###  A.2 Infinite dimensional filtering in Hilbert space

It is shown that the linear Kalman-Bucy filter can be generalized to linear
stochastic equations on a Hilbert space [10, 20]. Let $H$ be a Hilbert space.
Denote $L^{2}\left(\Omega,\mathcal{G},\mathbb{P};H\right)$ as the collection
of all $H$-valued $\mathcal{G}$-measurable square-integrable random variables.
The expectation of $u\in L^{2}\left(\Omega,\mathcal{G},\mathbb{P};H\right)$ is
denoted by

$\mathbb{E}\left[u\right]=\int_{\Omega}u\left(\omega\right)\mathrm{d}\mathbb{P}\left(\omega\right).$
(A5)

The inner produce for $u,v\in
L^{2}\left(\Omega,\mathcal{G},\mathbb{P};H\right)$ can be defined as
$\left\langle u,v\right\rangle_{2}=\mathbb{E}\left[\left\langle
u,v\right\rangle_{H}\right]=\int_{\Omega}\left\langle
u\left(\omega\right),v\left(\omega\right)\right\rangle_{H}\mathrm{d}\mathbb{P}\left(\omega\right)$.
With the above notations, the _covariance operator_ $\mathcal{C}$ can be
introduced as an element in the linear transformations
$\mathcal{L}\left(H,H\right)$.

######  Definition 13.

Let $u,v\in L^{2}\left(\Omega,\mathcal{G},\mathbb{P};H\right)$ be two
$H$-valued random variables. Then the covariance of $u$ and $v$ is given by

$\mathcal{C}\left(u,v\right)=\mathbb{E}\left[u\otimes
v\right]-\mathbb{E}\left[u\right]\otimes\mathbb{E}\left[v\right],$ (A6)

where $u\left(\omega\right)\otimes
v\left(\omega\right)\in\mathcal{L}\left(H,H\right)$ is a linear transformation
of $H$ into $H$ defined for any $f\in H$ as

$\left(u\otimes v\right)f=u\left\langle v,f\right\rangle_{H}.$

It is easy to check that the adjoint
$\mathcal{C}\left(u,v\right)^{*}=\mathcal{C}\left(v,u\right)$ and
$\mathcal{C}\left(u,u\right)^{*}=\mathcal{C}\left(u,u\right)$ is self-adjoint
since

$\left\langle f,\left(u\otimes v\right)g\right\rangle_{H}=\left\langle
u,f\right\rangle_{H}\left\langle v,g\right\rangle_{H}$

Notice that if $H=\mathbb{R}^{d}$ is finite-dimensional, for any $x,y\in
L^{2}\left(\Omega,\mathcal{G},\mathbb{P};\mathbb{R}^{d}\right)$, $x\otimes
y=xy^{\intercal}\in\mathbb{R}^{d\times d}$, then the covariance
$\mathcal{C}\left(x,y\right)\in\mathbb{R}^{d\times d}$ becomes the $d\times d$
matrix

$\mathcal{C}\left(x,y\right)=\mathbb{E}\left[xy^{\intercal}\right]-\mathbb{E}\left[x\right]\mathbb{E}\left[y\right]^{\intercal}=\mathbb{E}\left[\left(x-\mathbb{E}\left[x\right]\right)\left(y-\mathbb{E}\left[y\right]\right)^{\intercal}\right].$

Then, we call the $u_{t}\left(\omega\right)$ from
$\left[0,T\right]\times\Omega$ to $H$ an $H$-valued stochastic process. An
infinite-dimensional $H$-valued Wiener process $W_{t}$ can be defined
accordingly and the Itô integral can be generalized to infinite-dimensional
Hilbert space accordingly (see Chapter 2 of [44] with precise validations).
Therefore, the signal process $u_{t}$ of filtering can be given by the
following $H$-valued SDE

$\mathrm{d}u_{t}=\mathcal{A}_{t}u_{t}\mathrm{d}t+\mathcal{Q}_{t}\mathrm{d}W_{t},\quad
u_{t=0}=u_{0},$ (A7)

where $\mathcal{A}_{t}\in
L^{\infty}\left(\left[0,T\right];\mathcal{L}\left(H,H\right)\right)$ is a
regulated mapping of $\left[0,T\right]$ into $\mathcal{L}\left(H,H\right)$
(which is further generalized to unbounded operators in [11]) and
$\mathcal{Q}_{t}\in
L^{2}\left(\left[0,T\right];\mathcal{L}\left(H,H\right)\right)$, and $W_{t}$
is the $H$-valued Wiener process. It can be shown (Theorem 2.13 in [44] and
Theorem 5.1 in [20]) that (A7) has unique solution in
$L^{2}\left(\Omega,\mathcal{G},\mathbb{P};H\right)$ with initial value
$\mathbb{E}\left[\left|u_{0}\right|_{H}^{2}\right]<\infty$. The observation
stochastic process $y_{t}\in\mathbb{R}^{p}$ can be generated by the SDE

$\mathrm{d}y_{t}=\mathcal{H}_{t}u_{t}\mathrm{d}t+\mathrm{d}B_{t},\quad
y_{t=0}=y_{0},$ (A8)

where $\mathcal{H}_{t}$ is a linear mapping of $\left[0,T\right]$ into
$\mathcal{L}\left(H,\mathbb{R}^{p}\right)$, and $B_{t}$ is the Wiener process
in $\mathbb{R}^{p}$ independent of $W_{t}$. The infinite-dimensional filtering
problem can be then described as: given $y_{s},s\leq t$, determine the optimal
estimate $\hat{u}_{t}$ of $u_{t}$ that minimizes

$\mathbb{E}\left[\left|u_{t}-v\right|_{H}^{2}\right],\quad v\in
L^{2}\left(\Omega,\mathcal{G}_{t},\mathbb{P};H\right),$ (A9)

where $\mathcal{G}_{t}=\sigma\left\\{y_{s},s\leq t\right\\}$ is generated by
the observations up to time $t$.

Finally, in parallel to the Kalman-Bucy filter (A4) in finite-dimensional
space, similar result can be extended to the above infinite-dimensional
filtering problem (A7) and (A8). Below, we summarize the main results in
Theorem 7.10, 7.14 of [20].

######  Theorem 14.

The optimal filter solution $\hat{u}_{t}$ of (A9) exists and is unique, which
satisfies the following infinite-dimensional SDE

$\mathrm{d}\hat{u}_{t}=\mathcal{A}_{t}\hat{u}_{t}\mathrm{d}t+\mathcal{K}_{t}\left(\mathrm{d}y_{t}-\mathcal{H}_{t}\hat{u}_{t}\mathrm{d}t\right),\quad\hat{u}_{t=0}=u_{0},$
(A10)

where
$\mathcal{K}_{t}=\mathcal{C}_{t}\mathcal{H}_{t}^{*}\in\mathcal{L}\left(\mathbb{R}^{p},H\right)$
with $\mathcal{H}_{t}^{*}$ the adjoint of $\mathcal{H}_{t}$. And the
covariance operator
$\mathcal{C}_{t}=\mathbb{E}\left[\left(u_{t}-\hat{u}_{t}\right)\otimes\left(u_{t}-\hat{u}_{t}\right)\right]$
satisfies the following Riccati equation

$\dot{\mathcal{C}}_{t}=\mathcal{A}_{t}\mathcal{C}_{t}+\mathcal{C}_{t}\mathcal{A}_{t}^{*}-\mathcal{K}_{t}\mathcal{K}_{t}^{*}+\mathcal{Q}_{t}\mathcal{Q}_{t}^{*},\quad\mathcal{C}_{t=0}=\mathcal{C}\left(u_{0},u_{0}\right).$
(A11)

###  A.3 Kalman-Bucy filter with conditional Gaussian processes

The linear Kalman-Bucy filter can be generalized to nonlinear filtering
accepting the conditional Gaussian processes [36]. The conditional Gaussian
process $\left(v_{t},y_{t}\right),0\leq t\leq T$ is given by the solution of
the following coupled equations

$\displaystyle\mathrm{d}v_{t}$
$\displaystyle=\left[F_{t}\left(y_{t}\right)v_{t}+f_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Sigma_{t}\mathrm{d}W_{t},\quad
v_{t=0}=v_{0},$ (A12a) $\displaystyle\mathrm{d}y_{t}$
$\displaystyle=\left[H_{t}\left(y_{t}\right)v_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Gamma_{t}\mathrm{d}B_{t},\quad
y_{t=0}=y_{0},$ (A12b)

where $W_{t},B_{t}$ are mutually independent standard Gaussian white noise
processes, and the initial states $\left(v_{0},y_{0}\right)$ are random
variables independent of $W_{t},B_{t}$. In general, $v_{t}\in\mathbb{R}^{d}$
represents the signal process and $y_{t}\in\mathbb{R}^{p}$ represents the
observation process. The functions $f_{t},F_{t}$ and $h_{t},H_{t}$ are
globally Lipschitz continuous and uniformly bounded on the observed state
$y_{t}$ over the time interval $0\leq t\leq T$. Assume that the sequence
$\left(v_{t},y_{t}\right)$ is obtained from a realization $\omega$ and the
initial condition, $\left\\{v_{0},y_{0}\right\\}$. We can then define the
observation sequence $Y_{t}=\left\\{y_{s}\left(\omega\right),s\leq t\right\\}$
as well as the unobserved target signal $v_{t}=v_{t}\left(\omega\right)$. The
above system (A12) is called the conditional Gaussian process since the
conditional distribution $\mu_{t}=\mathbb{P}\left(v_{t}\in\cdot\mid
Y_{t}\right)$ given $Y_{t}$ becomes a Gaussian distribution a.s. if
$\mu_{0}=\mathbb{P}\left(v_{0}\in\cdot\mid y_{0}\right)$ is Gaussian (Theorem
12.6 of [36]).

Next, let $\mathcal{G}_{t}=\sigma\left\\{y_{s},s\leq t\right\\}$, and define
the mean $\hat{v}_{t}=\mathbb{E}\left[v_{t}\mid\mathcal{G}_{t}\right]$ and
covariance
$\hat{C}_{t}=\mathbb{E}\left[\left(v_{t}-\hat{v}_{t}\right)\left(v_{t}-\hat{v}_{t}\right)^{\intercal}\mid\mathcal{G}_{t}\right]$
w.r.t. the conditional Gaussian distribution
$\mu_{t}=\mathcal{N}\left(\hat{v}_{t},\hat{C}_{t}\right)$ . Then, it shows
that the explicit dynamical equations for
$\left(\hat{v}_{t},\hat{C}_{t}\right)$ can be derived based on the conditional
Gaussian process (A12). As a result, filtering equations from the linear
Kalman-Bucy filter (A4) can be directly applied to the conditional linear
system regardless of its essentially nonlinear dynamics. The equations of the
conditional Gaussian processes and the uniqueness of the solutions are proved
under suitable conditions for the model coefficients in Chapter 12 of [36]. We
summarize the results according to Theorem 12.7 of [36] for the nonlinear
conditional Gaussian filter.

######  Theorem 15.

The conditional distribution $\mu_{t}$ of the stochastic processes $v_{t}$
given $Y_{t}$ from (A12) is Gaussian,
$\mathcal{N}\left(\hat{v}_{t},\hat{C}_{t}\right)$. Then, with
$\Gamma_{t}\Gamma_{t}^{\intercal}\succ 0$ and the initial mean and covariance
$\hat{v}_{0},\hat{C}_{0}$, the solutions for the mean $\hat{v}_{t}$ and
covariance matrix $\hat{C}_{t}$ are uniquely given by the following closed
equations

$\displaystyle\mathrm{d}\hat{v}_{t}$
$\displaystyle=\left[F_{t}\left(y_{t}\right)\hat{v}_{t}+f_{t}\left(y_{t}\right)\right]\mathrm{d}t$
$\displaystyle+K_{t}\left(y_{t}\right)\left(\Gamma_{t}\Gamma_{t}^{\intercal}\right)^{-1}\left\\{\mathrm{d}y_{t}-\left[H_{t}\left(y_{t}\right)\hat{v}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\},$
(A13a) $\displaystyle\mathrm{d}\hat{C}_{t}$
$\displaystyle=\left[F_{t}\left(y_{t}\right)\hat{C}_{t}+\hat{C}_{t}F_{t}\left(y_{t}\right)^{\intercal}+\Sigma_{t}\Sigma_{t}^{\intercal}\right]\mathrm{d}t$
$\displaystyle-
K_{t}\left(y_{t}\right)\left(\Gamma_{t}\Gamma_{t}^{\intercal}\right)^{-1}K_{t}^{\intercal}\left(y_{t}\right)\mathrm{d}t,$
(A13b)

where
$K_{t}\left(y_{t}\right)=\hat{C}_{t}\left(y_{t}\right)H_{t}\left(y_{t}\right)^{\intercal}$.
The matrix $\hat{C}_{t}$ will remain positive-definite for all $0\leq t\leq T$
if $\hat{C}_{0}\succ 0$.

##  Appendix B Detailed proofs of theorems

###### Proof of Lemma 1.

First, consider the model (2.7) without the relaxation term in $R_{t}$, that
is, set $\epsilon^{-1}=0$. Applying Itô’s formula for $Z_{t}$ with any test
function $\varphi\in C_{b}^{2}\left(\mathbb{R}^{d}\right)$ gives

$\displaystyle\mathrm{d}\varphi\left(Z_{t}\right)=$
$\displaystyle\mathcal{L}\left(\bar{u}_{t},R_{t}\right)\varphi\left(Z_{t}\right)\mathrm{d}t+\nabla\varphi\left(Z_{t}\right)^{\intercal}\Sigma_{t}\mathrm{d}W_{t}$
$\displaystyle=$
$\displaystyle\nabla\varphi\left(Z_{t}\right)^{\intercal}\left[L\left(\bar{u}_{t}\right)Z_{t}+\varGamma\left(Z_{t}Z_{t}^{\intercal}-R_{t}\right)\right]\mathrm{d}t$
(B1)
$\displaystyle+\frac{1}{2}\Sigma_{t}\Sigma_{t}^{\intercal}:\nabla\nabla\varphi\left(Z_{t}\right)\mathrm{d}t+\nabla\varphi\left(Z_{t}\right)^{\intercal}\Sigma_{t}\mathrm{d}W_{t},$
(B2)

where $\mathcal{L}$ is the generator of $Z_{t}$. Given any statistical
solution $\left(\bar{u}_{t},R_{t}\right)$ and taking expectation using
$\varphi\left(Z\right)=Z$, the equation for the first moment of $Z_{t}$ can be
found as

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left[Z_{t}\right]=\left[L\left(\bar{u}_{t}\right)\mathbb{E}\left[Z_{t}\right]+\varGamma\left(\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]-R_{t}\right)\right].$
(B3)

Next by taking $\varphi\left(Z\right)=Z_{k}Z_{l}$, we have

$\displaystyle\mathrm{d}\left(Z_{k.t}Z_{l,t}\right)=$
$\displaystyle\sum_{m}\left[L_{km}\left(\bar{u}_{t}\right)Z_{m.t}Z_{l,t}+Z_{k.t}Z_{m,t}L_{lm}\left(\bar{u}_{t}\right)\right]\mathrm{d}t+\Sigma_{km,t}\Sigma_{lm,t}\mathrm{d}t$
$\displaystyle+\sum_{m,n}\gamma_{mnk}\left(Z_{m,t}Z_{n,t}Z_{l,t}-R_{mn,t}Z_{l,t}\right)\mathrm{d}t+\gamma_{mnl}\left(Z_{m,t}Z_{n,t}Z_{k,t}-R_{mn,t}Z_{k,t}\right)\mathrm{d}t$
$\displaystyle+\sum_{m}\Sigma_{km,t}Z_{l,t}\mathrm{d}W_{m,t}+\Sigma_{lm,t}Z_{k,t}\mathrm{d}W_{m,t}.$

This implies the second moment equation of $Z_{t}$ as

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]=$
$\displaystyle\sum_{m}\left[L_{km}\left(\bar{u}_{t}\right)\mathbb{E}\left[Z_{m.t}Z_{l,t}\right]+\mathbb{E}\left[Z_{k.t}Z_{m,t}\right]L_{lm}\left(\bar{u}_{t}\right)\right]\mathrm{d}t+\Sigma_{t}\Sigma_{t}^{\intercal}$
(B4)
$\displaystyle+\sum_{m,n}\gamma_{mnk}\left(\mathbb{E}\left[Z_{m,t}Z_{n,t}Z_{l,t}\right]-R_{mn,t}\mathbb{E}\left[Z_{l,t}\right]\right)\mathrm{d}t+\gamma_{mnl}\left(\mathbb{E}\left[Z_{m,t}Z_{n,t}Z_{k,t}\right]-R_{mn,t}\mathbb{E}\left[Z_{k,t}\right]\right)\mathrm{d}t.$

Assuming $\mathbb{E}\left[Z_{t}\right]=0$ and
$\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]=R_{t}$ for any time instant
$t$, first we see that the right hand side of (B3) will always stay zero.
Then, with the same statistics in the third moments of $Z_{t}$, the right hand
side of (B4) becomes equal to the right hand of the statistical equation of
$R_{t}$ in (2.7) with (2.8). Uniqueness of solution in the statistical
equations given the same initial values implies that the leading two moments
of $Z_{t}$ satisfy $\mathbb{E}\left[Z_{t}\right]=0$ and
$\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]=R_{t}$ for all $t>0$.

Finally, by adding the additional relaxation term,
$\epsilon^{-1}\left(\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]-R_{t}\right)$
in the dynamics of $R_{t},$ the statistical consistency condition (2.9)
guarantees $R_{t}=\mathbb{E}\left[Z_{t}Z_{t}^{\intercal}\right]$ for all the
time. Thus this term gives zero contribution, so the model (2.7) will also
produce the same statistical solution. ∎

###### Proof of Proposition 2.

First, consider $\mathbb{E}\varphi\left(u_{t}\right)$ with any test function
$\varphi\in C_{b}^{2}\left(\mathbb{R}^{d}\right)$ w.r.t. the PDF $p_{t}$ for
the state of the original system (1.1). Itô’s lemma shows that

$\frac{\mathrm{d}\mathbb{E}\varphi\left(u_{t}\right)}{\mathrm{d}t}=\mathbb{E}\left[\left(\Lambda
u_{t}+B\left(u_{t},u_{t}\right)+F_{t}\right)\cdot\nabla\varphi\left(u_{t}\right)\right]+\frac{1}{2}\mathbb{E}\left[\sigma_{t}\sigma_{t}^{\intercal}:\nabla\nabla\varphi\left(u_{t}\right)\right],$

where $A:\nabla\nabla f=\sum_{mn}a_{mn}\partial_{u_{m}}\partial_{u_{n}}f$. By
taking $\varphi=u$ and introducing the decomposition
$u^{\prime}=u-\mathbb{E}u=\sum_{k}u_{k,t}^{\prime}\hat{v}_{k}$ with
$u_{k,t}^{\prime}=\hat{v}_{k}\cdot u_{t}^{\prime}$, we have

$\displaystyle\frac{\mathrm{d}\mathbb{E}u_{t}}{\mathrm{d}t}$
$\displaystyle=\mathbb{E}\left[\Lambda u_{t}+B\left(u,u\right)+F_{t}\right]$
$\displaystyle=\Lambda\mathbb{E}u_{t}+B\left(\mathbb{E}u,\mathbb{E}u\right)+\sum_{k,l}\mathbb{E}\left[u_{k,t}^{\prime}u_{l,t}^{\prime}\right]B\left(\hat{v}_{k},\hat{v}_{l}\right)+F_{t}.$
(B5)

Above, we use the bilinearity of the operator $B$ and
$\mathbb{E}u^{\prime}=0$, such that

$\displaystyle\mathbb{E}B\left(u,u\right)$
$\displaystyle=\mathbb{E}B\left(\mathbb{E}u+u^{\prime},\mathbb{E}u+u^{\prime}\right)$
$\displaystyle=\mathbb{E}\left[B\left(\mathbb{E}u,\mathbb{E}u\right)+B\left(\mathbb{E}u,u^{\prime}\right)+B\left(u^{\prime},\mathbb{E}u\right)+B\left(u^{\prime},u^{\prime}\right)\right]$
$\displaystyle=B\left(\mathbb{E}u,\mathbb{E}u\right)+\mathbb{E}\sum_{k,l}B\left(u_{k,t}^{\prime}\hat{v}_{k},u_{l,t}^{\prime}\hat{v}_{l}\right)$
$\displaystyle=B\left(\mathbb{E}u,\mathbb{E}u\right)+\sum_{k,l}\mathbb{E}\left[u_{k,t}^{\prime}u_{l,t}^{\prime}\right]B\left(\hat{v}_{k},\hat{v}_{l}\right).$

Similarly, by taking $\varphi=\left(\hat{v}_{k}\cdot
u_{t}^{\prime}\right)\left(u_{t}^{\prime}\cdot\hat{v}_{l}\right)=u_{t}^{\prime\intercal}Au_{t}^{\prime}$
with $A=\hat{v}_{l}\hat{v}_{k}^{\intercal}$, we find

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left[u_{k,t}^{\prime}u_{l,t}^{\prime}\right]=$
$\displaystyle\mathbb{E}\left[\left(\Lambda
u_{t}+B\left(u_{t},u_{t}\right)+F_{t}\right)\cdot\left(A+A^{\intercal}\right)u_{t}^{\prime}\right]+\frac{1}{2}\sigma_{t}\sigma_{t}^{\intercal}:\left(A+A^{\intercal}\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(u_{t}^{\prime}\cdot\hat{v}_{l}\right)\left(\hat{v}_{k}^{\intercal}\Lambda
u_{t}^{\prime}\right)+\left(u_{t}^{\prime\intercal}\Lambda^{\intercal}\hat{v}_{l}\right)\left(\hat{v}_{k}\cdot
u_{t}^{\prime}\right)\right]$ $\displaystyle+\mathbb{E}\left[\hat{v}_{l}\cdot
B\left(\mathbb{E}u,u^{\prime}\right)u_{k,t}^{\prime}+\hat{v}_{k}\cdot
B\left(\mathbb{E}u,u^{\prime}\right)u_{l,t}^{\prime}\right]+\mathbb{E}\left[\hat{v}_{l}\cdot
B\left(u^{\prime},\mathbb{E}u\right)u_{k,t}^{\prime}+\hat{v}_{k}\cdot
B\left(u^{\prime},\mathbb{E}u\right)u_{l,t}^{\prime}\right]$
$\displaystyle+\mathbb{E}\left[B\left(u_{t}^{\prime},u_{t}^{\prime}\right)\cdot\left(\hat{v}_{l}u_{k,t}^{\prime}+\hat{v}_{k}u_{l,t}^{\prime}\right)\right]+\left(\sigma_{t}\cdot\hat{v}_{k}\right)\left(\hat{v}_{l}\cdot\sigma_{t}\right)$
$\displaystyle=$
$\displaystyle\sum_{m}\left(\hat{v}_{k}^{\intercal}\Lambda\hat{v}_{m}\right)\mathbb{E}\left[u_{m,t}^{\prime}u_{l,t}^{\prime}\right]+\mathbb{E}\left[u_{k,t}^{\prime}u_{m,t}^{\prime}\right]\left(\hat{v}_{l}^{\intercal}\Lambda^{\intercal}\hat{v}_{m}\right)+\left(\sigma_{t}\cdot\hat{v}_{k}\right)\left(\hat{v}_{l}\cdot\sigma_{t}\right)$
(B6)
$\displaystyle+\sum_{m}\left[\hat{v}_{k}^{\intercal}B\left(\mathbb{E}u,\hat{v}_{m}\right)\mathbb{E}\left[u_{m,t}^{\prime}u_{l,t}^{\prime}\right]+\mathbb{E}\left[u_{k,t}^{\prime}u_{m,t}^{\prime}\right]\hat{v}_{l}^{\intercal}B\left(\mathbb{E}u,\hat{v}_{m}\right)\right]$
$\displaystyle+\sum_{m}\left[\hat{v}_{k}^{\intercal}B\left(\hat{v}_{m},\mathbb{E}u\right)\mathbb{E}\left[u_{m,t}^{\prime}u_{l,t}^{\prime}\right]+\mathbb{E}\left[u_{k,t}^{\prime}u_{m,t}^{\prime}\right]\hat{v}_{l}^{\intercal}B\left(\hat{v}_{m},\mathbb{E}u\right)\right]$
$\displaystyle+\sum_{m,n}\left[\hat{v}_{k}^{\intercal}B\left(\hat{v}_{m},\hat{v}_{n}\right)\mathbb{E}\left[u_{l,t}^{\prime}u_{m,t}^{\prime}u_{n,t}^{\prime}\right]+\hat{v}_{l}^{\intercal}B\left(\hat{v}_{m},\hat{v}_{n}\right)\mathbb{E}\left[u_{k,t}^{\prime}u_{m,t}^{\prime}u_{n,t}^{\prime}\right]\right].$

Above in the second equality, we use the projection of modes
$\left(A+A^{\intercal}\right)u_{t}^{\prime}=\left(\hat{v}_{l}\hat{v}_{k}^{\intercal}+\hat{v}_{k}\hat{v}_{l}^{\intercal}\right)u_{t}^{\prime}=\hat{v}_{l}u_{k,t}^{\prime}+\hat{v}_{k}u_{k,t}^{\prime}$,
and the bilinearity of the quadratic operator $B$; and in the third equality,
we use the decomposition of fluctuation modes,
$u_{t}^{\prime}=\sum_{k}u_{k,t}^{\prime}\hat{v}_{k}$. We find the coupling
operator
$L_{km}=\left(\hat{v}_{k}^{\intercal}\Lambda\hat{v}_{m}\right)+\hat{v}_{k}^{\intercal}B\left(\hat{v}_{m},\mathbb{E}u\right)+\hat{v}_{k}^{\intercal}B\left(\hat{v}_{m},\mathbb{E}u\right)$,
the third-order coupling coefficients
$\gamma_{mnk}=\hat{v}_{k}^{\intercal}B\left(\hat{v}_{m},\hat{v}_{n}\right)$,
as well as the noise term
$\left(\sigma_{t}\cdot\hat{v}_{k}\right)\left(\hat{v}_{l}\cdot\sigma_{t}\right)$.

In addition, by subtracting the mean equation (B5) from the original system
(1.1), we find the SDE for the stochastic state

$\displaystyle\frac{\mathrm{d}u_{t}^{\prime}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\left(u_{t}-\mathbb{E}u_{t}\right)=$
$\displaystyle\Lambda u_{t}+B\left(u_{t},u_{t}\right)+\sigma_{t}\dot{W}_{t}$
$\displaystyle-\Lambda\mathbb{E}u_{t}-B\left(\mathbb{E}u,\mathbb{E}u\right)-\sum_{m,n}\mathbb{E}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}\right]B\left(\hat{v}_{m},\hat{v}_{n}\right)$
$\displaystyle=$
$\displaystyle\sum_{m}u_{m,t}^{\prime}\left[\Lambda\hat{v}_{m}+B\left(\mathbb{E}u,\hat{v}_{m}\right)+B\left(\hat{v}_{m},\mathbb{E}u\right)\right]$
$\displaystyle+\sum_{m,n}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}-\mathbb{E}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}\right]\right]B\left(\hat{v}_{m},\hat{v}_{n}\right)+\sigma_{t}\dot{W}_{t}.$

Again in the second equality, we use the spectral decomposition of the
fluctuation state $u_{k,t}^{\prime}=u_{t}^{\prime}\cdot\hat{v}_{k}$. By
projecting the state $u_{t}^{\prime}$ on the basis $\hat{v}_{k}$, we have

$\frac{\mathrm{d}u_{k,t}^{\prime}}{\mathrm{d}t}=\sum_{k}L_{km}\left(\mathbb{E}u_{t}\right)u_{m,t}^{\prime}+\sum_{m,n}\gamma_{mnk}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}-\mathbb{E}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}\right]\right]+\hat{v}_{k}^{\intercal}\sigma_{t}\dot{W}_{t},$
(B7)

with the same parameters $L_{km}$ and $\gamma_{mnk}$ defined before. The
generator $\mathcal{L}_{t}^{u}$ of $u_{t}^{\prime}$ can be written as

$\mathcal{L}_{t}^{u}\left(p_{t}\right)=\left[\sum_{k}L_{km}\left(\mathbb{E}u_{t}\right)u_{m,t}^{\prime}+\sum_{m,n}\gamma_{mnk}\left(u_{m,t}^{\prime}u_{n,t}^{\prime}-\mathbb{E}\left[u_{m,t}^{\prime}u_{n,t}^{\prime}\right]\right)\right]\cdot\nabla_{u^{\prime}}+\frac{1}{2}\Sigma_{t}\Sigma_{t}^{\intercal}:\nabla_{u^{\prime}}\nabla_{u^{\prime}},$

where $\mathcal{L}_{t}^{u}$ is dependent on $p_{t}$ in computing the
expectations.

Next, we consider the closure model (2.7) without the relaxation term

$\displaystyle\frac{\mathrm{d}\bar{u}_{t}}{\mathrm{d}t}=\Lambda\bar{u}_{t}+B\left(\bar{u}_{t},\bar{u}_{t}\right)$
$\displaystyle+Q_{m}\left(\mathbb{E}\left[Z_{t}\otimes
Z_{t}\right]\right)+F_{t},$ (B8)
$\displaystyle\frac{\mathrm{d}R_{t}}{\mathrm{d}t}=L\left(\bar{u}_{t}\right)R_{t}+R_{t}L\left(\bar{u}_{t}\right)$
$\displaystyle+Q_{v}\left(\mathbb{E}\left[Z_{t}\otimes Z_{t}\otimes
Z_{t}\right]\right)+\Sigma_{t}\Sigma_{t}^{\intercal},$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left[\varphi\left(Z_{t}\right)\right]=$
$\displaystyle\>\mathbb{E}\left[\mathcal{L}_{t}\left(\bar{u}_{t},R_{t}\right)\varphi\left(Z_{t}\right)\right].$

where $\mathcal{L}_{t}$ is the generator from (2.15) defined from the McKean-
Vlasov SDE of $Z_{t}$

$\mathcal{L}_{t}\left(\bar{u}_{t},R_{t}\right)=\left[L\left(\bar{u}_{t}\right)Z_{t}+\varGamma\left(Z_{t}Z_{t}^{\intercal}-R_{t}\right)\right]\cdot\nabla_{z}+\frac{1}{2}\Sigma_{t}\Sigma_{t}^{\intercal}:\nabla_{z}\nabla_{z}.$

The closure terms $Q_{m}$ and $Q_{v}$ in (2.8) have exactly the same structure
as that in the original system derived in (B5) and (B6). In addition, by
comparing the above SDEs (B7) and (2.5), it is realized that their generators,
$\mathcal{L}_{t}^{u}\left(p_{t}\right)$ and
$\mathcal{L}_{t}\left(\bar{u}_{t},R_{t}\right)$, share the same dynamical
structure with the dependence on the first two moments w.r.t. $p_{t}$ and the
statistical solutions $\bar{u}_{t},R_{t}$ in (B8). Therefore, at any time
instant $t$ if we assume consistent statistics

$\mathbb{E}_{p_{t}}\left[u_{t}\right]=\bar{u}_{t},\quad\mathbb{E}_{p_{t}}\left[\left(u_{t}^{\prime}\cdot\hat{v}_{k}\right)\left(u_{t}^{\prime}\cdot\hat{v}_{l}\right)\right]=R_{kl,t},$

as well as

$\mathbb{E}_{p_{t}}\left[\varphi\left(u_{t}^{\prime}\right)\right]=\mathbb{E}\left[\varphi\left(\sum_{k=1}^{d}Z_{k,t}\hat{v}_{k}\right)\right],$

the right hand sides of the original model (B5), (B6), and (B7) become the
same as that of the closure model (B8). Starting from the same initial
condition with uniqueness of the solution, it directly implies that the
statistical solutions of the two systems (1.1) and (B8) will remain the same
during the entire time evolution. ∎

###### Proof of Lemma 4.

We rewrite the filter model (3.18) for $\tilde{Z}_{t}$ by substituting the
explicit equation for the observation process,
$\mathrm{d}y_{t}=\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t+\Gamma_{t}\mathrm{d}B_{t}$,
in (3.3)

$\mathrm{d}\tilde{Z}_{t}=a_{t}\left(\tilde{Z}_{t}\right)\mathrm{d}t+K_{t}\left(\tilde{Z}_{t}\right)\left\\{\left[\mathcal{H}\rho_{t}-H\left(\tilde{Z}_{t}\right)\right]\mathrm{d}t+\Gamma_{t}\mathrm{d}B_{t}-\Gamma_{t}\mathrm{d}\tilde{B}_{t}\right\\}.$

By applying Itô’s formula on the above SDE, we have for $\varphi\in
C_{b}^{2}\left(\mathbb{R}^{d}\right)$

$\displaystyle\mathrm{d}\varphi\left(\tilde{Z}_{t}\right)=$
$\displaystyle\nabla\varphi\cdot\left[\left(a_{t}-K_{t}\left(H\left(\tilde{Z}_{t}\right)+h_{t}\left(y_{t}\right)\right)\right)\right]\mathrm{d}t$
$\displaystyle-\nabla\varphi\cdot
K_{t}\Gamma_{t}\mathrm{d}\tilde{B}_{t}+\nabla\varphi\cdot
K_{t}\mathrm{d}y_{t}+K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla\varphi\mathrm{d}t,$
(B9)

where we define
$A:\nabla\nabla\varphi=\sum_{i,j=1}^{d}A_{ij}\partial_{z_{i}z_{j}}\varphi$ and
take the convention $\left(\nabla f\right)_{ij}=\partial_{z_{i}}f_{j}$ for the
gradient of vector-valued functions $f\in
C^{1}\left(\mathbb{R}^{d};\mathbb{R}^{p}\right)$. Notice that above the
coefficient in the last term is $1$ considering the additional contributions
from the independent white noise process
$\Gamma_{t}\mathrm{d}B_{t}=\mathrm{d}y_{t}-\left[\mathcal{H}\rho_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t$
in the observation process besides the original $\mathrm{d}\tilde{B}_{t}$,
that is,

$\frac{1}{2}\nabla\nabla\varphi:\mathrm{d}\left\langle
K\Gamma\tilde{B},K\Gamma\tilde{B}\right\rangle_{t}+\frac{1}{2}\nabla\nabla\varphi:\mathrm{d}\left\langle
K\Gamma B,K\Gamma
B\right\rangle_{t}=\nabla\nabla\varphi:K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\mathrm{d}t,$

where we denote $\left\langle M,N\right\rangle_{t}$ as the Meyer’s process of
two martingales $M_{t}$ and $N_{t}$.

First, by taking $\varphi\left(z\right)=H\left(z\right)$ and taking
expectation $\tilde{\mathbb{E}}$ w.r.t. $\tilde{\rho}_{t}$ conditional on
$Y_{t}=\left\\{y_{s},s\leq t\right\\}\in\mathcal{G}_{t}$, we have

$\displaystyle\mathrm{d}\tilde{\mathbb{E}}H\left(\tilde{Z}_{t}\right)=$
$\displaystyle\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}a_{t}\right]\mathrm{d}t-\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}K_{t}\left(\bar{H}_{t}+H_{t}^{\prime}+h_{t}\left(y_{t}\right)\right)\right]\mathrm{d}t$
$\displaystyle+\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}K_{t}\right]\mathrm{d}y_{t}+\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla
H\left(\tilde{Z}_{t}\right)\right]\mathrm{d}t$ $\displaystyle=$
$\displaystyle\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}a_{t}\right]\mathrm{d}t+\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla
H\left(\tilde{Z}_{t}\right)\right]\mathrm{d}t-\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}K_{t}H_{t}^{\prime}\right]\mathrm{d}t$
(B10) $\displaystyle+\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}K_{t}\right]\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\}.$

In the first line above, we split
$H\left(\tilde{Z}_{t}\right)=\bar{H}_{t}+H_{t}^{\prime}$. Notice that the
observation process $y_{t}\in\mathcal{G}_{t}$ can be brought out of the
expectation
$\tilde{\mathbb{E}}\left[\cdot\right]=\mathbb{E}\left[\cdot\mid\mathcal{G}_{t}\right]$.
Using the first identity in (3.21) for $a_{t}$, there is

$\displaystyle\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}a_{t}\right]=$ $\displaystyle\int\nabla
H\left(z\right)^{\intercal}\left[\nabla\cdot\left(K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right)-K_{t}\Gamma_{t}^{2}\nabla\cdot
K_{t}^{\intercal}\right]\tilde{\rho}_{t}\left(z\right)\mathrm{d}z$
$\displaystyle=$ $\displaystyle\int\nabla
H\left(z\right)^{\intercal}\left[\nabla\cdot\left(\tilde{\rho}_{t}K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right)-K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\nabla\tilde{\rho}_{t}-K_{t}\Gamma_{t}^{2}\nabla\cdot
K_{t}^{\intercal}\tilde{\rho}_{t}\right]\mathrm{d}z$ $\displaystyle=$
$\displaystyle-\int\nabla\nabla
H\left(z\right):\left(K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right)\tilde{\rho}_{t}\mathrm{d}z-\int\nabla
H\left(z\right)^{\intercal}K_{t}\Gamma_{t}^{2}\nabla\cdot\left(\tilde{\rho}_{t}K_{t}^{\intercal}\right)\mathrm{d}z$
$\displaystyle=$
$\displaystyle-\tilde{\mathbb{E}}\left[\left(K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}\right):\nabla\nabla
H\left(\tilde{Z}_{t}\right)\right]-\int\nabla
H\left(z\right)^{\intercal}K_{t}\Gamma_{t}^{2}\nabla\cdot\left(\tilde{\rho}_{t}K_{t}^{\intercal}\right)\mathrm{d}z.$

Then using the second identity in (3.21) for $K_{t}$ and denoting
$H_{t}^{\prime}=H-\bar{H}_{t}$, the last term above gets simplified to

$-\int\nabla
H{}^{\intercal}K_{t}\Gamma_{t}^{2}\nabla\cdot\left(\tilde{\rho}_{t}K_{t}^{\intercal}\right)\mathrm{d}z=\int\nabla
H_{t}^{\prime\intercal}K_{t}\tilde{\rho}_{t}H_{t}^{\prime}\mathrm{d}z=\tilde{\mathbb{E}}\left[\nabla
H_{t}^{\prime\intercal}K_{t}H_{t}^{\prime}\right].$

With the above identities, first line of (B10) becomes zero. Further with the
second identity in (3.21) for $K_{t}$, there is

$\tilde{\mathbb{E}}\left(K_{t}^{\intercal}\nabla\psi\right)=\Gamma_{t}^{-2}\tilde{\mathbb{E}}\left[\left(H\left(\tilde{Z}_{t}\right)-\mathbb{E}H\right)\psi\left(\tilde{Z}_{t}\right)^{\intercal}\right],$

for any regular function $\psi$ with $\mathbb{E}\psi=0$. By taking
$\psi=H-\bar{H}_{t}$, there is

$\tilde{\mathbb{E}}\left[\nabla
H\left(\tilde{Z}_{t}\right)^{\intercal}K_{t}\right]=\tilde{\mathbb{E}}\left[\left(H\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right)\left(H\left(\tilde{Z}_{t}\right)-\bar{H}_{t}\right)^{\intercal}\right]\Gamma^{-2}=C_{t}^{H}\Gamma_{t}^{-2}.$

This gives the equation for
$\bar{H}_{t}=\tilde{\mathbb{E}}H\left(\tilde{Z}_{t}\right)$.

Next, we take $\varphi\left(z\right)=H_{k}\left(z\right)H_{l}\left(z\right)$.
For the convenience of computation, we separate the mean state $\bar{H}_{t}$
as

$\varphi\left(z\right)=\left[\bar{H}_{k,t}+H_{k}^{\prime}\left(z\right)\right]\left[\bar{H}_{l,t}+H_{l}^{\prime}\left(z\right)\right]=H_{k,t}^{\prime}\left(z\right)H_{l,t}^{\prime}\left(z\right)+\bar{H}_{k,t}H_{l}^{\prime}\left(z\right)+H_{k}^{\prime}\left(z\right)\bar{H}_{l,t}+\bar{H}_{k,t}\bar{H}_{l,t}.$

The last term above is independent of $z$, thus will vanish after applying
Itô’s formula (B9). We have for the first term on the right hand side

$\displaystyle\mathrm{d}\tilde{\mathbb{E}}\left[H_{k,t}^{\prime}\left(\tilde{Z}_{t}\right)H_{l,t}^{\prime}\left(\tilde{Z}_{t}\right)\right]$
$\displaystyle=\tilde{\mathbb{E}}\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}\left[\left(a_{t}-K_{t}\left(\bar{H}_{t}+H_{t}^{\prime}+h_{t}\left(y_{t}\right)\right)\right)\right]\mathrm{d}t$
$\displaystyle+\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\right]\mathrm{d}y_{t}+\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\right]\mathrm{d}t$
$\displaystyle=\tilde{\mathbb{E}}\left[\nabla\left(H_{k}^{\prime}H_{l}^{\prime}\right)^{\intercal}a_{t}\right]\mathrm{d}t+\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\right]\mathrm{d}t$
(B11)
$\displaystyle+\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\right]\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\}-\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}H_{t}^{\prime}\right].$

Using the identifies (3.21) for $a_{t}$ and $K_{t}$, again we can find

$\displaystyle\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}a_{t}\right]$
$\displaystyle=-\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\right]-\int\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\Gamma_{t}^{2}\nabla\cdot\left(\tilde{\rho}_{t}K_{t}^{\intercal}\right)\mathrm{d}z$
$\displaystyle=-\tilde{\mathbb{E}}\left[K_{t}\Gamma_{t}^{2}K_{t}^{\intercal}:\nabla\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\right]+\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}H_{t}^{\prime}\right].$

Therefore, we have
$\mathrm{d}\tilde{\mathbb{E}}\left[H_{k,t}^{\prime}H_{l,t}^{\prime}\right]=\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\right]\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\}$.
Further, using the identity for $K_{t}$, the coefficient becomes third moments
of $H_{t}^{\prime}$

$\displaystyle\tilde{\mathbb{E}}\left[\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\right]$
$\displaystyle=\int\nabla\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\tilde{\rho}_{t}\mathrm{d}z$
$\displaystyle=-\int\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\left[\nabla\cdot\left(\tilde{\rho}_{t}K_{t}^{\intercal}\right)\right]^{\intercal}\mathrm{d}z$
$\displaystyle=\int\left(H_{k,t}^{\prime}H_{l,t}^{\prime}\right)\left[\tilde{\rho}_{t}\Gamma_{t}^{-2}\left(H\left(z\right)-\bar{H}_{t}\right)\right]^{\intercal}\mathrm{d}z$
$\displaystyle=\tilde{\mathbb{E}}\left[H_{k,t}^{\prime}H_{l,t}^{\prime}H^{\prime\intercal}\right]\Gamma_{t}^{-2}.$

Similarly, by repeating the same procedure for
$\bar{H}_{k,t}H_{l}^{\prime}\left(z\right)$, we have

$\displaystyle\mathrm{d}\tilde{\mathbb{E}}\left[\bar{H}_{k,t}H_{l}^{\prime}\left(\tilde{Z}_{t}\right)\right]$
$\displaystyle=\tilde{\mathbb{E}}\left[\nabla\left(\bar{H}_{k,t}H_{l,t}^{\prime}\right)^{\intercal}K_{t}\right]\left\\{\mathrm{d}y_{t}-\left[\bar{H}_{t}+h_{t}\left(y_{t}\right)\right]\mathrm{d}t\right\\}$
$\displaystyle=\tilde{\mathbb{E}}\left[H_{l,t}^{\prime}H_{t}^{\prime\intercal}\right]\Gamma_{t}^{-2}\left[\mathrm{d}y_{t}-\left(\bar{H}_{t}+h_{t}\right)\mathrm{d}t\right]\bar{H}_{k,t}.$

And similar result can be achieved for
$\tilde{\mathbb{E}}\left[H_{k}^{\prime}\left(\tilde{Z}_{t}\right)\bar{H}_{l,t}\right]$.

Finally, applying Itô’s formula for $\bar{H}_{t}\bar{H}_{t}^{\intercal}$ where
$\mathrm{d}\bar{H}_{t}=C_{t}^{H}\Gamma_{t}^{-2}\left(\mathcal{H}\rho_{t}-\bar{H}_{t}\right)\mathrm{d}t+C_{t}^{H}\Gamma_{t}^{-1}\mathrm{d}B_{t}$
as we have derived, there is,

$\displaystyle\mathrm{d}\left(\bar{H}_{t}\bar{H}_{t}^{\intercal}\right)$
$\displaystyle=\left(\mathrm{d}\bar{H}_{t}\right)\bar{H}_{t}^{\intercal}+\bar{H}_{t}\left(\mathrm{d}\bar{H}_{t}^{\intercal}\right)+\mathrm{d}\left\langle
C^{H}\Gamma^{-1}B,C^{H}\Gamma^{-1}B\right\rangle_{t}$
$\displaystyle=C_{t}^{H}\Gamma_{t}^{-2}\left[\mathrm{d}y_{t}-\left(\bar{H}_{t}+h_{t}\right)\mathrm{d}t\right]\bar{H}_{t}^{\intercal}$
$\displaystyle+\bar{H}_{t}\left[\mathrm{d}y_{t}^{\intercal}-\left(\bar{H}_{t}+h_{t}\right)^{\intercal}\mathrm{d}t\right]\Gamma_{t}^{-2}C_{t}^{H}$
$\displaystyle+C_{t}^{H}\Gamma_{t}^{-2}C_{t}^{H}\mathrm{d}t.$

Notice again that the white noise process,
$C_{t}^{H}\Gamma_{t}^{-1}\mathrm{d}B_{t}$, gives the last term in the first
equality above. Putting all the above equations together, we get the equation
for
$\mathrm{d}C_{t}^{H}=\mathrm{d}\tilde{\mathbb{E}}\left[H\left(\tilde{Z}_{t}\right)H\left(\tilde{Z}_{t}\right)^{\intercal}\right]-\mathrm{d}\left(\bar{H}_{t}\bar{H}_{t}^{\intercal}\right)$
where
$C_{kl,t}^{H}=\tilde{\mathbb{E}}\left[H_{k,t}^{\prime}H_{l,t}^{\prime}\right]$.
∎

###### Proof of Proposition 10.

According to (5.4) with $K=\tilde{K}\Gamma^{-2}$, we need to show

$-\nabla\cdot\left(\tilde{K}^{\intercal}\tilde{\rho}\right)H^{\intercal}=\tilde{\rho}H^{\prime}H^{\intercal}\>\Rightarrow\>\tilde{\mathbb{E}}\left[\tilde{K}^{\intercal}\nabla
H\right]=\tilde{\mathbb{E}}\left[H^{\prime}\left(\bar{H}+H^{\prime}\right)^{\intercal}\right]=C^{H},$

according to the specific expressions $H=H^{m}$ and $H=H^{v}$. First, we can
compute

$\nabla_{z}H_{l}^{m}=2A_{l}z,\quad\nabla_{z}H_{pq}^{v}=2z_{q}A_{p}z+\left(z^{\intercal}A_{p}z\right)\delta_{qj}\hat{e}_{j}.$

Above in $H^{v}$ for simplicity, we only compute half of the symmetric
function and $\hat{e}_{j}$ is the unit vector with value $1$ in the $j$-th
entry.

From direct computations for $H^{m}$ and using
$H_{k}^{m}=z^{\intercal}A_{k}z$, we have

$\displaystyle\sum_{j}\tilde{K}_{j,k}^{m}\frac{\partial H_{l}^{m}}{\partial
z_{j}}$
$\displaystyle=\frac{1}{2}\left[\left(z^{\intercal}A_{k}z\right)-\bar{H}_{k}^{m}\right]\sum_{j}z_{j}2\left(A_{l}z\right)_{j}$
$\displaystyle=\left[\left(z^{\intercal}A_{k}z\right)-\bar{H}_{k}^{m}\right]\left(z^{\intercal}A_{l}z\right)=H_{k}^{m\prime}H_{l}^{m}.$

Similarly for $H^{v}$, we can compute

$\displaystyle\sum_{j}\tilde{K}_{j,kl}^{v}\frac{\partial H_{pq}^{v}}{\partial
z_{j}}$
$\displaystyle=\frac{1}{3}\left[\left(z^{\intercal}A_{k}z\right)z_{l}-\bar{H}_{kl}^{v}\right]\sum_{j}z_{j}\left[2z_{q}\left(A_{p}z\right)_{j}+\left(z^{\intercal}A_{p}z\right)\delta_{qj}\right]$
$\displaystyle=\frac{2}{3}\left[\left(z^{\intercal}A_{k}z\right)z_{l}-\bar{H}_{kl}^{v}\right]\left(z^{\intercal}A_{p}z\right)z_{q}+\frac{1}{3}\left[\left(z^{\intercal}A_{k}z\right)z_{l}-\bar{H}_{kl}^{v}\right]z_{q}\left(z^{\intercal}A_{p}z\right)$
$\displaystyle=\left[\left(z^{\intercal}A_{k}z\right)z_{l}-\bar{H}_{kl}^{v}\right]\left(z^{\intercal}A_{p}z\right)z_{q}=H_{kl}^{v\prime}H_{pq}^{v}.$

This finishes the proof. ∎

###### Proof of Proposition 12.

First, consider the mean equations from the same initial state, we have

$\displaystyle\bar{u}_{t}^{N,\tau}-\bar{u}_{t}=$
$\displaystyle\int_{0}^{t}\left[M\left(\bar{u}_{\sigma\left(s\right)}^{N,\tau}\right)-M\left(\bar{u}_{s}\right)\right]\mathrm{d}s+\int_{0}^{t}\left[F_{\sigma\left(s\right)}-F_{s}\right]\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\left[\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{s}\right\rangle\right]\mathrm{d}s.$

Therefore, using Hölder’s inequality and Lipschitz condition for $M$ there is

$\displaystyle\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{t}^{N,\tau}-\bar{u}_{t}\right|^{2}$ $\displaystyle\leq
3T\beta\mathbb{E}\int_{0}^{T}\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}\right|^{2}\mathrm{d}s$
$\displaystyle+3T\mathbb{E}\int_{0}^{T}\left|\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{s}\right\rangle\right|^{2}\mathrm{d}s+C_{T}\tau$
$\displaystyle\leq C_{1}\int_{0}^{T}\mathbb{E}\sup_{s^{\prime}\leq
s}\left|\bar{u}_{s^{\prime}}^{N,\tau}-\bar{u}_{s^{\prime}}\right|^{2}\mathrm{d}s+C_{2}\tau\left\|H^{m}\right\|_{\infty}^{2}+\frac{C_{3}}{N}\left\|H_{m}\right\|_{\infty}^{2}.$
(B12)

In the first term above on the right hand side, we follow the same procedure
to estimate the error in the corresponding continuous solution
$\bar{u}_{s}^{N,\tau}$ compared to the time discretization solution
$\bar{u}_{\sigma\left(s\right)}^{N,\tau}$

$\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}\right|^{2}\leq\left|\bar{u}_{\sigma\left(s\right)}^{N,\tau}-\bar{u}_{s}^{N,\tau}\right|^{2}+\left|\bar{u}_{s}^{N,\tau}-\bar{u}_{s}\right|^{2}\leq
C\tau\left\|H^{m}\right\|_{\infty}^{2}+\left|\bar{u}_{s}^{N,\tau}-\bar{u}_{s}\right|^{2}.$

And in the second line for the term related to $H^{m}$, combining the discrete
time estimate with (5.8) gives

$\mathbb{E}\left[\sup_{t\leq T}\left|\left\langle
H^{m},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{m},\tilde{\rho}_{t}\right\rangle\right|^{2}\right]\leq\left(C\tau^{2}+\frac{C_{T}}{N}\right)\left\|H^{m}\right\|_{\infty}^{2}.$

Therefore, applying Grönwall’s inequality to (B12), we get the mean estimate
in (5.11).

Next, under a similar fashion, we can compute from the covariance equation and
using the Lipschitz condition for $L$

$\displaystyle\mathbb{E}\sup_{t\leq T}\left|R_{t}^{N,\tau}-R_{t}\right|^{2}$
$\displaystyle\leq C_{1}\beta\mathbb{E}\sup_{t\leq
T}\left\|\bar{u}_{t}\right\|_{\infty}^{2}\int_{0}^{T}\left|R_{\sigma\left(s\right)}^{N,\tau}-R_{s}\right|^{2}\mathrm{d}s$
$\displaystyle+C_{2}\beta\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{\sigma\left(t\right)}^{N,\tau}-\bar{u}_{t}\right|^{2}\sup_{t\leq
T}\left\|R_{\sigma\left(t\right)}\right\|^{2}$
$\displaystyle+C_{3}\mathbb{E}\int_{0}^{T}\left|\left\langle
H^{v},\tilde{\rho}_{\sigma\left(s\right)}^{N}\right\rangle-\left\langle
H^{v},\tilde{\rho}_{s}\right\rangle\right|^{2}\mathrm{d}s+C_{T}\tau.$

Using the uniform boundedness of $\bar{u}_{t},R_{t}$ and (5.8) for $H^{v}$
together with the previous estimate for $\mathbb{E}\sup_{t\leq
T}\left|\bar{u}_{\sigma\left(t\right)}^{N,\tau}-\bar{u}_{t}\right|^{2}$, we
reach the final covariance estimate in (5.11). ∎

## References

  * [1] Eviatar Bach, Tim Colonius, Isabel Scherl, and Andrew Stuart.  Filtering dynamical systems using observations of statistics.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 34(3), 2024.
  * [2] Alan Bain and Dan Crisan.  Fundamentals of stochastic filtering, volume 3.  Springer, 2009.
  * [3] Abhishek Balakrishna and Animikh Biswas.  Determining map, data assimilation and an observable regularity criterion for the three-dimensional Boussinesq system.  Applied Mathematics & Optimization, 86(3):28, 2022.
  * [4] Kay Bergemann and Sebastian Reich.  An ensemble Kalman-Bucy filter for continuous data assimilation.  Meteorologische Zeitschrift, 21(3):213, 2012.
  * [5] Animikh Biswas and Michal Branicki.  A unified framework for the analysis of accuracy and stability of a class of approximate Gaussian filters for the Navier-Stokes equations.  arXiv preprint arXiv:2402.14078, 2024.
  * [6] Amarjit Budhiraja.  Asymptotic stability, ergodicity and other asymptotic properties of the nonlinear filter.  In Annales de l’IHP Probabilités et statistiques, volume 39, pages 919–941, 2003.
  * [7] Edoardo Calvello, Sebastian Reich, and Andrew M Stuart.  Ensemble Kalman methods: A mean field perspective.  Acta Numerica, 2025.
  * [8] Dan Crisan and Jie Xiong.  Approximate McKean–Vlasov representations for a class of SPDEs.  Stochastics An International Journal of Probability and Stochastics Processes, 82(1):53–68, 2010.
  * [9] Dan Crisan and Jie Xiong.  Numerical solution for a class of SPDEs over bounded domains.  Stochastics An International Journal of Probability and Stochastic Processes, 86(3):450–472, 2014.
  * [10] Ruth Curtain.  A survey of infinite-dimensional filtering.  SIAM Review, 17(3):395–411, 1975.
  * [11] Ruth F Curtain.  Infinite-dimensional filtering.  SIAM Journal on Control, 13(1):89–104, 1975.
  * [12] Jana de Wiljes and Xin T Tong.  Analysis of a localised nonlinear ensemble Kalman Bucy filter with complete and accurate observations.  Nonlinearity, 33(9):4752, 2020.
  * [13] Pierre Del Moral.  Nonlinear filtering: Interacting particle resolution.  Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 325(6):653–658, 1997.
  * [14] Patrick H Diamond, Sanae-I Itoh, and Kimitaka Itoh.  Modern Plasma Physics: Volume 1, Physical Kinetics of Turbulent Plasmas.  Cambridge University Press, 2010.
  * [15] Arnaud Doucet, Nando De Freitas, Neil James Gordon, et al.  Sequential Monte Carlo methods in practice, volume 1.  Springer, 2001.
  * [16] Paul Dupuis and Richard S Ellis.  A weak convergence approach to the theory of large deviations.  John Wiley & Sons, 2011.
  * [17] Oliver G Ernst, Björn Sprungk, and Hans-Jörg Starkloff.  Analysis of the ensemble and polynomial chaos Kalman filters in Bayesian inverse problems.  SIAM/ASA Journal on Uncertainty Quantification, 3(1):823–851, 2015\.
  * [18] Geir Evensen.  Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.  Journal of Geophysical Research: Oceans, 99(C5):10143–10162, 1994\.
  * [19] Geir Evensen.  The ensemble Kalman filter for combined state and parameter estimation.  IEEE Control Systems Magazine, 29(3):83–104, 2009.
  * [20] Peter L Falb.  Infinite-dimensional filtering: The Kalman-Bucy filter in Hilbert space.  1967\.
  * [21] Uriel Frisch.  Turbulence: the legacy of AN Kolmogorov.  Cambridge university press, 1995.
  * [22] Yuan Gao and Di Qi.  Mean field games for controlling coherent structures in nonlinear fluid systems.  arXiv preprint arXiv:2401.10356, 2024.
  * [23] Michael Ghil, S Cohn, John Tavantzis, K Bube, and Eugene Isaacson.  Applications of estimation theory to numerical weather prediction.  Dynamic meteorology: Data assimilation methods, pages 139–224, 1981\.
  * [24] Georg A Gottwald and Andrew J Majda.  A mechanism for catastrophic filter divergence in data assimilation for sparse observation networks.  Nonlinear Processes in Geophysics, 20(5):705–712, 2013.
  * [25] Carl Graham, Thomas G Kurtz, Sylvie Méléard, Philip E Protter, Mario Pulvirenti, Denis Talay, and Sylvie Méléard.  Asymptotic behaviour of some interacting particle systems; Mckean-Vlasov and Boltzmann models.  Probabilistic Models for Nonlinear Partial Differential Equations: Lectures given at the 1st Session of the Centro Internazionale Matematico Estivo (CIME) held in Montecatini Terme, Italy, May 22–30, 1995, pages 42–95, 1996.
  * [26] Peter L Houtekamer and Herschel L Mitchell.  Data assimilation using an ensemble Kalman filter technique.  Monthly Weather Review, 126(3):796–811, 1998.
  * [27] Rentian Hu, Thomas K Edwards, Leslie M Smith, and Samuel N Stechmann.  Initial investigations of precipitating quasi-geostrophic turbulence with phase changes.  Research in the Mathematical Sciences, 8:1–25, 2021.
  * [28] Rudolph E Kalman and Richard S Bucy.  New results in linear filtering and prediction theory.  1961\.
  * [29] Eugenia Kalnay.  Atmospheric modeling, data assimilation and predictability.  Cambridge university press, 2003.
  * [30] S Kullback and RA Leibler.  The annals of mathematical statistics.  On information and sufficiency, 22, 1951.
  * [31] Thomas G Kurtz and Jie Xiong.  Particle representations for a class of nonlinear SPDEs.  Stochastic Processes and their Applications, 83(1):103–126, 1999\.
  * [32] Thomas G Kurtz and Jie Xiong.  Numerical solutions for a class of SPDEs with application to filtering.  Stochastics in Finite and Infinite Dimensions: In Honor of Gopinath Kallianpur, pages 233–258, 2001.
  * [33] Kody Law, Andrew Stuart, and Kostas Zygalakis.  Data assimilation.  Cham, Switzerland: Springer, 214:52, 2015.
  * [34] Marcel Lesieur.  Turbulence in fluids: stochastic and numerical modelling, volume 488.  Nijhoff Boston, MA, 1987.
  * [35] Martin Leutbecher and Tim N Palmer.  Ensemble forecasting.  Journal of computational physics, 227(7):3515–3539, 2008.
  * [36] Robert S Liptser and Albert N Shiryaev.  Statistics of random processes II: Applications, volume 6.  Springer Science & Business Media, 2013.
  * [37] Andrew J Majda.  Introduction to turbulent dynamical systems in complex systems.  Springer, 2016.
  * [38] Andrew J Majda and Di Qi.  Strategies for reduced-order models for predicting the statistical responses and uncertainty quantification in complex turbulent dynamical systems.  SIAM Review, 60(3):491–549, 2018.
  * [39] Andrew J Majda and Di Qi.  Linear and nonlinear statistical response theories with prototype applications to sensitivity analysis and statistical control of complex turbulent dynamical systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(10), 2019.
  * [40] Andrew J Majda and Xin T Tong.  Performance of ensemble Kalman filters in large dimensions.  Communications on Pure and Applied Mathematics, 71(5):892–937, 2018\.
  * [41] Henry P McKean.  Propagation of chaos for a class of non-linear parabolic equations.  Stochastic Differential Equations (Lecture Series in Differential Equations, Session 7, Catholic Univ., 1967), pages 41–57, 1967\.
  * [42] Sergey Nazarenko and Sergei Lukaschuk.  Wave turbulence on water surface.  Annual Review of Condensed Matter Physics, 7:61–88, 2016.
  * [43] Karl Oelschlager.  A martingale approach to the law of large numbers for weakly interacting stochastic processes.  The Annals of Probability, pages 458–479, 1984.
  * [44] Étienne Pardoux.  Stochastic partial differential equations.  Fudan lecture notes, 87, 2007.
  * [45] Sahani Pathiraja, Sebastian Reich, and Wilhelm Stannat.  McKean–Vlasov SDEs in nonlinear filtering.  SIAM Journal on Control and Optimization, 59(6):4188–4215, 2021\.
  * [46] Francis J Pinski, Gideon Simpson, Andrew M Stuart, and Hendrik Weber.  Kullback–Leibler approximation for probability measures on infinite dimensional spaces.  SIAM Journal on Mathematical Analysis, 47(6):4091–4122, 2015.
  * [47] Di Qi and Jian-Guo Liu.  High-order moment closure models with random batch method for efficient computation of multiscale turbulent systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(10), 2023.
  * [48] Di Qi and Jian-Guo Liu.  A random batch method for efficient ensemble forecasts of multiscale turbulent systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(2), 2023.
  * [49] Di Qi and Andrew J Majda.  Rigorous statistical bounds in uncertainty quantification for one-layer turbulent geophysical flows.  Journal of Nonlinear Science, 28:1709–1761, 2018.
  * [50] Sebastian Reich and Colin Cotter.  Probabilistic forecasting and Bayesian data assimilation.  Cambridge University Press, 2015.
  * [51] Rick Salmon.  Lectures on geophysical fluid dynamics.  Oxford University Press, USA, 1998.
  * [52] Daniel Sanz-Alonso, Andrew Stuart, and Armeen Taeb.  Inverse problems and data assimilation, volume 107.  Cambridge University Press, 2023.
  * [53] Chris Snyder, Thomas Bengtsson, Peter Bickel, and Jeff Anderson.  Obstacles to high-dimensional particle filtering.  Monthly Weather Review, 136(12):4629–4640, 2008.
  * [54] Simone Carlo Surace, Anna Kutschireiter, and Jean-Pascal Pfister.  How to avoid the curse of dimensionality: Scalability of particle filters with and without importance weights.  SIAM review, 61(1):79–91, 2019.
  * [55] Shanyin Tong, Eric Vanden-Eijnden, and Georg Stadler.  Extreme event probability estimation using PDE-constrained optimization and large deviation theory, with application to tsunamis.  Communications in Applied Mathematics and Computational Science, 16(2):181–225, 2021.
  * [56] Jie Xiong.  An introduction to stochastic filtering theory, volume 18.  OUP Oxford, 2008.
  * [57] Tao Yang, Henk AP Blom, and Prashant G Mehta.  The continuous-discrete time feedback particle filter.  In 2014 American control conference, pages 648–653. IEEE, 2014\.
  * [58] Shing-Tung Yau and Stephen S-T Yau.  Finite dimensional filters with nonlinear drift, xi: Explicit solution of the generalized Kolmogorov equation in Brockett–Mitter program.  Advances in Mathematics, 140(2):156–189, 1998.


# Towards a phylogenetic measure to quantify HIV incidence

Pieter Libin equal contribution Artificial Intelligence Lab, Department of
computer science, Vrije Universiteit Brussel, Brussels, Belgium  Department of
Microbiology and Immunology, Rega Institute for Medical Research, KU Leuven -
University of Leuven, Leuven, Belgium  Nassim Versbraegen* Machine Learning
group, Université Libre de Bruxelles, Boulevard du Triomphe CP212, 1050
Bruxelles, Belgium  Interuniversity Institute for Bioinformatics in Brussels,
ULB-VUB, 1050 Brussels, Belgium  Ana B. Abecasis  Department of Microbiology
and Immunology, Rega Institute for Medical Research, KU Leuven - University of
Leuven, Leuven, Belgium  Global Health and Tropical Medicine, GHTM, Instituto
de Higiene e Medicina Tropical, IHMT, Universidade Nova de Lisboa, UNL,
Lisboa, Portugal  Perpetua Gomes  Laboratorio Biologia Molecular, LMCBM, SPC,
HEM, Centro Hospitalar Lisboa Ocidental  Tom Lenaerts  Artificial Intelligence
Lab, Department of computer science, Vrije Universiteit Brussel, Brussels,
Belgium  Machine Learning group, Université Libre de Bruxelles, Boulevard du
Triomphe CP212, 1050 Bruxelles, Belgium  Ann Nowé  Artificial Intelligence
Lab, Department of computer science, Vrije Universiteit Brussel, Brussels,
Belgium

###### Abstract

One of the cornerstones in combating the HIV pandemic is being able to assess
the current state and evolution of local HIV epidemics. This remains a complex
problem, as many HIV infected individuals remain unaware of their infection
status, leading to parts of HIV epidemics being undiagnosed and under-
reported. To that end, we firstly present a method to learn epidemiological
parameters from phylogenetic trees, using approximate Bayesian computation
(ABC). The epidemiological parameters learned as a result of applying ABC are
subsequently used in epidemiological models that aim to simulate a specific
epidemic. Secondly, we continue by describing the development of a tree
statistic, rooted in coalescent theory, which we use to relate epidemiological
parameters to a phylogenetic tree, by using the simulated epidemics. We show
that the presented tree statistic enables differentiation of epidemiological
parameters, while only relying on phylogenetic trees, thus enabling the
construction of new methods to ascertain the epidemiological state of an HIV
epidemic. By using genetic data to infer epidemic sizes, we expect to enhance
understanding of the portions of the infected population in which diagnosis
rates are low.

##  1 Introduction

About 37 million people are currently infected with HIV and an estimated 35
million people have died due to the effects of AIDS (the eventual result of
HIV infection) since the beginning of the epidemic at the start of the
twentieth century [1]. Global efforts have ensued to enhance the collection,
dissemination and accessibility of epidemiological data related to HIV. One of
the most burdensome aspects in curtailing the spread of HIV emerges from
infected individuals being unaware of their infection status. This stems from
the fact that a host can be infected for many years before noticing any
symptoms [2, 3, 4]. As a result, a significant fraction of the HIV infected
population remains undiagnosed, hampering effectiveness of interventions and
assessment of further developments of the epidemic. Consequently, methods that
deliver a well-founded estimate of the number of HIV infected individuals are
paramount [5]. Such an estimate enables deduction of the number of undiagnosed
infected individuals. State-of-the-art methods that aim to provide estimates
of the size of HIV epidemics generally consist of applying compartment models
to routine surveillance data to estimate the number of infected individuals
(i.a. number of new diagnoses over time and CD4+ cell counts) [6, 7].

An abundance of clinical data is available in the context of HIV epidemics, as
upon diagnosis a number of tests are performed and the results thereof
collected. One of those tests comprises of assessing the specific genotype of
the virus infecting a patient [8]. To that extent, the genetic sequence of the
virus is determined. As a result, a vast amount of HIV sequences have been
collected over the last decades.

The main benefit of developing a method to quantify a HIV epidemic that relies
on genetic data is to gain insight into the specific sub-populations that
contain a high rate of undiagnosed individuals, thus allowing for more
effective health policies, through diagnosis strategies that are directed
towards these particular sub-populations.
We validate our research on the HIV-1 epidemic in Portugal (see Section 3.1).
We therefore first present inference of the epidemiological parameters of said
epidemic by applying approximate Bayesian computation (ABC) [9]. We apply ABC
to fit a model that contains the epidemiological parameters in question (see
Section 3.2). We further show that calibrating simulations to specific
epidemics is essential, as the epidemiological dynamics has an important
impact on the shape of the phylogenetic tree (see Section 4). We then
construct a tree statistic that enables differentiation of epidemiological
parameters based on phylogenetic trees (see Section 3.4) and evaluate it on a
set of epidemiological simulations (see Section 4.3).

##  2 Background

###  2.1 Phylogenetic trees

Phylogenetic trees represent evolutionary relationships between organisms. A
rooted phylogenetic tree consists of a root, internal nodes, leaves and
branches interconnecting nodes with other nodes and leaves. The branches of a
phylogenetic tree indicate a measure of distance between the organisms
represented by their respective leaves. This distance can be based on the
amount of genetic change or can represent natural time, by using a molecular
clock [10]. A smaller path between two nodes (i.e. traversing the tree through
the nodes from one leaf to another) thus suggests a stronger evolutionary
relatedness [11].

###  2.2 Epidemiological models

Compartment models are one of the most popular concepts stemming from
mathematical epidemiology [12], Compartment models aim to capture population
dynamics by stratifying individuals into different compartments, through some
difference in state. They include transitions between compartments over time,
representing changes in state in the modelled population. We illustrate this
with the SIR Model [13], which consists of three coupled non-linear ordinary
differential equations, representing the change in each compartment over time;

$\dot{S}=-{\beta}SI$ (1) $\dot{I}={\beta}SI-{\gamma}I$ (2) $\dot{R}={\gamma}I$
(3)

With $\beta$ the infection rate, $\gamma$ the recovery rate and $t$ time.

Figure 1: SIR compartment model, compartments Susceptible ($S$), Infected
($I$) and Recovered ($R$) represented by circles, transitions between them by
arrows, denoted by their associated rates ($\beta$ and $\gamma$)

The model consists of 3 compartments; Susceptible ($S$), Infected ($I$) and
Recovered ($R$) (or removed) and transitions between those compartments (see
Figure 1). $R$, $S$ and $I$ represent the number of people in the respective
compartment at a certain point in time. The rates $\beta$ and $\gamma$ depend
on different factors (i.a. population density, infectiousness and pathology).

The model and extended versions thereof are especially relevant when trying to
answer health policy questions. In our work, such models will be used to
simulate specific HIV epidemics.

###  2.3 Detecting missed infections in phylogenetic trees.

The starting point for our tree statistic was the method presented in [14],
which uses the coalescent [15] to provide an indication of the extent to which
samples are missing, or overly present in a specific phylogenetic tree.

Equation 4 specifies how to calculate node probabilities for each node $j$ in
a tree [14].

$\pi^{j}(t,k,N)=\sum^{j}_{i=1}k^{s}_{i}Pr_{i}(C_{lineage})$ (4)

With $i$ the interval index, $t_{i}$ the interval length, $k$ the number of
lineages to coalesce in an interval, $N$ the population size and
$Pr_{i}(C_{lineage})$ the lineage coalescence probability.

$Z^{j}\;\propto\;n\pi^{j}-np^{j}$ (5)

To be able to gauge to what extent the population making up the tree under
investigation is under- or oversampled, one can convert the node probabilities
to z-scores, by using equation 5. The equation quantifies the relation between
number of leaves expected at a node and number of observed leaves for that
node [14].

###  2.4 Portuguese HIV-1 epidemic

The first reports of HIV-1 diagnoses in Portugal go back to 1983. By 2014,
53072 diagnoses had been reported [16]. As in most European HIV epidemics,
subtype B is the most prevalent HIV-1 subtype, followed by subtype G which is
atypical in a European setting [17]. Within the ”men who have sex with men”
risk group, the number of yearly diagnoses show a mild but steady increase,
while modes of transmission have transformed from intravenous drug use being
the main cause of new infections to infections stemming from heterosexual sex
in the period 2000-2014 [16]. This change is associated with lower diagnoses
rates overall (consistent with diagnoses in Europe in general), possibly
illustrating the beneficial results of a health policy implementation that was
particularly effective for intravenous drug-users (e.g. through providing
single-use needles). We apply our research on data stemming from this
epidemic.

##  3 Methods and Materials

###  3.1 Portuguese Data

The data used in the experiments was made available through a HIV-1 resistance
database from Hospital Egas Moniz. Henceforth, we will refer to the used data
as ’Portuguese data’ for the sake of brevity.
Data storage and querying was achieved through the RegaDB system [18]. Said
system allows for complex querying, which was key to assemble all the relevant
data in an efficient manner. We proceeded by querying genetic sequences of
HIV-1 belonging to distinct patients. In doing so, we assembled three genetic
sequence sets, differing in the HIV-1 subtypes they contain. The first set
only holds subtype B sequences (n= 2216) . The second set is made up of
subtype G sequences (n=1961). And a final set (n= 6079) that does not take
into account the specific subtype (and thus also includes other subtypes
besides B and G). HIV-1 sequences were classified using the Rega typing tool
[19, 20]. Each patient only contributed one genetic sequence to a set, if
multiple sequences were associated with a single patient, the sequence
established first was selected.
For each set, an alignment was created using MAFFT [21]. The resulting
alignments were then used to infer phylogenetic trees using maximum likelihood
trough RAxML [22]. RAxML was used with the GTR-$\gamma$ model. A maximum
likelihood tree was constructed and subsequently annotated through bootstraps.
Bootstrapping was halted automatically based on extended majority rule
consensus trees (i.e. autoMRE). In order to preserve the confidential nature
of the employed patient data, tree inference was carried out on local
computers exclusively and only anonymized patient data was used.

###  3.2 Simulation calibration

In order to validate our tree statistic on a real world epidemic, being able
to generate simulation data that was plausible with regard to the real world
epidemic was essential. We thus proceeded by inferring relevant
epidemiological parameters in order to calibrate subsequent simulations. To
that end, we opted to use ABC [9] to learn said parameters, this approach was
inspired by the work presented in [23]. ABC is closely related to Markov chain
Monte Carlo (MCMC), but unlike MCMC, does not require the calculation of exact
likelihoods, which can be intractable for complex models [23, 24]. Learning
the relevant epidemiological parameters in an ABC setting requires the
presence of some distance measure, as an alternative to the exact likelihoods
used in MCMC approaches. Taking into account that the available epidemic data
we want to infer parameters from exists in the form of a phylogenetic tree,
and the possibility of generating new phylogenetic trees through simulation,
we employ a kernel method developed by Poon [25] as a distance measure between
two trees. In concreto, we rely on a specified compartment model (see Figure 2
for the used model) that enables the generation of trees.

By using the aforementioned kernel method, we assume that correspondence in
trees reflect similarities between the model and the epidemic underlying the
observed tree. In order to explore the possible parameters of the specified
model efficiently, ABC is used. In essence, ABC varies parameter values in
order to simulate more data by using the proposed parameter values in a
specified model and aims to minimise the distance between the newly generated
data and the observed data (in this case using the kernel method as a distance
measure) [9, 23]. Table 1 shows the parameters used in our ABC application.

Parameter | Range | $\sigma$ | Initial | Prior
---|---|---|---|---
$N$ | [$10^{3}$ \- $10^{6}$] | $10^{3}$ | $10^{4}$ |  $X=e^{\mu+\sigma Z}$, $\mu=0.5$ and $\sigma=10000$
$\beta_{i}$ | [$10^{-3}$ \- $10$] | $0.05$ | $0.5$ |  $X=e^{\mu+\sigma Z}$, $\mu=1.0$ and $\sigma=0.01$
$\gamma_{i}$ | [$0$ \- $5$] | $0.01$ | $0.1$ |  $X=e^{\mu+\sigma Z}$, $\mu=1.0$ and $\sigma=0.01$
$\mu$ | [$0$ \- $1$] | $0.002$ | $0.02$ |  $X=e^{\mu+\sigma Z}$, $\mu=1.0$ and $\sigma=0.01$
Table 1: Parameters used in ABC kernel method SI model, $X=e^{\mu+\sigma Z}$
represents the log-normal distribution, $\beta_{i}$ being the infection rate,
$\gamma_{i}$ the mortality rate, $\mu$ the mortality rate from natural causes,
$N$ the population size.

###  3.3 Simulations

rcolgem [26, 27] was used to simulate epidemics, based on parameter ranges
obtained from application of the ABC-kernel method. Each simulation set
consisted of 1000 simulations, outputting phylogenetic trees, and a log of the
population dynamics over time. The used model (based on [28, 29]) is given by
the following equations, and is illustrated in Fig. 2.

$\dot{S}=bN-\mu S-(\beta_{0}I_{0}+\beta_{1}I_{1}+\beta_{2}I_{2})\frac{S}{N}$
$\dot{I_{0}}=(\beta_{0}I_{0}+\beta_{1}I_{1}+\beta_{2}I_{2})\frac{S}{N}-(\mu+\gamma_{0})I_{0}$
$\dot{I_{1}}=\gamma_{0}I_{0}-(\mu+\gamma_{1})I_{1}$
$\dot{I_{2}}=\gamma_{1}I_{1}-(\mu+\gamma_{2})I_{2}$

The model is an extension of a SIR model, and consists of three infection
stages and includes births and deaths (i.e. conceptual addition and removal of
simulated individuals) without a recovery state. Used parameter values,
determined from the results of ABC application, are as follows;
$\gamma_{0}=0.045$, $\gamma_{1}=0.14$, $\gamma_{2}=0.5$ $\mu=0.001$,
$\beta_{0}=0.12$, $\beta_{1}=0.03$ and $\beta_{2}=0.009$ while $S_{0}$ is
varied between $145000$ and $157000$ and the sample size (i.e. the number of
leaves in the tree) between $1000$ and $12000$. Parameters were sampled
between the specified ranges using Latin hypercube sampling. The goal being to
find parameters for the simulation engine that result in phylogenetic trees
that are similar to the ones inferred from the Portugal data, in order to
match the underlying epidemic.

Figure 2: extended SIR compartment model, compartments Susceptible ($S$),
three Infected stages ($I_{i}$) and Deceased ($\theta$) represented by
circles, transitions between them by arrows, denoted by their associated rates
(infections $\beta_{i}$, death through infection $\gamma_{i}$, death trough
natural causes $\mu$ and births $b$)

###  3.4 Constructing the $T_{z}$-score

The starting point in constructing our tree statistic was the method described
in [14]. The result of applying said method on a tree is an annotated tree,
which includes node probabilities for each node in the tree. In order to
construct our tree statistic from such an annotated tree, we devised a
procedure to infer information about the population on the basis of z-scores.
We proceeded by defining a statistic that demonstrates the overall extent to
which a tree is over- or undersampled. We call this statistic $T_{z}$. To
convert the obtained annotated tree to a single statistic value, we rely on
the z-score in the root of the tree, as node probabilities are by definition
propagated to the root node; $T_{z}=\frac{Z^{r}}{s}$, with $Z^{r}$ the z-score
of the root and $s$ the number of samples (i.e. leaves) making up the tree. In
order to apply the method, a tree and a $N$ need to be specified. Through
experimental analyses, we found that a large $N$ is necessary to obtain
informative results, we thus specified $N=10^{5}$.

##  4 Results

###  4.1 Approximate Bayesian computation

Figure 3: ABC kernel method result, $\mu$ over different iterations

Figure 3 presents the results of ABC application for the $\mu$ parameter. The
figure indicates that the ABC chain converged after about 1500 iterations and
had thus learned plausible values for said parameters with regard to the
specified model.

###  4.2 Phylogenetic tree assessment

We first present a visual comparison of tree topologies between the tree
inferred from the Portuguese dataset, a tree obtained through an ABC
application calibrated simulation and a PANGEA [30] simulation, which aims to
model the HIV epidemic in sub-Saharan Africa. If topologies show major
discrepancies, we assume this indicated the simulations are not well
calibrated with regard to the actual epidemic, while a relatively
corresponding topology would indicate simulations resembling the actual
epidemic. The PANGEA tree serves as an example of topology difference when
comparing different epidemics.

(a) Cladogram of tree inferred from Portuguese dataset

(b) Histogram of patristic distances in tree inferred from Portuguese dataset

(c) Cladogram of tree obtained through rcolgem simulation with simulation
parameters based on ABC results

(d) Histogram of patristic distances in tree obtained through rcolgem
simulation

(e) Cladogram of tree obtained through PANGEA simulation

(f) Histogram of patristic distances in tree obtained through PANGEA
simulation

Figure 4: Comparison of topology of inferred and simulated phylogenetic trees

Figures 4(a), 4(c) and 4(e) provide a visual representation of the relevant
phylogenetic trees. Figure 4(a) stems from the Portuguese dataset, and as such
offers a baseline of desired tree topology. Figure 4(c) presents the tree
obtained through a ABC calibrated simulation, and 4(e) shows a tree from a
PANGEA simulation. We demonstrate visually that the relevant tree topologies
display a relatively high level of correspondence. We further investigated
tree correspondence by using patristic distances [31]. The patristic distance
between two leaves $l_{0}$ and $l_{1}$ is the number of changes needed to
$l_{0}$ in order for it to become identical to $l_{1}$ [32]. Figures 4(b),
4(b) and 4(f) present a histogram of the patristic distances present in the
tree stemming from the Portuguese dataset the ABC calibrated simulated tree
and the PANGEA tree respectively. These show that the ABC calibrated
simulation tree is concordant with the tree stemming from the Portuguese
epidemic.

###  4.3 $T_{z}$-score distribution

Figure 5: rcolgem simulation analysis, sample size (i.e. number of leaves)
against $T_{z}$-score, with number of infected individuals in the simulations,
shown as a gradient from $2.56\times 10^{4}$ (red) to $4.64\times 10^{4}$
(blue), $N=10^{5}$

Figure 5 presents the result of application of our method on 1000 trees
obtained from rcolgem simulations. In the figure, obtained $T_{z}$ scores are
plotted against tree sample sizes and the number of infected individuals in
the simulations, shown as a gradient from $2.56\times 10^{4}$ (red) to
$4.64\times 10^{4}$ (blue). As a reminder, the number of infected individuals
is determined by $\beta_{i}$ and $S_{0}$. We can clearly observe a
distribution that allows distinction of the number of infected by $T_{z}$
scores and sample sizes. Indicating that we obtained an informative
distribution through application of our method. The figure shows that a lower
$T_{z}$ score correlates with a larger portion of the infected population not
being included in the phylogenetic tree. Additionally, as sample size goes
down the distribution becomes wider, indicating sufficiently large trees are
necessary to allow for meaningful inferences.

##  5 Discussion and Conclusions

The presented results show that the presented $T_{z}$ score enables
differentiation of epidemiological parameters based on phylogenetic trees. As
such, an interesting further development would be to learn a function
$f(S_{0},\beta,\gamma,N)\to T^{N}_{z}$, i.e. learning the effect of
epidemiological parameters on $T_{z}$ values in general, and from there,
$f(T^{N}_{z})\to\\{\beta,\gamma\\}$, i.e. constructing a function that relates
obtained $T_{z}$ scores to possible epidemiological parameters, which we
assume to be feasible, as $S_{0}$ should be ascertainable trough census data
and epidemiological studies.
We would further like to extend the $T_{z}$ score to include a measure of
uncertainty associated to inferences. A possible approach to accomplish this
would be to apply our method on multiple subtrees, resulting from pruning the
tree under investigation, and determining the extent to which results remain
coherent with regard to the number of pruned leaves. Additionally,
investigating the effect of simulating epidemics using current models
specifically tailored to HIV-1 (e.g. an approach where the currently prevalent
CD4+ models would be adapted to generate phylogenetic trees) would be an
interesting further development. Planned further research includes adaptation
of the method presented in [14] to draw coalescent probabilities from a
distribution that is specific to HIV evolutionary dynamics.

We have presented a tree statistic that can be employed to assort phylogenetic
trees on the basis of their underlying epidemiological parameters. By doing
so, we provide a first step towards a method to infer epidemiological
parameters from phylogenetic trees using coalescent theory, which would
additionally be able to indicate the specific subpopulations in which
diagnosis rates are low, providing a crucial tool for health policy
researchers.

## Acknowledgments

Pieter Libin was supported by a PhD grant of the FWO (Fonds Wetenschappelijk
Onderzoek Vlaanderen) and a grant of the VUB research council (VUB/OZR2714).

## References

  * [1] UNAIDS.  Fact sheet - latest global and regional statistics on the status of the aids epidemic., July 2017.
  * [2] Lisa B Hightow-Weidman, Carol E Golin, Kelly Green, Eva NP Shaw, Pia DM MacDonald, and Peter A Leone.  Identifying people with acute hiv infection: demographic features, risk factors, and use of health care among individuals with ahi in North Carolina.  AIDS and Behavior, 13(6):1075–1083, 2009.
  * [3] Gregory F Burton, Brandon F Keele, Jacob D Estes, Tyler C Thacker, and Suzanne Gartner.  Follicular dendritic cell contributions to hiv pathogenesis.  In Seminars in immunology, volume 14, pages 275–284. Elsevier, 2002\.
  * [4] Kristof Theys, Pieter Libin, Andrea-Clemencia Pineda-Pena, Ann Nowe, Anne-Mieke Vandamme, and Ana B Abecasis.  The impact of hiv-1 within-host evolution on transmission dynamics.  Current opinion in virology, 28:92–101, 2018.
  * [5] Elizabeth M Begier, Yussef Bennani, Lisa Forgione, Amado Punsalang, David B Hanna, Jeffrey Herrera, Lucia Torian, Maria Gbur, Kent A Sepkowitz, and Farah Parvez.  Undiagnosed hiv infection among New York City jail entrants, 2006: results of a blinded serosurvey.  JAIDS Journal of Acquired Immune Deficiency Syndromes, 54(1):93–101, 2010.
  * [6] Ard van Sighem, Fumiyo Nakagawa, Daniela De Angelis, Chantal Quinten, Daniela Bezemer, Eline Op de Coul, Matthias Egger, Frank de Wolf, Christophe Fraser, and Andrew Phillips.  Estimating hiv incidence, time to diagnosis, and the undiagnosed hiv epidemic using routine surveillance data.  Epidemiology (Cambridge, Mass.), 26(5):653, 2015.
  * [7] Alessia Mammone, Patrizio Pezzotti, Vincenza Regine, Laura Camoni, Vincenzo Puro, Giuseppe Ippolito, Barbara Suligoi, and Enrico Girardi.  How many people are living with undiagnosed hiv infection? an estimate for Italy, based on surveillance data.  AIDS (London, England), 30(7):1131, 2016.
  * [8] Madhu Vajpayee and Teena Mohan.  Current practices in laboratory monitoring of hiv infection.  The Indian journal of medical research, 134(6):801, 2011.
  * [9] Mikael Sunnåker, Alberto Giovanni Busetto, Elina Numminen, Jukka Corander, Matthieu Foll, and Christophe Dessimoz.  Approximate bayesian computation.  PLoS computational biology, 9(1):e1002803, 2013.
  * [10] Mark Ridley.  Evolution.  John Wiley & Sons Incorporated, 2009.
  * [11] Marco Salemi, Philippe Lemey, and Anne-Mieke Vandamme.  The phylogenetic handbook: a practical approach to phylogenetic analysis and hypothesis testing.  Cambridge University Press, 2009.
  * [12] Fred Brauer.  Compartmental models in epidemiology.  In Mathematical epidemiology, pages 19–79. Springer, 2008.
  * [13] William O Kermack and Anderson G McKendrick.  A contribution to the mathematical theory of epidemics.  In Proceedings of the Royal Society of London A: mathematical, physical and engineering sciences, volume 115, pages 700–721. The Royal Society, 1927.
  * [14] Stephanie Stacy, Allison Black, and Trevor Bedford.  Using the coalescent framework to detect missed infections in phylogenetic trees.  SURP, Fred Hutchinson, 2016.
  * [15] John FC Kingman.  The coalescent.  Stochastic processes and their applications, 13(3):235–248, 1982\.
  * [16] António Diniz, Raquel Duarte, Joana Bettencourt, Teresa Melo, Paulo Nogueira, Carla Farinha, Maria Isabel Alves, José Martins, and Ana Lisette Oliveira.  Portugal infeção vih, sida e tuberculose em números, 2015\.  Portugal Infeção VIH, SIDA e Tuberculose em números, 2015, pages 5–70, 2015.
  * [17] Ana B Abecasis, Annemarie MJ Wensing, Dimitris Paraskevis, Jurgen Vercauteren, Kristof Theys, David AMC Van de Vijver, Jan Albert, Birgitta Asjö, Claudia Balotta, Danail Beshkov, et al.  Hiv-1 subtype distribution and its demographic determinants in newly diagnosed patients in europe suggest highly compartmentalized epidemics.  Retrovirology, 10(1):7, 2013.
  * [18] Pieter Libin, Gertjan Beheydt, Koen Deforche, Stijn Imbrechts, Fossie Ferreira, Kristel Van Laethem, Kristof Theys, Ana Patricia Carvalho, Joana Cavaco-Silva, Giuseppe Lapadula, et al.  Regadb: community-driven data management and analysis for infectious diseases.  Bioinformatics, 29(11):1477–1480, 2013.
  * [19] Andrea-Clemencia Pineda-Peña, Nuno Rodrigues Faria, Stijn Imbrechts, Pieter Libin, Ana Barroso Abecasis, Koen Deforche, Arley Gómez-López, Ricardo J Camacho, Tulio de Oliveira, and Anne-Mieke Vandamme.  Automated subtyping of hiv-1 genetic sequences for clinical and surveillance purposes: performance evaluation of the new rega version 3 and seven other tools.  Infection, Genetics and Evolution, 19:337–348, 2013.
  * [20] Luiz Carlos Junior Alcantara, Sharon Cassol, Pieter Libin, Koen Deforche, Oliver G Pybus, Marc Van Ranst, Bernardo Galvao-Castro, Anne-Mieke Vandamme, and Tulio de Oliveira.  A standardized framework for accurate, high-throughput genotyping of recombinant and non-recombinant viral sequences.  Nucleic acids research, 37(suppl_2):W634–W642, 2009.
  * [21] Kazutaka Katoh, Kazuharu Misawa, Kei-ichi Kuma, and Takashi Miyata.  Mafft: a novel method for rapid multiple sequence alignment based on fast fourier transform.  Nucleic acids research, 30(14):3059–3066, 2002.
  * [22] Alexandros Stamatakis.  Raxml version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies.  Bioinformatics, 30(9):1312–1313, 2014.
  * [23] Art FY Poon.  Phylodynamic inference with kernel abc and its application to hiv epidemiology.  Molecular biology and evolution, 32(9):2483–2495, 2015.
  * [24] Pieter Libin, Laurens Hernalsteen, Kristof Theys, Perpetua Gomes, Ana Abecasis, and Ann Nowe.  Bayesian inference of set-point viral load transmission models.  In 30th Benelux Conference on Artificial Intelligence, pages 107–121, 2018.
  * [25] Art FY Poon, Lorne W Walker, Heather Murray, Rosemary M McCloskey, P Richard Harrigan, and Richard H Liang.  Mapping the shapes of phylogenetic trees from human and zoonotic rna viruses.  PLoS one, 8(11):e78122, 2013.
  * [26] David A Rasmussen, Erik M Volz, and Katia Koelle.  Phylodynamic inference for structured epidemiological models.  PLoS computational biology, 10(4):e1003570, 2014.
  * [27] Erik M Volz.  Complex population dynamics and the coalescent under neutrality.  Genetics, 190(1):187–201, 2012.
  * [28] Erik M Volz.  Estimating hiv transmission rates with rcolgem, 2014.  accessed: August 2017.
  * [29] Erik M Volz.  Simulating genealogies with an epidemiological coalescent model using rcolgem, 2015.  accessed: August 2017.
  * [30] Oliver Ratmann, Emma B. Hodcroft, Michael Pickles, Anne Cori, Matthew Hall, Samantha Lycett, Caroline Colijn, Bethany Dearlove, Xavier Didelot, Simon Frost, A.S. Md Mukarram Hossain, Jeffrey B. Joy, Michelle Kendall, Denise Kühnert, Gabriel E. Leventhal, Richard Liang, Giacomo Plazzotta, Art F.Y. Poon, David A. Rasmussen, Tanja Stadler, Erik Volz, Caroline Weis, Andrew J. Leigh Brown, Christophe Fraser, and .  Phylogenetic tools for generalized hiv-1 epidemics: Findings from the pangea-hiv methods comparison.  Molecular Biology and Evolution, 34(1):185–203, 2017.
  * [31] Pieter Libin, Ewout Vanden Eynden, Francesca Incardona, Ann Nowé, Antonia Bezenchek, EucoHIV Study Group, Anders Sönnerborg, Anne-Mieke Vandamme, Kristof Theys, and Guy Baele.  Phylogeotool: interactively exploring large phylogenies in an epidemiological context.  Bioinformatics, 33(24):3993–3995, 2017.
  * [32] Tod F Stuessy and Christiane König.  Patrocladistic classification.  Taxon, 57(2):594–601, 2008.